随 (rand) (dp，原根，矩阵/倍增优化)

随 (rand):

　　题干：

　　题解：

　　　　对于$10\%$的数据mod=2,a[i]必定为1，直接输出1.

　　　　对于另$10\%$的数据n=1,输出$ a[1]^m\%mod $ 注意模mod而非$1e9+7$。

　　　　对于$50\%$的数据，我们考虑dp，

　　　　　　设$ f[i][j] $为i次操作后结果为j的方案数，

　　　　　　易得转移方程为：$ f[i][j*a[k]\%mod]=\sum\limits f[i-1][j] $(考场上想到的sb DP)

　　　　　　但$ O(n*m*mod) $的复杂度令人望而却步，并不能多得分。

　　　　　　于是考虑优化，观察得$mod<=300$而$n<=100000$,想到类似离散化的思路，

　　　　　　用一个数组$c[i]$记录第$i$数出现的次数

　　　　　　转移方程式变为：$ f[i][j*a[k]\%mod]=\sum\limits f[i-1][j]*c[k] $

　　　　　　时间复杂度变为$ O(m*mod^2) $,期望得分50pts。

　　　　对于$ 100\% $的数据，我们考虑优化，

　　　　　　方法一：矩阵快速幂+原根优化（然而蒟蒻并不会）附上wd大神博客

　　　　　　方法二：倍增优化，

　　　　　　　　回过头来再看刚才的DP式子，我们可得：$ f[i+q][j*a[k]\%mod]=f[i][j]*f[q][k] $

　　　　　　　　同样可得：$ f[2*i][j*k\%mod]=f[i][j]*f[i][k] $

　　　　　　　　于是我们可以在$ O(logm) $的时间求出 $f[2^1],f[2^2],f[2^3]\cdots,f[2^n]  (2^n<= m< 2^{n+1}) $

　　　　　　　　所以我们可以使用快速幂的思想，将m二进制拆分，下附普通快速幂模板：

ll qpow(ll x,ll y,ll mod)
{
	ll ans=1;
	while(y){
		if(y&1)ans=ans*x%mod;
		x=x*x%mod;y>>=1;
	}
	return ans;
}

　　　　　　　　将其中y变为m，x变为f数组，ans为g数组。

　　　　　　　　$f[i][j]$数组含义也变为操作$2^i$次，结果为$j$的方案数。

　　　　　　　　而$g[i][j]$数组为操作$2^1+2^2+2^3+\cdots+2^i$次，结果为$j$的方案数。

　　　　　　　　所以$ f[i][j*k\%mod]=\sum\limits_{j=1}^{mod-1}\sum\limits_{k=1}^{mod-1}f[i-1][j]*f[i-1][k] $

　　　　　　　　　　$ g[i][j*k\%mod]=\sum\limits_{j=1}^{mod-1}\sum\limits_{k=1}^{mod-1}g[i-1][j]*f[i][k] $

　　　　　　　　还有最后一个问题：内存

　　　　　　　　当前我们的空间复杂度为：$ O(logm*mod) $，而f，g数组都只与上一个状态有关，显然可用滚动数组,最后的空间复杂度为$O(2*mod)$，时间复杂度为$O(mod^2logm)$。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define R register
#define ll long long
inline ll read(){
    ll aa=0;R int bb=1;char cc=getchar();
    while(cc<'0'||cc>'9')
        {if(cc=='-')bb=-1;cc=getchar();}
    while(cc>='0'&&cc<='9')
        {aa=(aa<<1)+(aa<<3)+(cc^48);cc=getchar();}
    return aa*bb;
}
const int p=1e9+7;
const int N=1005;
ll qpow(ll x,ll y,ll mod)
{
    ll ans=1;
    while(y){
        if(y&1)ans=ans*x%mod;
        x=x*x%mod;y>>=1;
    }
    return ans;
}
ll n,mod,m,tmp=0,cur=0;
ll ans,f[2][N],g[2][N];
inline void split(int y)
{
    while(y){
        if(y&1){
            memset(g[cur^1],0,sizeof(g[cur^1]));
            for(R int i=1;i<mod;++i)
                for(R int j=1;j<mod;++j)
                    g[cur^1][i*j%mod]=(g[cur^1][i*j%mod]+f[tmp][i]*g[cur][j]%p)%p;
            cur^=1;
        }
        y>>=1;
        memset(f[tmp^1],0,sizeof(f[tmp]));
        for(R int i=1;i<mod;++i){
            for(R int j=1;j<mod;++j)
                f[tmp^1][i*j%mod]=(f[tmp^1][i*j%mod]+f[tmp][i]*f[tmp][j])%p;
        }
        tmp^=1;
    }
}
int main()
{
    n=read();m=read();mod=read();
    if(mod==2){puts("1");return 0;}
    if(n==1){
        int x=read();
        ans=qpow(x,m,mod);
        printf("%lld\n",ans%p);
        return 0;
    }
    for(R int i=1;i<=n;++i)
        ++f[0][read()];
    g[0][1]=1;
    split(m);
    for(R int i=1;i<=mod;++i)
        ans=(ans+g[cur][i]*i%p)%p;
    ans=1ll*ans*qpow(n,m*(p-2),p)%p;
    printf("%lld\n",ans);
    return 0;
}