Unknown Treasure(Lucas+中国剩余定理）

2018-08-11

Unknown Treasure

题意

给你n，m，和K个素数$p_1,p_2 …… p_k$，求$C_n^mmod\prod{p_i}$。

n,m<=1e18。

思路

计算组合数用Lucas即可，每次取模Pi，然后由于最终是要模$\prod{p_i}$，所以用中国剩余定理。

令 $tot = \prod{p_i} $，$X=C_n^mmod\prod{p_i}$，$K=\lfloor C_n^m/\prod{p_i}\rfloor$

有：$C_n^m=K*tot+X$

两边同对$pi$取模，

$$
X % p_1 = C_n^m % p_1 = Lucas(n,m,p_1)
$$

$$
X % p_2 = C_n^m % p_2 = Lucas(n,m,p_2)
$$

$$
…
$$

$$
X % p_n = C_n^m % p_n = Lucas(n,m,p_n)
$$

$$
X % p_1 = C_n^m % p_1 = Lucas(n,m,p_1)
$$

利用中国剩余定理求出X即可。

因为做乘法过程中可能会超long long ，所以我们要用快速乘法去处理，方便取模

看！代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 100005;

ll fac[N];

ll modpow(ll a, ll b, ll MOD)
{
    ll ret = 1;
    while(b)
    {
        if(b & 1) ret = (ret * a)%MOD;
        a = (a * a)%MOD;
        b >>= 1;
    }
    return ret;
}

ll modmul(ll a, ll b, ll MOD)  //快速乘
{
    ll ret = 0;
    while(b)
    {
        if(b & 1) ret = (ret + a) % MOD;
        a = (a + a) % MOD;
        b >>= 1;
    }
    return ret;
}

ll getFactor(ll p) 
{
    fac[0] = 1;
    for(int i = 1; i <= p; i++)
        fac[i] = (fac[i-1]*i) % p;
}

ll Lucas(ll n, ll m, ll p)
{
    ll ret=1;
    while(n && m)
    {
        ll a = n%p, b = m%p;
        if(a < b) return 0;
        ret = (ret * fac[a] * modpow(fac[b]*fac[a-b]%p, p-2, p)) % p;
        n/=p;
        m/=p;
    }
    return ret;
}

ll exgcd (ll a, ll b, ll &x, ll &y)
{
    if (!b)
    {
        x = 1, y = 0;
        return a;
    }
    int ans = exgcd ( b, a % b, y, x );
    y -= a / b * x;
    return ans;
}

//中国剩余定理，a[i]存放余数，m[i]存放两两互质的数
ll remainder(ll a[], ll m[], int len)
{
    ll d, x, y, ret = 0;
    ll M = 1;
    for (int i = 0; i < len; i++) M *= m[i];
    for (int i = 0; i < len; i++)
    {
        ll w = M / m[i];
        d = exgcd(m[i], w, x, y);
        ret = (ret + modmul(modmul(y, w, M), a[i], M) ) % M; //这里 y*w*a[i]直接乘会超long long，所以要用快速乘，M为每次的MOD。
    }
    return (ret + M) % M;
}

ll n, m;
int K;
ll p[15], lucas[15];

int main()
{
    int t;
    scanf("%d", &t);
    while(T--)
    {
        scanf("%lld%lld%d", &n, &m, &K);
        for(int i = 0; i < K; i++)
        {
            scanf("%lld", &p[i]);
            getFactor(p[i]);
            lucas[i] = Lucas(n, m, p[i]); //对于每一个pi计算一次C(n,m)
        }
        ll ans = remainder(lucas, p, K);
        printf("%lld\n", ans);
    }
    return 0;
}

展开全文 >>

在比赛中遇到的稀奇古怪的结论

2018-08-09

将多边形对角线相连，其对角线最多被分割成多少段？

$$
max=(n^4-6n^3+17n^2-24n)/12
$$

一个圆，一个圆的内接n边形，连接多边形对角线，求整个圆被分割成几块?

$$
max=(n^4-6n^3+23n^2-18n+24)/24
$$

根号函数的性质：
$$
\sqrt{a+b+c+…}\le\sqrt{a}+\sqrt{b}+\sqrt{c}+…
$$
求$(\sum_{i=1}^mx_i )= k (0≤x_i < n)$的整数解的组数(Character Encoding )

分析后发现答案为$(1+x+x^2+x^3+…+x^n)^m$的展开式中项的次数为k的项的系数，对于n,m,k，有答案：
$$
\sum_{i=0}^{\lfloor{\frac{k}{n}}\rfloor}(-1)^iC_m^iC_{m+k-1-n*i}^{m-1}
$$
多边形公式

$$S=\frac{1}{2}\left|\sum\limits_{i=0}^{n-1}(X_iY_{i+1}-X_{i+1}Y_i)\right|$$

其中是相邻的点，对于n-1的i+1是0代表这第一个点，由于顺时针逆时针会使正负号改变，所以取绝对值。

摆线与圆面积关系

摆线留下的面积是圆面积的三倍。
它的长度等于旋转圆直径的 4 倍，是一个不依赖于π的有理数。

摆线定义：摆线(cycloid)是数学中众多的迷人曲线之一．它是这样定义的：一个圆沿一直线缓慢地滚动，则圆上一固定点所经过的轨迹称为摆线．

两条直线夹角

$${\left|\frac{A_{1}B_{2}-A_{2}B_{1}}{A_{1}A_{2}+B_{1}B_{2}}\right|}$$

素数定理
$$\lim\limits_{n->\infty}\frac{\pi(n)}{n/\ln{n}}$$

描述素数的分布，素数分布规律，以36N（N+1）为单位，随着N的增大，素数的个数发波浪形式渐渐增多。素数定理可以给出第n个素数p(n)的渐近估计：它也给出从整数中抽到素数的概率。从不大于n的自然数随机选一个，它是素数的概率大约是1/ln n。 π(n)为小于等于n的素数的个数。

GCD(1)
a> 0；m,n > 0，那么有$$\gcd(a^m-1,a^n-1)=a^{gcd(m,n)}-1$$

GCD(2)
a>b;gcd(a,b)=1,那么$\gcd(a^m-b^m,a^n-b^n)=a^{\gcd(m,n)}-{b^{\gcd(m,n)}} $

GCD(3)
设$$\gcd(a^m-1,a^n-1)=a^{gcd(m,n)}-1$$,那么G的值为
(1)n为素数，那么答案就是n。
(2)n有多个素因子，那么答案就是1。
(3)n只有一个素因子，那么答案就是该因子。

GCD(4)

$$\sum\limits_{i=1}^{N}\gcd(i,N)=\sum\limits_{d|N}dφ(\frac{N}{d})=\sum\limits_{d|N}\frac{N}{d}φ(d)$$

φ(x)为欧拉函数，d是N的质因子

欧拉函数

$\sum_{d|n}\varphi(d)=n $

Fibonacci(1)
设Fn为Fib数，那么有$$\gcd(F_m,F_n) = F_{gcd(m,n)}$$

Fibonacci(2)
给定两个互素的正整数A和B，那么它们最大不能组合的数为A∗B−A−B不能组合的个数为$$num=\frac{(A-1)(B-1)}{2}$$

lcm(1)
$$(n+1)lcm(C_n^0,C_n^1,\cdots,C_n^n)=lcm(1,2,\cdots,n+1)$$

lcm(2)
给一个正整数n，求lcm(1,2,3,⋯,n)lcm(1,2,3,⋯,n)的值，(1<=n<=10^8)

定义L(n)L(n)为1,2,3,……,n的最小公倍数。则可以发现：
$$
L(n+1)=L(n)*p \ \ \ \ \ \ \ \ \ \ 如果(n+1)是p的次方数
$$

$$
= L(n) \ \ \ \ \ \ \ \ \ \ \ \ 否则
$$

所以有:

L(1) = 1
L(2) = 1 x 2
L(3) = 1 x 2 x 3
L(4) = 1 x 2 x 3 x 2
L(5) = 1 x 2 x 3 x 2 x 5
L(6) = 1 x 2 x 3 x 2 x 5
L(7) = 1 x 2 x 3 x 2 x 5 x 7

φ(ab)和gcd(a,b)的关系
$$
\frac{φ(ab)}{φ(a)φ(b)}=\frac{gcd(a,b)}{φ(gcd(a,b))}
$$

展开全文 >>

Harvest of Apples (莫队)

2018-08-09

Harvest of Apples

题意

题意是从n个苹果里面选取最多m个苹果的方案数量。m、n、t<=1e5

思路

莫队分块+组合预处理，原来莫队中的区间[l,r]在这里用n、m代替了。

接下来是推公式部分：

S(n,m)表示$\sum_{i=1}^n\sum_{j=1}^mC_i^j$

我们可以得出$S(n,m)=S(n,m-1)+C_n^m$ ,——-①

又因为有$C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$（组合数学的一个基本公式），因此又可以推出：$S(n,m)=S(n-1,m)+S(n-1,m-1)=2S(n-1,m)-C_{n-1}^m$——–②

有了①②两个式子，我们就可以得到S(n,m)和S(n-1,m)、S(n+1,m)、S(n,m-1)、S(n,m+1)之间的关系，先预处理出n！和相应的需要的逆元，然后把输入的查询分块，按照代码中所示的“套路”排序法，排一排，然后对于每个n，m开始处理，分四种情况这样。

看！代码

#include <bits/stdc++.h>
using namespace std;
const int N = 100010;
const long long MOD =1e9+7;
typedef long long ll;
ll multi[N];
ll devide[N];
ll ans[N];
ll pos[N];
ll er;
int D=(int)sqrt(1.0*N);
struct node{
    int n,m;
    int id;
    bool operator <(const node& p) const{
        if(pos[n]==pos[p.n]) return m<p.m;
        return n<p.n;                           //N大的排前面
    }
};
node qujian[N];
ll ex_gcd(ll a,ll b,ll &x,ll &y)//扩展欧几里得（扩展gcd）
{
	if (a==0&&b==0) return -1;
	if (b==0){x=1;y=0;return a;}
	ll d=ex_gcd(b,a%b,y,x);
	y-=a/b*x;
	return d;
}

ll mod_inverse(ll a,ll n)//乘法逆元,n取模
{
	ll x,y;
	ll d = ex_gcd(a,n,x,y);
	return (x%n+n)%n;
}

void init(int n)
{
    multi[0]=1;
    multi[1]=1;
    for(int i=2;i<=n;i++)
        multi[i]=multi[i-1]*i%MOD;
    devide[n]=mod_inverse(multi[n],MOD);
    for(int i=n-1;i>=0;i--)
        devide[i]=devide[i+1]*(i+1)%MOD;
}
ll cal(int n,int m)
{
    return ((multi[n]*devide[m])%MOD*devide[n-m])%MOD;
}
void operation1(int n,int m,ll &res)
{
    res=((2ll*res)%MOD-cal(n,m)%MOD+MOD)%MOD;
}
void operation2(int n,int m,ll &res)
{
    res=(res+cal(n-1,m))%MOD*er%MOD;
}
void operation3(int n,int m,ll &res)
{
    res=(res+cal(n,m+1))%MOD;
}
void operation4(int n,int m,ll &res)
{
    res=(res-cal(n,m)+MOD)%MOD;
}
int main()
{
    init(100000);
    int t;
    int n,m;
    er=mod_inverse(2ll,MOD);
    scanf("%d",&t);
    for(int i=1;i<=t;i++)
    {
        scanf("%d%d",&qujian[i].n,&qujian[i].m);
        qujian[i].id=i;
        pos[i]=i/D;
    }
    sort(qujian+1,qujian+t+1);
    int curl,curr;
    curl=curr=1;
    ll res=2;
    for(int j=1;j<=t;j++)
    {
        while(curl<qujian[j].n) operation1(curl++,curr,res);
        while(curr<qujian[j].m) operation3(curl,curr++,res);
        while(curl>qujian[j].n) operation2(curl--,curr,res);
        while(curr>qujian[j].m) operation4(curl,curr--,res);
        ans[qujian[j].id]=res;
    }
    for(int i=1;i<=t;i++)
        printf("%lld\n",ans[i]);
    return 0;
}

展开全文 >>

Glad You Came（线段树||ST表）

2018-08-09

Glad You Came

题意

不晓得咋说……反正就是把m个区间对应m个vi，对于每一个区间的数字，小于vi的都变成vi，最后求最终结果的i*a[i]的异或。

1、线段树：肯定要维护一个最小值，而且这题由于初始的ai都是0，因此min也是0，而且暴力单点更新，我们只需要建一个一维的数组就好了。然后一个update函数，一个最终结果的计算函数，注意剪枝：” $if(tree[index].val>=v) return “$。

2、反向ST表：我的理解是，对于一个区间$dp[i][d]$，表示的是从i位置开始长度为$2^d$的区间,也就是$[i , i+2^d-1]$。

比如对于区间[1,10]，长度为10，那么我们找到最大的那个<10的$2^d$，即8，d=3，可将区间分成两个：[1,8],[3,10]，用dp[],[]表示就是$dp[1][3]$和$dp[10-2^3+1][3]$。这样对于每一次询问L,R，从一个大区间分出来的两个小区间，维护他们的max值，取max(dp ,val ) 。这样，我们就从最大的区间开始往小了维护，直到$dp[i][0]$这样。

看！代码

/*线段树莽过的做法*/
#include <bits/stdc++.h>
#define MAXN 100010
#define inf 0x3f3f3f3f
using namespace std;
typedef unsigned int uint;

long long tree[MAXN<<2];//一定要开到4倍多的空间
int cnt;
long long ans;

uint mod=1<<30;
uint fx(uint &x,uint &y,uint &z)
{
    x=x^(x<<11);
    x=x^(x>>4);
    x=x^(x<<5);
    x=x^(x>>14);
    uint w=x^(y^z);
    x=y;
    y=z;
    z=w;
    return z;
}

int L,R;
long long v;
void update(int l,int r,int index)
{
    if(tree[index]>=v)
        return;
    if(l==r&&tree[index]<v)
    {
        tree[index]=v;
        return ;
    }
    int mid=l+r>>1;
    if(L<=mid)update(l,mid,index<<1);
    if(R>mid)update(mid+1,r,index<<1|1);

    tree[index]=min(tree[index<<1],tree[index<<1|1]);
}
long long cal(int l,int r,int index)
{
    if(l==r)
        return tree[index];
    int mid=l+r>>1;
    if(cnt<=mid)
        return cal(l,mid,index<<1);
    else
        return cal(mid+1,r,index<<1|1);

}
int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        memset(tree,0,sizeof(tree));
        int n,m;
        cnt=1;
        long long ans=0;
        scanf("%d%d",&n,&m);
        uint x,y,z;
        scanf("%u%u%u",&x,&y,&z);
        for(int i=1;i<=m;i++)
        {
            L=fx(x,y,z)%n+1;
            R=fx(x,y,z)%n+1;
            v=fx(x,y,z)%mod;
            if(L>R) swap(L,R);
            update(1,n,1);
        }
        for(cnt=1;cnt<=n;cnt++)
            ans=ans^(cnt*cal(1,n,1));
        printf("%lld\n",ans);
    }
    return 0;
}

