[BZOJ2432][Noi2011]兔农 矩阵乘法+exgcd

Posted 2020-10-01 LadyLex

2432: [Noi2011]兔农

Time Limit: 10 Sec  Memory Limit: 256 MB

Description

 农夫栋栋近年收入不景气，正在他发愁如何能多赚点钱时，他听到隔壁的小朋友在讨论兔子繁殖的问题。
问题是这样的：第一个月初有一对刚出生的小兔子，经过两个月长大后，这对兔子从第三个月开始，每个月初生一对小兔子。新出生的小兔子生长两个月后又能每个月生出一对小兔子。问第n个月有多少只兔子？
聪明的你可能已经发现，第n个月的兔子数正好是第n个Fibonacci(斐波那契)数。栋栋不懂什么是Fibonacci数，但他也发现了规律：第i+2个月的兔子数等于第i个月的兔子数加上第i+1个月的兔子数。前几个月的兔子数依次为：
1 1 2 3 5 8 13 21 34 …
栋栋发现越到后面兔子数增长的越快，期待养兔子一定能赚大钱，于是栋栋在第一个月初买了一对小兔子开始饲养。
每天，栋栋都要给兔子们喂食，兔子们吃食时非常特别，总是每k对兔子围成一圈，最后剩下的不足k对的围成一圈，由于兔子特别害怕孤独，从第三个月开始，如果吃食时围成某一个圈的只有一对兔子，这对兔子就会很快死掉。
我们假设死去的总是刚出生的兔子，那么每个月的兔子数仍然是可以计算的。例如，当k=7时，前几个月的兔子数依次为：
1 1 2 3 5 7 12 19 31 49 80 …
给定n，你能帮助栋栋计算第n个月他有多少对兔子么？由于答案可能非常大，你只需要告诉栋栋第n个月的兔子对数除p的余数即可。

Input

输入一行，包含三个正整数n, k, p。

Output

输出一行，包含一个整数，表示栋栋第n个月的兔子对数除p的余数。

Sample Input

6 7 100

Sample Output

7

HINT

1<=N<=10^18

2<=K<=10^6

2<=P<=10^9

 

题解：

这题怎么说呢……又爱又恨……

暴力分给的特别足，有75分，作为noi题目简直太良心了！

可是，如果想拿到剩下的25分特别困难……思路必须十分清晰

在5分钟打完暴力分之后，我们不妨打个表观察一下。

 拿样例的%7作为例子

按照题意计算的数列，%7后大概长这样（为了美观，我们每出现一个0换一次行）：

1, 1, 2, 3, 5, 0,
5, 5, 3, 0,
3, 3, 6, 2, 0,
2, 2, 4, 6, 3, 2, 5, 0,
5, 5, 3, 0,
3, 3, 6, 2, 0,

……

那么这又有什么规律可循呢？

性质1：我们发现，每段开头必为相同的两数，并且它们恰是上一段的最末一位非0数；由于总共只有k?1种余数，

所以不超过k段就会出现循环（如果有的话），比如上面k=7时的上面的数列，

5, 5, 3, 0,
3, 3, 6, 2, 0,
2, 2, 4, 6, 3, 2, 5, 0,

就是一个循环节。

然后我们可以发现一些新的性质：

性质2：设开头的两个数为a，那么由于这是斐波那契数列，后面的数会变成2*a%k,3*a%k,5*a%k……

这其实很容易证明。如果上一行出现了“……a 0 ”，那么这一行第一位是a+0=a，第二位是0+a=a，后面一样递推

第i项的系数，就是斐波那契数列第i项的值。

那么，在我们减1使最后一个余数变为0之前，应该有：

a* f[len] ≡ 1 （mod k）

不难发现，f[len]就是a在mod k意义下的逆元

所以，我们可以通过exgcd或者扩展欧拉定理，来快速求出f[len]。

如果这个逆元不存在，就证明后面的操作中没有循环节了，我们直接矩阵快速幂递推即可

由上面的性质1可以看出，下一行的第一个值b=a* f[len-1]%k

因此我们还需要知道len是多少

那么我们再定义一个vis数组，其中vis[i]表示f值等于i的len最小值（有点拗口……）

那么显然，对于某一个f[len]，vis[f[len]]=len（循环节中肯定是最短的啊……）

有了vis数组以后，我们就可以快速由f[len]求出len值了。

这样，我们的处理流程就是：

(1)根据a* f[len] ≡ 1 （mod k）,求出f[len]
(2)根据f[len]和vis数组反推出len
(3)由x*f[len-1]得到下一段的开头。

