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

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中

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

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstdlib>
#include <ctime>
using namespace std;
typedef long long LL;
const int N=1000010;
LL n,k,p,vis[N],f[N*6],ni[N],len[N];
bool ext[N];
void exgcd(LL a,LL b,LL&x,LL &y,LL &d)
{
    if(b==0)d=a,x=1,y=0;
    else exgcd(b,a%b,y,x,d),y-=x*(a/b);
}
inline LL inv(LL a,LL mod)
{
    LL d,x,y;exgcd(a,mod,x,y,d);
    return d==1?(x+mod)%mod:-1;
}
struct marx
{
    LL m[4][4];
    marx(){}
    inline LL *operator [](LL x){return m[x];}
    inline void clear(){memset(m,0,sizeof(m));}
    marx operator * (const marx &b) const
    {
        marx c;c.clear();
        for(int k=1;k<=3;k++)
            for(int i=1;i<=3;i++)
                for(int j=1;j<=3;j++)
                    c.m[i][j]=(c.m[i][j]%p+m[i][k]%p*b.m[k][j]%p)%p;
        return c;
    }
}mat[N],A,B,C,ans;
inline marx poww(marx x,LL len)
{
    marx ret;ret.clear();
    ret[1][1]=ret[2][2]=ret[3][3]=1;
    while(len)
    {
        if(len&1)ret=ret*x;
        len>>=1;x=x*x;
    }
    return ret;
}
int main()
{
    scanf("%lld%lld%lld",&n,&k,&p);
    f[1]=f[2]=1;
    for(int i=3;;i++)
    {
        f[i]=(f[i-1]+f[i-2])%k;
        if(!vis[f[i]])vis[f[i]]=i;
        if(f[i]==f[i-1]&&f[i]==1)break;
    }
    A[1][1]=A[1][2]=A[2][1]=A[3][3]=1;
    B[1][1]=B[2][2]=B[3][3]=1,B[3][1]=-1;
    ans[1][2]=ans[1][3]=1;
    LL x=1;bool flag=0;
    while(n)
    {
        if(!ni[x])ni[x]=inv(x,k);
        if(ni[x]==-1)ans=ans*poww(A,n),n=0;
        else
        {
            if(!ext[x]||flag)
            {
                ext[x]=1;
                if(!vis[ni[x]])
                    ans=ans*poww(A,n),n=0;
                else
                {
                    len[x]=vis[ni[x]];
                    if(n>=len[x])
                    {
                        n-=len[x];
                        mat[x]=poww(A,len[x])*B;
                        ans=ans*mat[x],x=x*f[len[x]-1]%k;
                    }
                    else ans=ans*poww(A,n),n=0;
                }
            }
            else
            {
                LL cnt=0;
                C.clear();C[1][1]=C[2][2]=C[3][3]=1;
                for(LL i=x*f[len[x]-1]%k;i!=x;i=i*f[len[i]-1]%k)
                    cnt+=len[i],C=C*mat[i];
                cnt+=len[x],C=mat[x]*C;
                ans=ans*poww(C,n/cnt);
                n%=cnt,flag=1;
            }
        }
    }
    printf("%lld",(ans[1][1]%p+p)%p);
}