#include <bits/stdc++.h>
#define MAXN 100010
#define inf 0x3f3f3f3f
using namespace std;
typedef unsigned int uint;
uint f[15000010];
uint mod=1<<30;
uint fx(uint &x,uint &y,uint &z)
{
    x=x^(x<<11);
    x=x^(x>>4);
    x=x^(x<<5);
    x=x^(x>>14);
    uint w=x^(y^z);
    x=y;
    y=z;
    z=w;
    return z;
}

int l[5000010],r[5000010];
long long v[5000010];
long long dp[100010][18];
int lg[100010];
int main()
{
    int t;
    for(int i=1;i<=100000;i++) lg[i]=log2(i);
    scanf("%d",&t);
    while(t--)
    {
        int n,m;
        scanf("%d%d",&n,&m);
        uint x,y,z;
        scanf("%u%u%u",&x,&y,&z);
        for(int i=1;i<=3*m;i++)
            f[i]=fx(x,y,z);
        for(int i=1;i<=m;i++)
        {
            l[i]=min(f[i*3-2]%n+1,f[i*3-1]%n+1);
            r[i]=max(f[i*3-2]%n+1,f[i*3-1]%n+1);
            v[i]=f[i*3]%mod;
        }
        memset(dp,0,sizeof(dp));
        int d,MAX=-1;
        for(int i=1;i<=m;i++)
        {
            d=lg[r[i]-l[i]+1];
            MAX=max(MAX,d);//记录最长的长度
            dp[l[i]][d]=max(v[i],dp[l[i]][d]); //分为两部分，每一部分取覆盖的最小值和v比较
            dp[r[i]-(1<<d)+1][d]=max(v[i],dp[r[i]-(1<<d)+1][d]);
        }
        for(int j=MAX;j>0;j--)
        {
            for(int i=1;i+(1<<j)-1<=n;i++)
            {
                dp[i][j-1]=max(dp[i][j],dp[i][j-1]);
                dp[i+(1<<(j-1))][j-1]=max(dp[i][j],dp[i+(1<<(j-1))][j-1]);
            }
        }
        long long ans=0;
        for(int i=1;i<=n;i++)
            ans=ans^(1ll*i*dp[i][0]);
        printf("%lld\n",ans);
    }
    return 0;
}

