随(rand) （概率dp）

9.23

思路：
考察概率和期望的求解,矩阵乘法,原根的性质,循环矩阵的性质,倍增优化DP.
首先需要注意到虽然n可以达到10^5,但相同数字可以合并考虑,只需要考虑mod个不同的数字及选择它们的概率.
第1个测试点:mod=2,则n个数字都是1,直接输出1即可.
第2个测试点:每次乘上去的数字只有一种选择,快速幂即可.
第3,4,5个测试点:定义f[i][j]表示i次操作后x的数值为j的概率.直接转移,复杂度O(m*mod^2)
第6,7,8个测试点:第3,4,5个测试点中的DP转移可以转化为矩阵乘法形式,利用矩阵快速幂进行优化,复杂度O(mod^3*logm)
第9,10个测试点(标算):利用原根进行转化,则乘法转化为加法,f[i][j]表示i次操作后x取模后等于原根的j次方的概率.指数需要对(mod-1)取模.这样转化一下我们发现转移还是矩阵的形式,而且是循环矩阵的形式.循环矩阵快速幂,复杂度O(mod^2*logm)
另一种标算:定义f[i][j]表示i次操作后变成原根的j次方的概率.求出g[i][j]表示2^i次操作后变成原根的j次方的概率.倍增的思想求出f[m][]这个数组.也是O(mod^2*logm)
更加优越的算法:本质上我们要做的是循环卷积,可以使用fft.但本题的模数使得fft较为不方便..

#include <cstdio>
#include <iostream>
#define MOD 1000000007
#define N 1005
using namespace std;

int n, m, mod, mod2, rt;

int pw[N], glog[N];
int cnt[N];
int P[33][N], f[2][N];

int qpow(int a, int x, int mod){
    int ans = 1;
    for( ; x; x>>=1, a=a*1ll*a%mod)
        if(x & 1) ans = ans * 1ll * a % mod;
    return ans;
}

int getroot(int x){
    for(int i=1; i<=x; ++i){
        int tmp = 1;
        bool flag = true;
        for(int j=1; j<x-1; ++j){
            tmp = tmp * 1ll * i % mod;
            if(tmp == 1){
                flag = false; break;
            }
        }
        if( flag ) return i;
    }
}

int main(){
    freopen ("rand.in", "r", stdin);
    freopen ("rand.out", "w", stdout);
    scanf("%d%d%d", &n, &m, &mod);
    rt = getroot(mod);
    pw[0] = 1;
    for(int i=1; i<=mod; ++i){
        pw[i] = pw[i-1] * 1ll * rt % mod;
    }//pw[i]记录 rt^i mod 的 ans 
    for(int i=0; i<mod-1; ++i) glog[pw[i]] = i;
    mod2 = mod - 1;
    for(int i=1; i<=n; ++i){
        int x; scanf("%d", &x);
        cnt[glog[x]]++;
    }
    int invn = qpow(n, MOD-2, MOD);//a*(b^(10^9+5))模10^9+7  10^9+7是P  10^9+5是P-2  所以也就是求a * invb 
    for(int i=0; i<mod2; ++i){
        P[0][i] = invn * 1ll * cnt[i] % MOD;//2^0次操作之后 同余于rt^i的概率 
    }
    for(int i=1; i<=32; ++i){
        for(int j=0; j<mod2; ++j){
            for(int k=0; k<mod2; ++k){
                P[i][(j+k)%mod2] = ( P[i][(j+k)%mod2] + P[i-1][j] * 1ll * P[i-1][k] % MOD ) % MOD;
            }
        }
    }
    int p = 0; f[0][0] = 1;//初始值为1 
    for(int i=0; m; m>>=1, ++i){
        if(m & 1){
            p ^= 1;
            for(int j=0; j<mod2; ++j) f[p][j] = 0;//滚动数组清零 
            for(int j=0; j<mod2; ++j){
                for(int k=0; k<mod2; ++k){
                    f[p][(j+k)%mod2] = ( f[p][(j+k)%mod2] + f[p^1][j] * 1ll * P[i][k] % MOD ) % MOD;
                }//m次操作之后 同余于rt^i的概率 
            }
        }
    }
    for(int j=0; j<mod2; ++j) f[0][j] = f[p][j];
    int ans = 0;
    for(int i=0; i<mod2; ++i){
        ans = (ans + f[0][i] * 1ll * pw[i] % MOD) % MOD;//概率乘权值 
    }
    printf("%d\n", ans);
    return 0;
}