HDU 6120 All Kill（数论+FFT+莫比乌斯反演）

Description
给一个非负序列x1,2,...,n，值小于32677，令yi,j=xi∗xj mod 32677，问有多少六元组(a,b,c,d,e,f)满足gcd(ya,b,yc,d)=gcd(yc,d,ye,f)=gcd(ye,f,ya,b)=1，结果模230
Input
第一行一整数T表示用例组数，每组用例首先输入一整数n表示序列长度，之后输入n个非负整数x1,2,...,n(1≤T≤5，1≤n≤2×105，0≤x1,2,...,n<32677)
Output
对于每组用例，输出满足条件的六元组数，结果模230
Sample Input
2
1
1
5
1 2 3 4 5
Sample Output
1
1087
Solution

令nx[k]​为xi=k​的个数，ny[k]​为yi,j=k​的个数，则有ny[k]=∑i×j≡k mod 32677nx[i]×nx[j]​
考虑到p=32677=41×797，p1=41和p2=797都是素数，用中国剩余定理拆开

ny[k]=b[k%p1][k%p2]

b[s][t]=∑[i×j≡s mod p1,k×l≡t mod p2]a[i][k]×a[j][l]

a[i][j]=|{xk|xk≡i mod p1,xk≡j mod p2}|

单独处理下标为0的a[i][j]，对于非零下标用原根的幂次表示即可转化为普通的二维循环卷积，做

一遍二维FFT即可得到ny数组，令sum(x)=∑i|j,1≤j<pny[j]，对原式反演并化简得

ans====∑1≤a,b,c<p[(a,b)=(a,c)=(b,c)=1]ny[a]ny[b]nz[c]∑1≤a,b,c<pny[a]ny[b]nz[c]∑x|(a,b)μ(x)∑y|(a,c)μ(y)∑z|(b,c)μ(z)∑1≤x,y,z<pμ(x)μ(y)μ(z)∑[x,y]|any[a]∑[x,z]|bny[b]∑[y,z]|cny[c]∑1≤x,y,z<pμ(x)μ(y)μ(z)sum([x,y])sum([x,z])sum([y,z])

注意到μ(x)≠0,μ(y)≠0,[x,y]<p的(x,y)对不是很多，而且三元组(x,y,z)需两两满足该条件，

故对满足该条件的(x,y)建边去找三元环即可，注意几种特殊情况

1.(a,b,c)中有一个为0的情况，对答案贡献是3×ny[0]×ny[1]×ny[1]

2.(a,b,c)全部相等的情况，对答案贡献是∑1≤i<pμ3(i)sum3(i)

3.(a,b,c)有两个相等的情况，在找满足条件的二元组时可以处理

Code