展开全文 >>

莫比乌斯题集（有丶意思）

2018-08-09

[YY的GCD](https://www.luogu.org/problemnew/show/P2257)

题意

给你一个n,m，求gcd(i,j)=质数的(x,y)有多少对？(1<=x<=N,1<=y<=M)

思路

题目实际上是要求：
$$
ans=\sum_{i=1}^N\sum_{j=1}^Mgcd(i,j)=prim
$$
由莫比乌斯反演的基本公式：
$$
F(n)=\sum_{n|d}f(d)
$$

$$
f(n)=\sum_{n|d}μ(\frac{d}{n})F(n)
$$

听说一般这种题都是设f(d)= gcd(i,j)==d的个数，F(n) =gcd(i,j)== d和d的倍数的个数 (那么对于上界N和M，就有$F(n)=\lfloor\frac{N}{n}\rfloor \lfloor\frac{M}{n}\rfloor$，因为对于(i,j)，它们的上限为(N,M)，gcd(i,j)=n,而F(n)求的是d和d的倍数的个数，所以对于i，j来说，这样的搭配有N/n*M/n个）。

所以我们可以得到下面的式子：
$$
f(d)=\sum_{i=1}^N\sum_{j=1}^M[gcd(i,j)==d]
$$

$$
F(n)=\sum_{n|d}f(d) =\lfloor\frac{N}{n}\rfloor \lfloor\frac{M}{n}\rfloor
$$

带入莫比乌斯公式：

$$
f(n)=\sum_{n|d}μ(\frac{d}{n})F(n)=\lfloor\frac{N}{n}\rfloor \lfloor\frac{M}{n}\rfloor\sum_{n|d}μ(\frac{d}{n})
$$

$$
ans=\sum_{i=1}^{N}\sum_{j=1}^{M}[gcd==p]
$$

$$
=\sum_{p=prim}f(p)
$$

$$
=\sum_{p=prim}\sum_{p|d}μ(\frac{d}{p})\lfloor\frac{N}{p}\rfloor \lfloor\frac{M}{p}\rfloor
$$

枚举$\lfloor\frac{d}{p}\rfloor$,或者说令d=dp，那么d就变成了p的倍数，所以sum枚举的值就变成了1~d(d最大可以取到N/p和M/p中的最小值，因为p是gcd(i,j)==p，那么就不可能取i，j中大的值)：
$$
ans=\sum_{p=prim}\sum_{d=1}^{min(\lfloor\frac{N}{p}\rfloor ,\lfloor\frac{M}{p}\rfloor)}μ(d)\lfloor\frac{N}{dp}\rfloor \lfloor\frac{M}{dp}\rfloor
$$
将dp=T，则p∈{t|T%t==0&&t==prim}：
$$
ans=\sum_{T=1}^{min(N,M)}\sum_{t|T,t∈prim}μ(\lfloor\frac{T}{t}\rfloor)\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor
$$

$$
=\sum_{T=1}^{min(N,M)}\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor\sum_{t|T,t∈prim}μ(\lfloor\frac{T}{t}\rfloor)
$$

就可以做了。

看！代码

第一份O(n)直接求和的代码，当然会T。

#include<bits/stdc++.h>
using namespace std;
const int N=10000010;
int prim[N];//316067
bool vis[N];
int mu[N];
int cnt=0;
//埃式筛法，这里没有用到
//求n以内的质数的个数
void get_list(int n)
{
    cnt=0;
    memset(vis,false,sizeof(vis));
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]) prim[++cnt]=i;
            for(int j=1;j<=cnt&&i*prim[j]<=n;j++){
                   vis[i*prim[j]]=1;
                   //合数标为1，同时，prime[j]是合数i*prime[j]的最小素因子
                   if(i%prim[j]==0) break;
                    //即比一个合数大的质数和该合数的乘积可用一个更大的合数和比其小的质数相乘得到
             }
       }
}
/** 对于μ函数的线筛 */
void get_mu(int n)
{
	mu[1]=1;
	memset(vis,0,sizeof(vis));
	for(int i=2;i<=n;i++)
	{
    	if(!vis[i]){prim[++cnt]=i;mu[i]=-1;}
   	 	for(int j=1;j<=cnt&&prim[j]*i<=n;j++)
        {
        	vis[prim[j]*i]=1;
        	if(i%prim[j]==0)break;
        	else mu[i*prim[j]]=-mu[i];
    	}
	}
}
int main()
{
    get_mu(10000000);
    int t,n,m;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d",&n,&m);
        int Min=min(n,m);
        long long ans=0;
        for(int i=1;i<=Min;i++)
        {
            long long sum=0;
            for(int j=1;j<=cnt&&prim[j]<=i;j++)
            {
                if(i%prim[j]==0)
                    sum+=mu[i/prim[j]];
            }
            sum*=(n/i)*(m/i);
            ans+=sum;
        }
        printf("%lld\n",ans);
    }
    return 0;
}

