随(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 1005using 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;}