数学 随rand

问题 A: 随
时间限制: 2 Sec 内存限制: 512 MB
题目描述
给出n个正整数a1,a2…an和一个质数mod.一个变量x初始为1.进行m次操作.每次在n个数中随机选一个ai,然后x=x*ai%mod.问m次操作之后x的取值的期望.
答案一定可以表示成a/b的精确分数形式.a和b可能很大,所以只需要输出a*(b^(10^9+5))模10^9+7的结果.
输入
第一行三个整数n,m,mod.
接下来一行n个空格隔开的正整数a1,a2…an
输出
一行一个整数表示答案
样例输入
2
1 3 1 2
样例输出
3
提示
第1个测试点:mod=2
第2个测试点:n=1
第3,4,5个测试点:m<=1000,1<=mod<=300.
第6,7,8个测试点: 1<=mod<=300
对于全部测试点: 1<=ai <=mod,mod为质数1<=mod<=1000,
对于全部测试点:1<=n<=10^5,1<=m<=10^9
【孙金宁教你学数学】
质数P的原根g满足1<=rt<=P,且rt的1次方,2次方…(P-1)次方在模P意义下可以取遍1到(P-1)的所有整数.
欧拉定理:对于质数P,1<=x<=P的任意x的P-1次方在模P意义下都为1.
显然,原根的1次方,2次方…（P-2）次方在模P意义下都不为1,只有(P-1)次方在模P意义下为1.这也是一个数成为原根的充分必要条件.

蛮有意思，尤其是这个原根。如果这个数的1～p-2次方%p全不等一，那他就是原根了。所以利用这个，只要考虑x乘了i次之后是原根的j次方的概率，当幂的次数大于p-1时，就从头来（g^p-1%p==1）知道了是原根的j次方，那一定对应不同的数，统计出来每个数出现的概率（最多有一千个数，大大削减了n=10^5）就方便多了。
转移需要矩阵乘+滚动数组。
f[i][(j+k)%(p-1)]+=f[i-1][j]*f[i-1][k];
在这里我学了学倍增法。预处理出乘2^i时为j次幂的概率，就可以在矩阵乘使用了。
原先快速幂要把一个数不断平方，而倍增因为需要乘的数较多，所以预处理好了。

#include<cstdio>#include<cstdlib>#include<iostream>#include<algorithm>#include<cstring>#define Mod 1000000007#define ll long longusing namespace std;ll n,m,mod1,mod2,g,b=1;ll p[34][1005],f[2][1005],xp[1005],a[1005],cnt[1005];ll get(ll x){    for(int i=1;i<=x;i++)    {        ll p=1;bool flag=1;        for(int j=1;j<x-1;j++)        {            p=p*i%mod1;            if(p==1){flag=0;break;}        }        if(flag)return i;     }}int main(){    //freopen("rand.in","r",stdin);    //freopen("rand.out","w",stdout);    scanf("%lld%lld%lld",&n,&m,&mod1);mod2=mod1-1;    for(ll x=Mod-2,y=n;x;x/=2,y=y*y%Mod)if(x&1)b=b*y%Mod;    g=get(mod1);xp[0]=1;int x;    for(int i=1;i<=mod1;i++)xp[i]=xp[i-1]*g%mod1;    for(int i=0;i<mod2;i++)a[xp[i]]=i;    for(int i=1;i<=n;i++)scanf("%d",&x),cnt[a[x]]++;    for(int i=0;i<mod2;i++)p[0][i]=b*cnt[i]%Mod;    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-1][j]*p[i-1][k]%Mod)%=Mod;    f[0][0]=1;    for(int k=0;m;m/=2,k++)    {        if(m&1)        {            for(int i=0;i<mod2;i++)f[1][i]=0;            for(int i=0;i<mod2;i++)                for(int j=0;j<mod2;j++)                   (f[1][(i+j)%mod2]+=f[0][i]*p[k][j]%Mod)%=Mod;             for(int i=0;i<mod2;i++)f[0][i]=f[1][i];        }    }       ll ans=0;       for(int i=0;i<mod2;i++)(ans+=f[0][i]*xp[i]%Mod)%=Mod;    cout<<ans;}