第二份用了整数分块优化，写的我头真大，比如多加一个memset，T，int开了long long T，

#include<bits/stdc++.h>
using namespace std;
const int N=10000010;

bool vis[N];
int prim[N],pre[N],mu[N];
int cnt;
long long sum[N];

/** 对于μ函数的线筛*/
void get_mu(int n)
{
   // memset(vis,false,sizeof(vis));
    cnt=0;
	mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]){mu[i]=-1;prim[++cnt]=i;}
        for(int j=1;j<=cnt&&prim[j]*i<=n;j++)
        {
            vis[i*prim[j]]=1;
            if(i%prim[j]==0)break;
            else mu[prim[j]*i]=-mu[i];
        }
    }
    sum[0]=0;
	for(int j=1;j<=cnt;j++)
        for(int i=1;i*prim[j]<=n;i++)
        pre[i*prim[j]]+=mu[i];   //预处理对于某一个数，它的质数因子的mu（T/i）和
    //为了得到区间i的mu和，通过求出前缀和来做
    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+(long long)pre[i];
}
int main()
{
    get_mu(10000000);
    int t,n,m;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d",&n,&m);
        int Min=min(n,m);
        long long ans=0;
        for(int l=1,r;l<=Min;l=r+1)
        {
            r=min(n/(n/l),m/(m/l)); //这里要注意
            ans+=1ll*(m/l)*(n/l)*(sum[r]-sum[l-1]); 
        }
        printf("%lld\n",ans);
    }
    return 0;
}

