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

Description
给一个非负序列$x_{1,2,...,n}$，值小于32677，令$y_{i,j}=x_{i}*x_{j} \ mod\ 32677$，问有多少六元组$(a,b,c,d,e,f)$满足$gcd(y_{a,b},y_{c,d})=gcd(y_{c,d},y_{e,f})=gcd(y_{e,f},y_{a,b})=1$，结果模$2^{30}$
Input
第一行一整数T表示用例组数，每组用例首先输入一整数n表示序列长度，之后输入n个非负整数$x_{1,2,...,n}(1\leq T \leq5，1\leq n \leq2\times10^{5}，0\leq x_{1,2,...,n} <32677$)
Output
对于每组用例，输出满足条件的六元组数，结果模$2^{30}$
Sample Input
2
1
1
5
1 2 3 4 5
Sample Output
1
1087
Solution

令$nx[k]​$为$x_{i}=k​$的个数，$ny[k]​$为$y_{i,j}=k​$的个数，则有$ny[k]=\sum\limits_{i\times j \equiv k\ mod\ 32677 }nx[i]\times nx[j]​$
考虑到$p=32677=41\times 797$，$p1=41$和$p2=797$都是素数，用中国剩余定理拆开

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

$b[s][t]=\sum [i\times j\equiv s\ mod\ p1,k\times l\equiv t\ mod\ p2]a [i][k]\times a[j][l]$

$a[i][j]=|\{x_{k}|x_{k}\equiv i\ mod\ p1,x_{k}\equiv j\ mod\ p2\}|$

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

一遍二维FFT即可得到$ny$数组，令$sum(x)=\sum\limits_{i|j,1\leq j<p}ny[j]$，对原式反演并化简得

$$\begin{eqnarray} ans &= &\sum\limits_{1\leq a,b,c <p} [(a,b)=(a,c)=(b,c)=1]ny[a]ny[b]nz[c]\\ &=&\sum\limits_{1\leq a,b,c<p}ny[a]ny[b]nz[c]\sum\limits_{x|(a,b)}\mu(x)\sum\limits_{y|(a,c)}\mu(y)\sum\limits_{z|(b,c)}\mu(z)\\ &=&\sum\limits_{1\leq x,y,z<p}\mu(x)\mu(y)\mu(z)\sum\limits_{[x,y]|a}ny[a]\sum\limits_{[x,z]|b}ny[b]\sum\limits_{[y,z]|c}ny[c]\\ &=&\sum\limits_{1\leq x,y,z<p}\mu(x)\mu(y)\mu(z)sum([x,y])sum([x,z])sum([y,z])\\ \end{eqnarray}$$
注意到 $\mu(x)\neq0,\mu(y)\neq0,[x,y]<p$的 $(x,y)$对不是很多，而且三元组 $(x,y,z)$需两两满足该条件，

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

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

2.$(a,b,c)$全部相等的情况，对答案贡献是$\sum\limits_{1\leq i<p}\mu^{3}(i)sum^{3}(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+5
const 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;
}