zoukankan      html  css  js  c++  java
  • 随(rand):原根,循环矩阵,dp

    20分特判,一个puts("1")一个快速幂,不讲。

    50%算法:

    上次就讲了,可是应该还是有像 xuefen某 或 Dybal某 一样没听的。

    用a×inv(b)%mod来表示分数的时候,这个分数值可加可乘(有空证明)

    像是一个dp题啊。

    初状态是1方案数为1,然后做乘法转移不就好了嘛?

    设dp[i][j]表示进行了i次操作后所得的值为j

    dp[i][j*a[k]%mod]+=dp[i-1][j];

    复杂度O(mod2×m)

     1 #include<iostream>
     2 #include<cstdio>
     3 #include<cmath>
     4 #include<algorithm>
     5 #include<queue>
     6 #include<vector>
     7 #include<string>
     8 #include<cstring>
     9 #define int long long
    10 #define m(a) memset(a,0,sizeof(a))
    11 #define AA cout<<"Alita"<<endl
    12 using namespace std;
    13 const int mod=1e9+7;
    14 int ans,S,n,m,p,v[1050],a[100500],f[1050][1050],g[1050][1050],tmp[1050][1050];
    15 int poww(int x,int y,int z)
    16 {
    17         int sum=1;
    18         while(y)
    19         {
    20                 if(y&1) sum=sum*x%z;
    21                 y>>=1;
    22                 x=x*x%z;
    23         }
    24         return sum;
    25 }
    26 void X_g()
    27 {
    28         memset(tmp,0,sizeof(tmp));
    29         for(int i=1;i<p;i++)
    30         {
    31                 for(int j=1;j<p;j++)
    32                 {
    33                         (tmp[1][i]+=g[1][j]*f[j][i]%mod)%=mod;
    34                 }
    35         }
    36         for(int i=1;i<p;i++) 
    37         {
    38                 g[1][i]=tmp[1][i];
    39         }
    40 }
    41 void X_f()
    42 {
    43         memset(tmp,0,sizeof(tmp));
    44         for(int i=1;i<p;i++)
    45         {
    46                 for(int j=1;j<p;j++)
    47                 {
    48                         for(int k=1;k<p;k++)
    49                         {
    50                                 (tmp[i][j]+=f[i][k]*f[k][j]%mod)%=mod;
    51                         }
    52                 }
    53         }
    54         for(int i=1;i<p;i++)
    55         {
    56                 for(int j=1;j<p;j++)
    57                 {
    58                         f[i][j]=tmp[i][j];
    59                 }
    60         }
    61 }
    62 void work(int y)
    63 {
    64         while(y)
    65         {
    66                 if(y&1) X_g();
    67                 y>>=1;
    68                 X_f();
    69         }
    70 }
    71 signed main()
    72 {
    73         //freopen("1.in","r",stdin);
    74         //freopen("1.out","w",stdout);
    75         scanf("%lld%lld%lld",&n,&m,&p);
    76         if(p==2){puts("1");return 0;}
    77         for(int i=1;i<=n;i++) scanf("%lld",&a[i]);
    78         if(n==1){printf("%lld",poww(a[1],m,p));return 0;}
    79         S=poww(n,mod-2,mod);
    80         for(int i=1;i<=n;i++)
    81         {
    82                 (v[a[i]%p]+=S)%=mod;
    83         }
    84         g[1][1]=1;
    85         for(int j=1;j<p;j++)
    86         {
    87                 for(int k=1;k<p;k++)
    88                 {
    89                         f[j][j*k%p]+=v[k];
    90                 }
    91         }
    92         work(m);
    93         for(int i=1;i<p;i++)
    94         {
    95                 (ans+=g[1][i]*i%mod)%=mod;
    96         }
    97         printf("%lld",ans);
    98         return 0;
    99 }
    Dybala的50分代码(已征得同意后转载)

    因为我直接没想这个这么暴力的dp,看到题目的m已经破1e8了那么要么复杂度与之无关要么把它log掉

    理论80%算法:

    这题m又不是一个计算参数,复杂度不太可能与之无关,尝试log?

    带log的dp。。。矩阵快速幂啊!

    可以发现dp的第二维并不大,是mod级别,矩阵乘mod3貌似勉强可以接受

    O(mod3×log2m)的复杂度。我也不知道题解为什么说能得80分。

    略微一算,知道复杂度已经略微超过极限,卡常?没有用,还是50分。

    个人认为把矩阵快速幂和暴力dp的得分设置成一样不太道德。

    到现在还有人不会用矩阵快速幂优化dp啊。。。这。。。我该怎么讲?

    ans矩阵用来更新答案,base矩阵是快速幂的基底。

    普通快速幂是这个样子的:for(;t;t>>=1,base=base*base){if(t&1)ans*=base;base*=base;}

    换成矩阵也是一样的啊:for(;t;t>>=1){if(t&1)mult_ans();mult_base();}

    只不过是两个矩阵乘法函数而已。

    提醒:包括用重载运算符的矩阵乘,传参时要带上‘&’符号,不然有时会TLE/RE

    base矩阵的构造也很简单,对于每一个可能被选到的数a[i],枚举乘之前的数j。

    base[j][j*a[i]%p]+=1;(这里是求解总方案书数,如果是概率打法要把1改成概率)

    就是能从a转移到b的话,base[a][b]就有值,是方案的话就是1,算概率的话就是概率。

    ans矩阵可以弄成一维数组,能略快一些。这题初始为1,那么ans[1]=1就好了。

    注意:代码中的p是原题中的mod,而代码中的mod是1e9+7。我是概率打法

     1 #include<cstdio>
     2 #define int long long
     3 #define mod 1000000007ll
     4 int n,m,p,base[1003][1003],ans[1003],res[1003][1003],a[100005],ANS;
     5 int pow(int bbase,int t,int modd,int aans=1){
     6     for(;t;t>>=1,bbase=bbase*bbase%modd)if(t&1)aans=aans*bbase%modd;
     7     return aans;
     8 }
     9 void mult_base(){
    10     for(int i=0;i<p;++i)for(int j=0;j<p;++j)for(int k=0;k<p;++k)(res[i][j]+=base[i][k]*base[k][j])%=mod;
    11     for(int i=0;i<p;++i)for(int j=0;j<p;++j)base[i][j]=res[i][j],res[i][j]=0;
    12 }
    13 void mult_ans(){
    14     for(int i=0;i<p;++i)for(int j=0;j<p;++j)(res[0][j]+=ans[i]*base[i][j])%=mod;
    15     for(int i=0;i<p;++i)ans[i]=res[0][i],res[0][i]=0;
    16 }
    17 signed main(){
    18     scanf("%lld%lld%lld",&n,&m,&p);
    19     if(n==1){
    20         scanf("%lld",&a[1]);
    21         printf("%lld
    ",pow(a[1],m,p));return 0;
    22     }
    23     ans[1]=1;const int P=pow(n,mod-2,mod);//概率,即1/n
    24     for(int i=1;i<=n;++i)scanf("%lld",&a[i]),a[i]%=p;
    25     for(int i=1;i<=n;++i)for(int j=0;j<p;++j)(base[j][j*a[i]%p]+=P)%=mod;
    26     for(;m;m>>=1){if(m&1)mult_ans();mult_base();}
    27     for(int i=0;i<p;++i)(ANS+=ans[i]*i)%=mod;
    28     printf("%lld
    ",ANS);
    29 }
    卡死常数还是50分

    100%算法:

    想象一下如果题目改一下,你是随机选择数加上它而不是乘,该怎么做?

    原始dp方程:dp[i][(j+a[k])%mod]+=dp[i-1][j];

    同理,构造矩阵。

    base[i][(i+a[j])%mod]=1;//或者概率

    写几个不太大的样例,把矩阵写出来,你会发现,它是一个循环矩阵。

    这个加法类型的转移矩阵基本上都是循环矩阵(对于不同的被转移数,加数都一样)

    回到这道题。

    首先从复杂度上考虑,也许可以猜出来这是个循环矩阵。

    可是那个诡异的乘法形式并不是循环的,虽然它对于不同的被转移数,乘数都一样。

    如果它也是一个加法,该多好啊。

    没办法,它是个乘法,我们要接受这个现实。

    但是,这并不能阻碍我们把它强制弄成加法。

    看看孙金宁老师的叮嘱:啊,什么原根,啊,几次方,啊什么取到所有值。。。(困)

    次方?乘法?加法?欧拉定理?

    xa × xb=xa+b

    诶,这里有加法了!这样可以把乘法换成加法。

    动用一下欧拉定理:

    xa × xb=x(a+b)%(mod-1)

    如果题目的所有输入都是x的几次方,就很好做了。

    考虑:刚开始一个数是a0,你可以从a4,a7,a23里面随机选数把它乘起来,值对mod取模,得到ans,求最后ans的期望

    这个问题就相当于:考虑:刚开始一个数是0,你可以从4,7,23里面随机选数把它加起来,值对mod-1取模,求最后aans%mod的期望

    好做好做!就是普通的加法dp,循环矩阵肝它!

    但是。。。这个a是多少呢?

    继续听孙金宁的数学课。原根啥啥啥的。。。

    原根?啥?原根?哦。好吧。原根就原根吧。

    如果你比我聪明,你可能会直接发现,既然原根的次方值能够取遍1~mod-1的所有数,那么这些数就都可以用原根唯一确定的表示出来。那么就可以把题目的所有数都用原根表示,然后当成加法dp来做。你就A了。

    如果你没我聪明,那么就相当与你和我一样聪明。

    那么如果咱们一样聪明,诶,那好啊,咱们都想不到。既然它一直在念叨原根,那就试试拿原根表示呗。

    那么首先我们需要求出原根。用题目给出的那个充要条件很好求。

    for(int i=1;i<p;i++){
        int now=1,j;
        for(j=0;j<p-1;++j)
            if(al[now]==-1)//像它描述里说的一样,不能出现重复的,因为它要在mod-1次里取到mod-1个不同值,故不能重复
                al[now]=j,//记录下now这个值唯一确定的对应着i的几次幂
                qpow[j]=now,//记录下来i的j次幂是几,以后会用到
                now=now*i%p;//now是i的j次幂,j++,更新
            else break;
        if(j==p-1){g=i;break;}//如果j成功的走完了,没有被中途break掉,那么i就是原根。
        else for(int i=1;i<p;++i)al[i]=-1;//清空
    }

    然而我这么求是用了它说的第一条,这样的话可以顺便求出qpow和al数组里面的值。

    al数组不仅是用来标记是否出现过的。如果它出现过,还记录了它对应的是原根g的几次幂。

    接下来,我们就拥有了1~mod-1里面每个值和g的几次幂间的双向映射。

    那么就把原问题映射成加法dp,循环矩阵干他,再映射回来统计答案即可。

     1 #include<cstdio>
     2 #define int long long
     3 #define mod 1000000007ll
     4 int n,m,p,base[1003],ans[1003],res[1003],a[100005],ANS;
     5 int al[1003],g,qpow[1003];
     6 int pow(int bbase,int t,int modd,int aans=1){
     7     for(;t;t>>=1,bbase=bbase*bbase%modd)if(t&1)aans=aans*bbase%modd;
     8     return aans;
     9 }
    10 void mult_base(){
    11     for(int i=0;i<p;++i)for(int j=0;j<p;++j)(res[j]+=base[i]*base[(j-i+p)%p])%=mod;
    12     for(int i=0;i<p;++i)base[i]=res[i],res[i]=0;
    13 }
    14 void mult_ans(){
    15     for(int i=0;i<p;++i)for(int j=0;j<p;++j)(res[j]+=ans[i]*base[(j-i+p)%p])%=mod;
    16     for(int i=0;i<p;++i)ans[i]=res[i],res[i]=0;
    17 }
    18 signed main(){
    19     scanf("%lld%lld%lld",&n,&m,&p);
    20     ans[0]=1;const int P=pow(n,mod-2,mod);
    21     for(int i=1;i<=1000;++i)al[i]=-1;
    22     for(int i=1;i<p;i++){
    23         int now=1,j;
    24         for(j=0;j<p-1;++j)
    25             if(al[now]==-1)al[now]=j,qpow[j]=now,now=now*i%p;
    26             else break;
    27         if(j==p-1){g=i;break;}
    28         else for(int i=1;i<p;++i)al[i]=-1;
    29     }p--;
    30     for(int i=1;i<=n;++i)scanf("%lld",&a[i]),a[i]=al[a[i]];
    31     for(int i=1;i<=n;++i)(base[a[i]]+=P)%=mod;
    32     for(;m;m>>=1){if(m&1)mult_ans();mult_base();}
    33     for(int i=0;i<p;++i)(ANS+=ans[i]*qpow[i])%=mod;
    34     printf("%lld
    ",ANS);
    35 }
    于倬浩:然后你就愉快的拿到了100分的好成绩
  • 相关阅读:
    传统金融和互联网金融
    集团培训
    Javascript和JQuery之间的联系
    this和$(this)区别
    原生JavaScript支持6种方式获取元素
    绩效考核
    web服务端安全之分布式拒绝服务攻击
    web服务端安全之暴力破解
    web服务端安全之权限漏洞
    web服务端安全之文件上传漏洞
  • 原文地址:https://www.cnblogs.com/hzoi-DeepinC/p/11256096.html
Copyright © 2011-2022 走看看