POI2007[ZAP-Queries](https://www.luogu.org/problemnew/show/P3455)

题意

求$\sum_{i=1}^{a}\sum_{j=1}^{b}[gcd(x,y)=d] $，多组，$1\le d\le a,b\le 50000$

思路

$$ ans=\sum_{i=1}^{a}\sum_{j=1}^{b}[gcd(x,y)=d] $$

我们假设：
$$
f(d)=\sum_{i=1}^{a}\sum_{j=1}^{b}[gcd(x,y)=d]
$$

$$
F(n)=\sum_{d|n}f(d)=\lfloor\frac{a}{n}\rfloor\lfloor\frac{b}{n}\rfloor
$$

$$
f(d)=\sum_{d|n}μ(\lfloor\frac{n}{d}\rfloor)F(n)
$$

发现原来f(d)=ans，那么有：
$$
ans=\sum_{d|n}μ(\lfloor\frac{n}{d}\rfloor)F(n)=\sum_{d|n}μ(\lfloor\frac{n}{d}\rfloor)\lfloor\frac{a}{n}\rfloor\lfloor\frac{b}{n}\rfloor
$$
枚举$\lfloor{\frac{n}{d}\rfloor}$=t：
$$
ans=\sum_{t=1}^{min(a,b)}μ(t)\lfloor\frac{a}{dt}\rfloor\lfloor\frac{b}{dt}\rfloor
$$
就能做了。

因为前面已经做了一道比这题难的，所以就不写O(n)的写法了，直接上代码，一遍过，爽歪歪。

看！代码

#include<bits/stdc++.h>
using namespace std;
const int N=50010;

bool vis[N];
int prim[N],pre[N],mu[N];
int cnt;
long long sum[N];

/** 对于μ函数的线筛*/
void get_mu(int n)
{
   // memset(vis,false,sizeof(vis));
    cnt=0;
	mu[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i]){mu[i]=-1;prim[++cnt]=i;}
        for(int j=1;j<=cnt&&prim[j]*i<=n;j++)
        {
            vis[i*prim[j]]=1;
            if(i%prim[j]==0)break;
            else mu[prim[j]*i]=-mu[i];
        }
    }
    sum[0]=0;
    for(int i=1;i<=n;i++)
        pre[i]+=mu[i]; 
    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+(long long)pre[i];
}
int main()
{
    get_mu(50000);
    int t,a,b,d;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d%d",&a,&b,&d);
        int Min=min(a,b);
        long long ans=0;
        for(int l=1,r;l<=Min;l=r+1)
        {
            r=min(a/(a/l),b/(b/l));
            ans+=1ll*(a/(l*d))*(b/(l*d))*(sum[r]-sum[l-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}

约数个数和

题意

设d(x)为x的约数个数，给定N、M，求$\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)$

多组，N,M。1<=N, M<=50000

思路

关于这个约数的个数d(x)，它其实有下面这样的函数式：
$$
d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]
$$
接下来套路的设：
$$
f(d)=\sum_{i=1}^N\sum_{j=1}^{M}[gcd(i,j)==d]
$$

$$
F(n)=\sum_{n|d}f(d)=\lfloor\frac{N}{n}\rfloor\lfloor\frac{M}{n}\rfloor –①
$$

由莫比乌斯反演可以得到：
$$
f(n)=\sum_{n|d}μ(\frac{d}{n})F(n)
$$
因为：
$$
ans=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]
$$
由于$∑_{d|n}μ(d)$的性质，即d=1的时候=1，其他都是0，可以将式子化为：

$$
ans=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j} \sum_{d|gcd(x,y)}μ(d)
$$
枚举d|gcd(x,y)的情况，d从1取到min(N,M)
$$
ans=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j} \sum_{d=1}^{min(N,M)}μ(d)*[d|gcd(x,y)]
$$
由于μ(d)与x，y，i，j无关，所以可以提到最前面
$$
ans=\sum_{d=1}^{min(N,M)}μ(d) \sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j} [d|gcd(x,y)]
$$
对于x|i，y|j的情况，可以由①得到$=\lfloor\frac{N}{x}\rfloor\lfloor\frac{M}{y}\rfloor$

由枚举i,j和它们的约数改变为枚举它们的约数再直接乘上这些约数的倍数的个数。因为每一个约数都会对它的倍数产生贡献。

$$
c
$$