#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#include<vector>using namespace std;typedef long long ll;#define maxfft 65536+5const double pi=acos(-1.0);struct cp{    double a,b;    cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}    cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}    cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}    cp operator *(const double &o)const {return (cp){a*o,b*o};}    cp operator !() const{return (cp){a,-b};}}w[maxfft];int pos[maxfft];void fft_init(int len){    int j=0;    while((1<<j)<len)j++;    j--;    for(int i=0;i<len;i++)        pos[i]=pos[i>>1]>>1|((i&1)<<j);}void fft(cp *x,int len,int sta){    for(int i=0;i<len;i++)        if(i<pos[i])swap(x[i],x[pos[i]]);    w[0]=(cp){1,0};    for(unsigned i=2;i<=len;i<<=1)    {        cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};        for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];        for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;        for(int j=0;j<len;j+=i)        {            cp *a=x+j,*b=a+(i>>1);            for(int l=0;l<i>>1;l++)            {                cp o=b[l]*w[l];                b[l]=a[l]-o;                a[l]=a[l]+o;            }        }    }    if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;}cp x[maxfft],y[maxfft];void FFT(int *a,int n,ll *c)//自卷积 {    int len=1;    while(len<n)len<<=1;    fft_init(len);    for(int i=n/2;i<len;i++)x[i].a=x[i].b=0;    for(int i=0;i<n;i++)(i&1?x[i>>1].b:x[i>>1].a)=a[i];    fft(x,len,1);    for(int i=0;i<len/2;i++)    {        int j=len-1&len-i;        y[i]=x[i]*x[i]-(x[i]-!x[j])*(x[i]-!x[j])*(w[i]+(cp){1,0})*0.25;    }    for(int i=len/2;i<len;i++)    {        int j=len-1&len-i;        y[i]=x[i]*x[i]-(x[i]-!x[j])*(x[i]-!x[j])*((cp){1,0}-w[i^len>>1])*0.25;    }    fft(y,len,-1);    for(int i=0;i<2*n;i++)        if(i&1)c[i]=(ll)(y[i>>1].b+0.5);        else c[i]=(ll)(y[i>>1].a+0.5);}#define P 32677#define P1 41#define P2 797#define mod (1<<30)int r1,g1[P1],e1[P1];int r2,g2[P2],e2[P2];//r为模p原根,e[i]=r^i%p,i=r^g[i] void Deal(int p,int &r,int *e,int *g) {    r=2;    while(1)    {        for(int i=1;i<p;i++)g[i]=0;        e[0]=1;        for(int i=1;i<p;i++)        {            e[i]=e[i-1]*r%p;            if(g[e[i]])            {                e[0]=-1;                break;            }            g[e[i]]=i;        }        if(e[0]==1)break;        r++;    }    g[1]=0;}bool check[P];int prime[P];unsigned mu[P];void Moblus(int n){    memset(check,0,sizeof(check));    mu[1]=1;    int tot=0;    for(int i=2;i<=n;i++)    {        if(!check[i])prime[tot++]=i,mu[i]=-1;        for(int j=0;j<tot;j++)        {            if(i*prime[j]>n)break;            check[i*prime[j]]=1;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            mu[i*prime[j]]=-mu[i];        }    }}int T,n,a[P1][P2],A[maxfft],c[P];ll b[P1][P2],B[maxfft<<1];int du[P];unsigned sum[P],ans;typedef pair<int,unsigned>PP;vector<PP>v[P];int main(){    Moblus(P-1);    Deal(P1,r1,e1,g1),Deal(P2,r2,e2,g2);    scanf("%d",&T);    while(T--)    {        scanf("%d",&n);        memset(a,0,sizeof(a));        memset(b,0,sizeof(b));        while(n--)        {            int temp;            scanf("%d",&temp);            a[temp%P1][temp%P2]++;        }        memset(A,0,sizeof(A));        for(int i=0;i<P1;i++)            for(int j=0;j<P2;j++)                A[2*P2*g1[i]+g2[j]]=a[i][j];        FFT(A,2*P,B);        memset(b,0,sizeof(b));        for(int i=0;i<4*P;i++)            b[e1[i/(2*P2)%(P1-1)]][e2[i%(2*P2)%(P2-1)]]+=B[i];        for(int i=0;i<P1;i++)        {            ll t=a[i][0];            for(int j=1;j<P2;j++)t+=2ll*a[i][j];            for(int j=0;j<P1;j++)b[i*j%P1][0]+=t*a[j][0];        }        for(int i=1;i<P2;i++)        {            ll t=a[0][i];            for(int j=1;j<P1;j++)t+=2ll*a[j][i];            for(int j=1;j<P2;j++)b[0][i*j%P2]+=t*a[0][j];        }        for(int i=0;i<P;i++)c[i]=b[i%P1][i%P2]%mod;        ans=3*c[0]*c[1]*c[1];        memset(sum,0,sizeof(sum));        memset(du,0,sizeof(du));        for(int i=1;i<P;i++)            for(int j=i;j<P;j+=i)sum[i]+=c[j];        for(int i=1;i<P;i++)ans+=mu[i]*mu[i]*mu[i]*sum[i]*sum[i]*sum[i];        for(int gcd=1;gcd<P;gcd++)            for(int lcm=gcd;lcm<P;lcm+=gcd)                if(mu[lcm])                {                    int m=lcm/gcd;                    for(int j=1;j*j<m;j++)                        if(m%j==0)                        {                            int x=j*gcd,y=m/j*gcd;                            du[x]++,du[y]++;                            ans+=3*mu[x]*mu[x]*mu[y]*sum[x]*sum[lcm]*sum[lcm];                            ans+=3*mu[x]*mu[y]*mu[y]*sum[y]*sum[lcm]*sum[lcm];                        }                }        for(int i=1;i<P;i++)v[i].clear();        for(int gcd=1;gcd<P;gcd++)            for(int lcm=gcd;lcm<P;lcm+=gcd)                if(mu[lcm])                {                    int m=lcm/gcd;                    for(int j=1;j*j<m;j++)                        if(m%j==0)                        {                            int x=j*gcd,y=m/j*gcd;                            if(du[x]>du[y])swap(x,y);                            v[x].push_back(PP(y,sum[lcm]));                        }                }        for(int i=1;i<P;i++)sort(v[i].begin(),v[i].end());        for(int i=1;i<P;i++)            for(int l=0;l<v[i].size();l++)            {                int j=v[i][l].first,p1=0,p2=0;                unsigned sum1=v[i][l].second,sum2,sum3;                while(p1<v[i].size()&&p2<v[j].size())                {                    if(v[i][p1].first==v[j][p2].first)                    {                        int k=v[i][p1].first;                        sum2=v[i][p1].second,sum3=v[j][p2].second;                         ans+=6*mu[i]*mu[j]*mu[k]*sum1*sum2*sum3;                        p1++,p2++;                    }                    else if(v[i][p1].first<v[j][p2].first)p1++;                    else p2++;                }            }        printf("%d\n",ans%mod);    }    return 0;}