首先推一下式子,变成求
1nm∑(ak1ak2..akm mod P)=1nm∑kk∑[ak1ak2..akm mod P==k]
注意模 P 是在
容易想到 DP , fij %p←f′i×cntj 。做 m 次。其中
如果借助原根,转化成 g 的次幂,就变成
其实不用原根转化,也可以直接类似快速幂的搞,叫倍增也行吧…
如果考虑用矩阵乘法优化 DP ,直接做会T。但发现转移矩阵是个循环矩阵。就可以只维护一行,每次乘 O(P2) 。本质上和上面的做法是一样的(写出来也一样),只是如果直接矩阵乘就做了很多无用功。
#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;
}