为了将d|gcd(x,y)这个条件消去，我们更换枚举项x->dx，y->dy：
$$
ans=\sum_{d=1}^{min(N,M)}μ(d) \sum_{x=1}^{\frac{N}{d}}\sum_{y=1}^{\frac{M}{d}}\lfloor\frac{N}{dx}\rfloor\lfloor\frac{M}{dy}\rfloor
$$

$$
ans=\sum_{d=1}^{min(N,M)}μ(d) \sum_{x=1}^{\frac{N}{d}}\lfloor\frac{N}{dx}\rfloor\sum_{y=1}^{\frac{M}{d}}\lfloor\frac{M}{dy}\rfloor
$$

就可以做了。

#include <bits/stdc++.h>
using namespace std;
int mu[50010];
int prim[50010];
bool vis[50010];
long long sum[50010],chu[50010];
int cnt;
/** 对于μ函数的线筛 */
void get_mu(int n)
{
	mu[1]=1;
	for(int i=2;i<=n;i++)
	{
    	if(!vis[i]){prim[++cnt]=i;mu[i]=-1;}
   	 	for(int j=1;j<=cnt&&prim[j]*i<=n;j++)
        {
        	vis[prim[j]*i]=1;
        	if(i%prim[j]==0)break;
        	else mu[i*prim[j]]=-mu[i];
    	}
	}
	//这里是key
    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+mu[i]; //第一个sum
    for(int i=1;i<=n;i++)
    {

        long long ans=0;
        for(int l=1,r;l<=i;l=r+1)
        {
            r=i/(i/l); 
            ans+=1ll*(r-l+1)*(i/l);
        }
        chu[i]=ans;//第2、3个sum，枚举n（n/d）来算
    }

}
int main()
{
    get_mu(50000);
    int t,n,m;
    scanf("%d",&t);
    while(t--)
    {
        long long ans=0;
        scanf("%d%d",&n,&m);
        int Min=min(n,m);
        for(int l=1,r;l<=Min;l=r+1)
        {
            r=min(n/(n/l),m/(m/l)); //右边的数应当取min，这里只有一个n就没关系
            ans+=(sum[r]-sum[l-1])*(chu[n/l])*(chu[m/l]);
            //三个sum相乘
        }
        printf("%lld\n",ans);
    }
    return 0;
}

problem b

题意

对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。
输入格式：
第一行一个整数n，接下来n行每行五个整数，分别表示a、b、c、d、k
输出格式：
共n行，每行一个整数表示满足要求的数对(x,y)的个数
1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000

思路

这里用了容斥的思想：
$$
ans=\sum_{i=a}^{b}\sum_{j=c}^{d}[gcd(i,j)==k]
$$
实际上，（画个图就明白了）
$$
ans=sum{(1,b)(1,d)}-sum(1,b)(1,c-1)-sum(1,a-1)(1,d)+sum(1,a-1)(1,c-1)
$$
所以还是老套路，ans的每一部分可以分别套入公式算出来，所以只要求出一个公式：
$$
ans=\sum_{t=1}^{min(X,Y)}μ(t)\lfloor\frac{X}{kt}\rfloor\lfloor\frac{Y}{kt}\rfloor
$$
将x、y分别用b、d；b、c-1；a-1、d；a-1，c-1带入公式加加减减就好了。

p.s.现在的题目都T的很。。。无聊。。。我把计算的式子用ans来接收返回的函数然后输出ans,T掉，直接输出结果。。。。就快一点。。话是这样说。。。但是这也太。。。。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int mu[50010];
int prim[50010];
bool vis[50010];
long long sum[50010],chu[50010];
int cnt;
/** 对于μ函数的线筛 */
void get_mu(int n)
{
    mu[1]=1;
    for(int i=2;i<=n;i++)
    {
    	if(!vis[i]){prim[++cnt]=i;mu[i]=-1;}
   	 	for(int j=1;j<=cnt&&prim[j]*i<=n;j++)
        {
        	vis[prim[j]*i]=1;
        	if(i%prim[j]==0)break;
        	else mu[i*prim[j]]=-mu[i];
    	}
    }
    //这里是key
    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+mu[i];

}
ll cal(int n,int m,int k)
{
    ll ans=0;
    int Min=min(n,m);
    for(int l=1,r;l<=Min;l=r+1)
    {
        r=min(n/(n/l),m/(m/l)); //右边的数应当取min，这里只有一个n就没关系
        ans+=(sum[r]-sum[l-1])*(n/(l*k))*(m/(l*k));
            //这里根据情况改啊！r-（l-1）是因为这里的积性函数可以看作是μ(i)=1，n/l是后面的那个n/i
    }
   // printf("%lld\n",ans);
    return ans;
}
int main()
{
    get_mu(50000);
    int t,n,m;
    scanf("%d",&t);
    while(t--)
    {
        int a,b,c,d,k;
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        printf("%lld\n",cal(b,d,k)-cal(a-1,d,k)-cal(b,c-1,k)+cal(a-1,c-1,k));
    }
    return 0;
}

crash的数字表格

题意

求$Ans=\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) $

输入格式：
输入的第一行包含两个正整数，分别表示N和M。
输出格式：
输出一个正整数，表示表格中所有数的和mod20101009的值。

思路

