• 2018秦皇岛ccpc-camp Steins;Gate (原根+FFT)


    因为给定的模数P保证是素数,所以P一定有原根.
    根据原根的性质,若(g)(P)的原根,则(g^k)能够生成([1,P-1])中所有的数,这样的k一共有P-2个.
    (a_i*a_j(mod P)=a_k) 就可以转化为(g^i*g^j(mod P) = g^{i+j}(mod P)=g^k).
    问题转化为了求有多少对有序的<i,j>满足 ((i+j)(mod (P-1)) = k).
    求出原根后,对([1,P-1])中的每个数编号, 统计每个编号出现的次数,然后FFT求卷积
    要特判0,因为原根不会生成0.所以用总的有序对数-其他不含0的有序对数得到含0的有序对,这是0的答案;超过P-1的数肯定没有符合的有序对.

    #include <bits/stdc++.h>
    using namespace std;
    typedef long long LL;
    const int maxn = 2e5+5;
    bool notpri[maxn];
    int pri[maxn],zyz[maxn];
    typedef long long LL;
    
    void pre(int N){
        notpri[1]=1;
        for(int i=2;i<N;++i){
            if(!notpri[i]){
                pri[++pri[0]]=i;
            }
            for(int j=1;j<=pri[0] && (LL)i*(LL)pri[j]<N;++j){
                notpri[i*pri[j]]=1;
                if(i%pri[j]==0){
                    break;
                }
            }
        }
    }
    LL Quick_Pow(LL x,LL p,LL mod){
        if(!p){
            return 1ll;
        }
        LL res=Quick_Pow(x,p>>1,mod);
        res=res*res%mod;
        if((p&1ll)==1ll){
            res=(x%mod*res)%mod;
        }
        return res;
    }
    int FindRoot(int x){/*求素奇数的最小原根,倘若x不是奇数,但是也有原根的话,将质
    因子分解改成对phi(x)即可。倘若要求多个原根,直接接着暴力验证即可*/
        int tmp=x-1;
        for(int i=1;tmp && i<=pri[0];++i){
            if(tmp%pri[i]==0){
                zyz[++zyz[0]]=pri[i];
                while(tmp%pri[i]==0){
                    tmp/=pri[i];
                }
            }
        }
        for(int g=2;g<=x-1;++g){
            bool flag=1;
            for(int i=1;i<=zyz[0];++i){
                if(Quick_Pow((LL)g,(LL)((x-1)/zyz[i]),(LL)x)==1){
                    flag=0;
                    break;
                }
            }
            if(flag){
                return g;
            }
        }
        return 0;
    }
    
    const int MAXN = 4e5 + 10;
    const double PI = acos(-1.0);
    struct Complex{
        double x, y;
        inline Complex operator+(const Complex b) const {
            return (Complex){x +b.x,y + b.y};
        }
        inline Complex operator-(const Complex b) const {
            return (Complex){x -b.x,y - b.y};
        }
        inline Complex operator*(const Complex b) const {
            return (Complex){x *b.x -y * b.y,x * b.y + y * b.x};
        }
    } va[MAXN * 2 + MAXN / 2], vb[MAXN * 2 + MAXN / 2];
    int lenth = 1, rev[MAXN * 2 + MAXN / 2];
    int N, M;   // f 和 g 的数量
        //f g和 的系数
        // 卷积结果
        // 大数乘积
    int f[MAXN],g[MAXN];
    vector<LL> conv;
    vector<LL> multi;
    //f g
    void init()
    {
        int tim = 0;
        lenth = 1;
        conv.clear(), multi.clear();
        memset(va, 0, sizeof va);
        memset(vb, 0, sizeof vb);
        while (lenth <= N + M - 2)
            lenth <<= 1, tim++;
        for (int i = 0; i < lenth; i++)
            rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (tim - 1));
    }
    void FFT(Complex *A, const int fla)
    {
        for (int i = 0; i < lenth; i++){
            if (i < rev[i]){
                swap(A[i], A[rev[i]]);
            }
        }
        for (int i = 1; i < lenth; i <<= 1){
            const Complex w = (Complex){cos(PI / i), fla * sin(PI / i)};
            for (int j = 0; j < lenth; j += (i << 1)){
                Complex K = (Complex){1, 0};
                for (int k = 0; k < i; k++, K = K * w){
                    const Complex x = A[j + k], y = K * A[j + k + i];
                    A[j + k] = x + y;
                    A[j + k + i] = x - y;
                }
            }
        }
    }
    void getConv(){             //求多项式
        init();
        for (int i = 0; i < N; i++)
            va[i].x = f[i];
        for (int i = 0; i < M; i++)
            vb[i].x = g[i];
        FFT(va, 1), FFT(vb, 1);
        for (int i = 0; i < lenth; i++)
            va[i] = va[i] * vb[i];
        FFT(va, -1);
        for (int i = 0; i <= N + M - 2; i++)
            conv.push_back((LL)(va[i].x / lenth + 0.5));
    }
    
    LL vz[MAXN];
    int id[MAXN];
    int cnt[MAXN];
    int rnk[MAXN];
    LL ans[MAXN];
    
    int main()
    {
        #ifndef ONLINE_JUDGE
            freopen("in.txt","r",stdin);
            freopen("out.txt","w",stdout);
        #endif
        pre(maxn-1);
        int n,P ;
        scanf("%d %d",&n, &P);
        int rt = FindRoot(P);
        memset(id,0,sizeof(id));
    
    
        LL idx = 1;
        for(int i=0;i<P-1;++i){         //一共只有[0,P-2],P-1个id
            id[idx] = i;
            rnk[i]= idx;
            idx = idx*rt%P;
        }
    
        for(int i=1;i<=n;++i){
            LL tmp;
            scanf("%lld",&tmp);
            vz[i] = tmp;
            tmp%=P;
            if(tmp==0) continue;        //卷积中不考虑0的贡献
            cnt[id[tmp]]++;
        }
    
        N = M = P;
        for(int i=0;i<P-1;++i){
            f[i] = g[i] = cnt[i];
        }
    
        getConv();
        int sz = conv.size();
    
        for(int i=0;i<sz;++i){
            LL tmp = conv[i];
            ans[rnk[i%(P-1)]] += tmp;
        }
    
        LL tot = (LL)n*n;           //全部的枚举可能-不选0之外的组合 = 包含0的组合
        for(int i=1;i<P;++i){
            tot -= ans[i];
        }
        ans[0] = tot;
        for(int i=1;i<=n;++i){
            if(vz[i]>=P) printf("0
    ");
            else{
                printf("%lld
    ",ans[vz[i]]);
            }
        }
        return 0;
    }
    
    
  • 相关阅读:
    Git配置
    第一次作业
    第二次作业
    python目前最好用的IDE——pycharm
    九九乘法表
    python语言的优点和缺点
    Python高效编程的19个技巧
    Python中 filter | map | reduce | lambda的用法
    Python 代码优化常见技巧
    求逆序对
  • 原文地址:https://www.cnblogs.com/xiuwenli/p/9735804.html
Copyright © 2020-2023  润新知