现在我们的问题就变成了如何求vis数组？

 然后据说有一个很强的结论：斐波那契数列在模k意义下一定是以0 1 1……为开头的循环，并且循环节长度<=6*k（但也可能很长，比n大），因此不用担心算不完，只需要暴力递推记录每一位%k的f值和vis值即可

还剩下一点就是转移矩阵：在行内使用下面的A矩阵，在行间转移使用B矩阵。答案矩阵一开始长C这样

1 1 0    1 0 0  0 1 1

1 0 0    0 1 0

0 0 1   -1 0 1

   A           B　　C

这样，我们的处理过程就是：递推记录vis数组和f数组->建立转移矩阵->计算每一行的转移矩阵->如果出现循环节，用矩阵乘法快速处理->将剩余的转移用A矩阵乘到C中

这样，本题就被我们切掉了……真的是很不错的一道题。具体实现时还有一些小点需要注意。代码见下：

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #include <cstdlib>
 5 #include <ctime>
 6 using namespace std;
 7 typedef long long LL;
 8 const int N=1000010;
 9 LL n,k,p,vis[N],f[N*6],ni[N],len[N];
10 bool ext[N];
11 void exgcd(LL a,LL b,LL&x,LL &y,LL &d)
12 {
13     if(b==0)d=a,x=1,y=0;
14     else exgcd(b,a%b,y,x,d),y-=x*(a/b);
15 }
16 inline LL inv(LL a,LL mod)
17 {
18     LL d,x,y;exgcd(a,mod,x,y,d);
19     return d==1?(x+mod)%mod:-1;
20 }
21 struct marx
22 {
23     LL m[4][4];
24     marx(){}
25     inline LL *operator [](LL x){return m[x];}
26     inline void clear(){memset(m,0,sizeof(m));}
27     marx operator * (const marx &b) const
28     {
29         marx c;c.clear();
30         for(int k=1;k<=3;k++)
31             for(int i=1;i<=3;i++)
32                 for(int j=1;j<=3;j++)
33                     c.m[i][j]=(c.m[i][j]%p+m[i][k]%p*b.m[k][j]%p)%p;
34         return c;
35     }
36 }mat[N],A,B,C,ans;
37 inline marx poww(marx x,LL len)
38 {
39     marx ret;ret.clear();
40     ret[1][1]=ret[2][2]=ret[3][3]=1;
41     while(len)
42     {
43         if(len&1)ret=ret*x;
44         len>>=1;x=x*x;
45     }
46     return ret;
47 }
48 int main()
49 {
50     scanf("%lld%lld%lld",&n,&k,&p);
51     f[1]=f[2]=1;
52     for(int i=3;;i++)
53     {
54         f[i]=(f[i-1]+f[i-2])%k;
55         if(!vis[f[i]])vis[f[i]]=i;
56         if(f[i]==f[i-1]&&f[i]==1)break;
57     }
58     A[1][1]=A[1][2]=A[2][1]=A[3][3]=1;
59     B[1][1]=B[2][2]=B[3][3]=1,B[3][1]=-1;
60     ans[1][2]=ans[1][3]=1;
61     LL x=1;bool flag=0;
62     while(n)
63     {
64         if(!ni[x])ni[x]=inv(x,k);
65         if(ni[x]==-1)ans=ans*poww(A,n),n=0;
66         else
67         {
68             if(!ext[x]||flag)
69             {
70                 ext[x]=1;
71                 if(!vis[ni[x]])
72                     ans=ans*poww(A,n),n=0;
73                 else
74                 {
75                     len[x]=vis[ni[x]];
76                     if(n>=len[x])
77                     {
78                         n-=len[x];
79                         mat[x]=poww(A,len[x])*B;
80                         ans=ans*mat[x],x=x*f[len[x]-1]%k;
81                     }
82                     else ans=ans*poww(A,n),n=0;
83                 }
84             }
85             else
86             {
87                 LL cnt=0;
88                 C.clear();C[1][1]=C[2][2]=C[3][3]=1;
89                 for(LL i=x*f[len[x]-1]%k;i!=x;i=i*f[len[i]-1]%k)
90                     cnt+=len[i],C=C*mat[i];
91                 cnt+=len[x],C=mat[x]*C;
92                 ans=ans*poww(C,n/cnt);
93                 n%=cnt,flag=1;
94             }
95         }
96     }
97     printf("%lld",(ans[1][1]%p+p)%p);
98 }