由公式：$lcm(i,j)=\frac{ij}{gcd(i,j)}$可以得到，
$$
ans=\sum_{i=1}^{N}\sum_{j=1}^{M}\frac{ij}{gcd(i,j)}
$$
gcd在分母中不好，枚举gcd=d，将其放到前面来：
$$
ans=\sum_{d=1}^{min(N,M)}\sum_{i=1}^{N}\sum_{j=1}^{M}\frac{ij}{d}[gcd(i,j)=d]
$$
d在分母中又不爽，还有gcd，枚举gcd（让i=di，j=dj）：
$$
ans=\sum_{d=1}^{min(N,M)}d\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}ij[gcd(i,j)=1]
$$
由于性质$\sum_{x|n}μ(x) = (n=1)$
$$
ans=\sum_{d=1}^{min(N,M)}d\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}ij \sum_{x|gcd(i,j)}μ(x)
$$
要消去x|gcd(i,j)，因此枚举x：
$$
ans=\sum_{d=1}^{min(N,M)}d\sum_{x=1}^{min(\lfloor\frac{N}{d}\rfloor,\lfloor\frac{M}{d}\rfloor)}μ(x)\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}ij[x|gcd(i,j)]
$$
为了消去那个[]的条件，我们让i=xi，j=xj：
$$
ans=\sum_{d=1}^{min(N,M)}d\sum_{x=1}^{min(\lfloor\frac{N}{d}\rfloor,\lfloor\frac{M}{d}\rfloor)}x^2μ(x)\sum_{i=1}^{\lfloor\frac{N}{dx}\rfloor}i\sum_{j=1}^{\lfloor\frac{M}{dx}\rfloor}j
$$

展开全文 >>

Color it (简单扫描)

2018-08-09

Color it

题意

给你一个表格样，再给出m个圆，满足式子：$\sqrt{(i-x_c)^2+(j-y_c)^2}\le{r}$ 的，都要染成黑色，问你最后没有被染色的格子数。

思路就是简单模拟一下扫描线的操作，可以先对每一个圆，处理出圆内列对应的扫描线的长度，存在相应的列中。然后再遍历每一列，对于列中存在的扫描线，按照l小r大排序，然后如果这个区间的l要大于上个区间的r，ans就减去this.r-this.l+1，否则如果this.l<last.r,并且this.r>last.r，就this.r-last.r（不加1因为last.r已经算过），如果this.r<last.r，则不必再计算。

看！代码

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1e4+10;
struct node{
    int up,down;
    node(int x,int y){
        up=x;down=y;
    }
};
bool cmp(node a,node b)
{
    if(a.up==b.up) return a.down>b.down;
    return a.up<b.up;
}

int main()
{
    int t;
    int n,m,q;
    scanf("%d",&t);
    while(t--)
    {
        int x,y,r,ux,dx,ly,ry;
        scanf("%d%d%d",&n,&m,&q);
        vector <node> v[MAXN*2];
        for(int i=0;i<q;i++)
        {
            scanf("%d%d%d",&x,&y,&r);
            ly=max(0,y-r);
            ry=min(y+r,m-1);
            for(int j=ly;j<=ry;j++) //先处理出每列对应的扫描线，存下来
            {
                int derta=sqrt(r*r-(y-j)*(y-j));
                ux=max(0,x-derta);
                dx=min(n-1,x+derta);
                v[j].push_back(node(ux,dx));
            }
        }
        int len;
        int last,ans=n*m;
        for(int i=0;i<m;i++) //对应每一列存在的扫描线，先从l小，r大排序，然后判断计算
        {
            len=v[i].size();
            if(len)
            {
                last=-1;
                sort(v[i].begin(),v[i].end(),cmp);
                for(int j=0;j<len;j++)
                {
                    if(v[i][j].down<=last)
                        continue;
                    if(v[i][j].up>last)
                        ans-=v[i][j].down-v[i][j].up+1;
                    else if(v[i][j].down>last)
                        ans-=v[i][j].down-last;
                    last=v[i][j].down;
                }
            }
        }
        printf("%d\n",ans);
    }
    return 0;

}

展开全文 >>

中国剩余定理模板

2018-08-08

#include <iostream>
using namespace std;
 
int Extended_Euclid(int a,int b,int &x,int &y)    //扩展欧几里得算法
{
	int d;
	if(b==0)
	{
		x=1;y=0;
		return a;
	}
	d=Extended_Euclid(b,a%b,y,x);
	y-=a/b*x;
	return d;
}
 
int Chinese_Remainder(int a[],int w[],int len)    //中国剩余定理  a[]存放余数  w[]存放两两互质的数
{
	int i,d,x,y,m,n,ret;
	ret=0;
	n=1;
	for (i=0;i<len;i++)
		n*=w[i];
	for (i=0;i<len;i++)
	{
		m=n/w[i];
		d=Extended_Euclid(w[i],m,x,y);
		ret=(ret+y*m*a[i])%n; //注意这里如果会超出longlong的范围就要用快速乘
	}
	return (n+ret%n)%n;
}
 
 
int main()
{
	int n,i;
	int w[15],b[15];
	while (scanf("%d",&n),n)   
	{
		for (i=0;i<n;i++)
		{
			scanf("%d%d",&w[i],&b[i]);
		}
		printf("%d/n",Chinese_Remainder(b,w,n));
	}
	return 0;
}

展开全文 >>

Lucas 模板

2018-08-08

#include <bits/stdc++.h>
using namespace std;
 
typedef long long ll;
 
ll quick_mod(ll a,ll b,ll m)
{
    ll ans = 1;
    a %= m;
    while(b)
    {
        if(b&1)
            ans = ans * a % m;
        b >>= 1;
        a = a * a % m;
    }
    return ans;
}
 
ll getC(ll n, ll m,ll mod)
{
    if(m > n)
        return 0;
    if(m > n-m)
        m = n-m;
    ll a = 1,b = 1;
    while(m){
        a = (a*n)%mod;
        b = (b*m)%mod;
        m--;
        n--;
    }
    return a*quick_mod(b,mod-2,mod)%mod;
}
 
ll Lucas(ll n,ll k,ll mod)
{
    if(k == 0)
        return 1;
    return getC(n%mod,k%mod,mod)*Lucas(n/mod,k/mod,mod)%mod;
}
 
