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 }