随 (rand) by liu_runda

擦汗。。。最后是颓代码过的

尽管之前已经讲过什么矩阵优化，除了优化费波数（早忘得一干二净），就没有码过了（线性代数就做了一道题）。

考虑30%（数据3，4，5）的dp[i][j]表示i次操作后，x为j的方案数

for(int i=1;i<mod;++i)dp[1][i]=cnt[i]%MOD;
//1预处理
for(int i=2;i<=m;++i)
	for(int s=1;s<mod;++s)
		for(int k=1;k<mod;++k)
			dp[i][k*s%mod]=(dp[i][k*s%mod]+dp[i-1][k]*cnt[s]%MOD)%MOD;
for(int i=1;i<mod;++i)
	ans=(ans+dp[m][i])%MOD;
ans=ans*qpow(n,m*(MOD-2)%(MOD-1))%MOD;

为O(mod^2*m)，mod<=1000还是很友善的，然鹅m暴到了10^9，优化这个m没什么好方法，往二进制上靠就完事了。快速幂可以有一些启示，任意一个数a都可以拆成2^k1+2^k2+...+2^kn。考虑将m也拆成这种形式，相当于式子中预处理cnt[k] s 次选了后的方案数。此时cnt[]->cnt[][]。

对于矩阵乘，可以省去一位，因为这一维一直是1

倍增也能过

从lockey那里颓来的代码（捂脸）

#include<cstdio>
#include<iostream>
#include<cstring>
#define ll long long
#define MOD 1000000007
using namespace std;
ll n,m,mod;
ll cnt[1010];
ll dp[1010];
ll qpow(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=ans*a%MOD;
		b>>=1;
		a=a*a%MOD;
	}
	return ans%MOD;
}
ll tmp[1010];
void x_dp(){
	for(ll i=1;i<mod;++i)
		for(ll j=1;j<mod;++j)
			tmp[i*j%mod]=(tmp[i*j%mod]+dp[i]*cnt[j]%MOD)%MOD;
	memcpy(dp,tmp,sizeof(tmp));
	memset(tmp,0,sizeof(tmp));
}
void x_self(){
	for(ll i=1;i<mod;++i)
		for(ll j=1;j<mod;++j)
			tmp[i*j%mod]=(tmp[i*j%mod]+cnt[i]*cnt[j]%MOD)%MOD;
	memcpy(cnt,tmp,sizeof(tmp));
	memset(tmp,0,sizeof(tmp));
}
int main(){
//	freopen("da.in","r",stdin);
	scanf("%lld%lld%lld",&n,&m,&mod);
	ll rn;
	for(int i=1;i<=n;++i){
		scanf("%lld",&rn);
		++cnt[rn];
	}
	memcpy(dp,cnt,sizeof(cnt));
	ll s=m;
	--m;
	while(m){
		if(m&1)x_dp();
		x_self();
		m>>=1;
	}
	ll res=0;
	for(ll i=1;i<mod;++i)
		res=(res+i*dp[i]%MOD)%MOD;
	printf("%lld\n",res*qpow(n,s*(MOD-2)%(MOD-1))%MOD);
}