int main()
{
    int T;
    scanf("%d",&T);
    while(T--){
        ll n,m,mod;
        scanf("%lld%lld%lld",&n,&m,&mod);
        printf("%lld\n",Lucas(n+m,m,mod));
    }
    return 0;
}

展开全文 >>

Ponds （并查集||拓扑+深搜）

2018-08-07

Ponds

题意

给你n个点，m条无向边，每个点有vi的价值，你必须尽可能删去点，要求是删去的点的度数必须小于2。问你最后每个连通图中剩余奇数个点的总价值？

思路

1、一开始想是并查集吧，因为要计算每个环的点的数量，先把所有能删的点删去之后，再把剩下的边union起来，计算每个父亲所拥有的子的个数，奇数就加上去……

2、拓扑+深搜。应该还是把能删掉的点先删掉，然后深搜也是找环的点数。

看！代码

/**
	并查集做法
	*/
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const long long INF=0x3f3f3f3f;
const int N = 1e4+10;
struct node{
    vector<int>next;
    int deg;
    int id;
    node()
    {
        deg=0;
    }
};
int id[N];
int sz[N];

ll poor[N];
int vis[N];
int cnt[N];
int edge[N][2];
int Find(int x)  ///找跟根结点
{
    while(x!=id[x])
    {
        id[x]=id[id[x]];
        x=id[x];
    }
    return x;
}

int un(int p,int q)
{
    int pr=Find(p);
    int qr=Find(q);
    if(pr==qr) return -1;
    if(sz[pr]<=sz[qr])
    {
        sz[qr]+=sz[pr];id[pr]=qr;
    }
    else
    {
        sz[pr]+=sz[qr];id[qr]=pr;
    }
}

void clear()
{
    for(int i=1;i<=N;i++)
    {
        id[i]=i;sz[i]=1;
    }
}

int main()
{
    int t;
    int p,m,u,v;
    scanf("%d",&t);
    while(t--)
    {
        node x[N];
        memset(vis,0,sizeof(vis));
        memset(cnt,0,sizeof(cnt));
        clear();

        scanf("%d%d",&p,&m);
        for(int i=1;i<=p;i++)
            scanf("%lld",&poor[i]);

        for(int i=1;i<=p;i++)
            x[i].id=i;

        for(int i=0;i<m;i++)
        {
            scanf("%d%d",&u,&v);
            x[u].next.push_back(v);
            x[u].deg++;
            x[v].next.push_back(u);
            x[v].deg++;
            edge[i][0]=u;
            edge[i][1]=v;
        }
        queue<node>q;
        for(int i=1;i<=p;i++)
            if(x[i].deg==1)
                q.push(x[i]);
            else if(x[i].deg==0)
                vis[i]=1;
        node temp;
        while(!q.empty())
        {
            temp=q.front();
            q.pop();
            vis[temp.id]=1;
            for(int i=0;i<temp.next.size();i++)
                if(!vis[temp.next[i]])
            {
                x[temp.next[i]].deg--;
                if(x[temp.next[i]].deg==1)
                    q.push(x[temp.next[i]]);
                else if(x[temp.next[i]].deg==0)
                    vis[x[temp.next[i]].id]=1;
            }
        }

        for(int i=0;i<m;i++)
        {
            int x=edge[i][0],y=edge[i][1];
            if(vis[x]==1||vis[y]==1)continue;
            un(x,y);
        }
        for(int i=1;i<=p;i++)
        {
            if(vis[i]==0)
                cnt[Find(i)]++;
        }

        long long ans=0;
        for(int i=1;i<=p;i++)
        {
            if(vis[i]==0&&cnt[Find(i)]>=3&&cnt[Find(i)]%2==1)
                ans+=poor[i];
        }
        printf("%lld\n",ans);
    }
    return 0;
}

展开全文 >>

Fruit Ninja(头一次接触的随机化算法)

2018-08-05

Fruit Ninja

题意

屏幕上有n个水果，你可以画一条直线切掉其中一些，给你N和x，M/N>=x，M为切掉的水果数，x为0~1的纯小数，指的是切掉的水果比吧。问你是否有满足这个式子的这样一条直线，有输出yes，没有输出no。

N<=1e4 , x小数点后一位（0,1）。

思路

随机化，巧妙巧妙，厉害厉害，佩服佩服，不过要分析一波随机的数量。由于m是切掉的水果数，n是总水果数，x是两者之比，也就是说在这条直线上的水果数是m个，意味着随机选一个水果，刚好在直线上的概率是x，因此随机选两个刚好组成这条直线的概率是x^2，而x∈[0.1,0.9]，因此概率最小为0.01，也就是说取100次就能取到那条我们需要的直线，反之如果没有满足的，就说明“NO”。为了保险起见，我们将随机范围设置为200。

看！代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
struct node{
    ll x,y;
};
node point[10010];
int main()
{
    int t;
    srand((unsigned)time(NULL));
    scanf("%d",&t);
    while(t--)
    {
        int n;
        double x;
        scanf("%d%lf",&n,&x);
        for(int i=0;i<n;i++)
        {
            scanf("%lld%lld",&point[i].x,&point[i].y);
        }
        int flag=0;
        for(int i=0;i<200;i++)
        {
            int a=rand()%n;
            int b=rand()%n;
            while(a==b)
                a=rand()%n;

            int cnt=0;
            for(int j=0;j<n;j++)
                if((point[a].x-point[b].x)*(point[a].y-point[j].y)==(point[a].y-point[b].y)*(point[a].x-point[j].x))
                    cnt++;

            if(cnt>=n*x)
            {
                flag=1;break;
            }
        }
        if(flag)    puts("Yes");
        else puts("No");
    }
    return 0;
}

展开全文 >>