[循环卷积 || 循环矩阵] HHHOJ #120. 随

首先推一下式子，变成求

$$\frac{1}{n^m}\sum(a_{k_1}a_{k_2}..a_{k_m} \text{ mod } P)=\frac{1}{n^m}\sum_{k} k\sum[a_{k_1}a_{k_2}..a_{k_m} \text{ mod } P==k]$$
注意模 $P$ 是在 $\sum$ 里面的。现在只需求 $F_k=\sum[a_{k_1}a_{k_2}..a_{k_m} \text{ mod } P==k]$。
容易想到 $DP$, $f_{ij\text{ %}p} \leftarrow f'_{i} \times cnt_j$ 。做 $m$ 次。其中 $cnt_i$ 表示 $a$ 中 $i$ 的个数。
如果借助原根，转化成 $g$ 的次幂，就变成 $f_{(i+j)\text{%}(P-1)} \leftarrow f'_{i} \times cnt_j$。本质上就是做 $m$ 次循环卷积。快速幂+暴力卷就可以 $O(P^2 \log m)$。如果用 $FFT$ 求卷积就可以 $O(P \log P\log m)$。
其实不用原根转化，也可以直接类似快速幂的搞，叫倍增也行吧…
如果考虑用矩阵乘法优化 $DP$，直接做会T。但发现转移矩阵是个循环矩阵。就可以只维护一行，每次乘 $O(P^2)$。本质上和上面的做法是一样的(写出来也一样)，只是如果直接矩阵乘就做了很多无用功。

#include<cstdio>
#include<algorithm>
using namespace std;
const int MOD=1e9+7,maxp=1005;
typedef long long LL;
LL n,m,P,g,ans,gp[maxp],n_gp[maxp];
LL Pow(LL a,LL b){
    LL res=1;
    for(;b;b>>=1,a=a*a%MOD) if(b&1) res=(res*a)%MOD;
    return res;
}
int vis[maxp],clk;
bool check(int g){
    clk++;
    for(LL i=1,now=g;i<=P-1;i++,now=(now*g)%P){
        if(vis[now]==clk) return false;
        gp[i]=now; n_gp[now]=i; vis[now]=clk;
    } 
    return true;
}
const int maxn=305;
int c[maxn];
inline void Mul(int a[],int b[],int n){
    for(int i=0;i<n;i++){
        c[i]=0;
        for(int j=0;j<n;j++) (c[i]+=(LL)a[j]*b[(i-j+n)%n]%MOD)%=MOD;
    }
    for(int i=0;i<n;i++) a[i]=c[i];
}
int cnt[maxp],T[maxn],res[maxn],T_n,T_m,res_n,res_m;
int main(){
    freopen("A.in","r",stdin);
    freopen("A.out","w",stdout); 
    scanf("%lld%lld%lld",&n,&m,&P);
    for(g=2;g<=P-1;g++) if(check(g)) break;
    gp[0]=1; n_gp[1]=0;
    for(int i=1;i<=n;i++){
        int x; scanf("%d",&x);
        cnt[n_gp[x]]++;
    }
    for(int i=0;i<=P-2;i++) T[i]=cnt[i];
    res[0]=1; 
    for(int b=m;b;b>>=1,Mul(T,T,P-1)) if(b&1) Mul(res,T,P-1);
    for(int i=0;i<=P-2;i++) (ans+=(gp[i]*res[i])%MOD)%=MOD;
    printf("%lld\n",((ans*Pow(Pow(n,m),MOD-2)%MOD)+MOD)%MOD);
    return 0;
}