[bzoj2432][矩阵乘法]兔农

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

题解

这个真心可爱的一逼..
对我这种懒得打表的人就是不和谐
手玩一下样例膜K
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
….
显然有一个性质
分段之后每段开头一定是两个相同的数
借着这个可以推出第二个性质
循环节不会超过K段（如果有循环节的话）
设上一段结尾的数是x
下一段的数可以写为
$1x,1x,2x,3x,5x,....$
发现是一个斐波那契数列
当我们减一个1使得最后一个数成0的时候，应该有
$x*t \equiv 1 (\mod K)$
发现$t$就是$x$在模意义下的逆元
预处理小于K的数的逆元在斐波那契数列中第一个出现的位置
这个时候引入一个结论
斐波那契数列在模K意义下一定是以0 1 1 开头的循环数列 且循环长度不超过6*K
暴力预处理即可
剩下就是转移了
行之间转移用A矩阵 列之间转移用B矩阵 找到循环节后直接C矩阵处理转移
代码还挺好看的..

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#define mod K
#define LL long long
using namespace std;
LL n,K,P;
struct matrix
{
    LL m[4][4];
    matrix(){memset(m,0,sizeof(m));}
    friend matrix operator *(matrix u,matrix v)
    {
        matrix ret;
        for(int i=1;i<=3;i++)
            for(int j=1;j<=3;j++)
                for(int k=1;k<=3;k++)
                    ret.m[i][k]=(ret.m[i][k]+u.m[i][j]*v.m[j][k])%P;
        return ret;
    }
}tp[1100000],ans,A,B,C;int len[1100000];
LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if(a==0)
    {
        x=0;y=1;
        return b;
    }
    else
    {
        LL tx,ty;
        LL d=exgcd(b%a,a,tx,ty);
        x=ty-(b/a)*tx;
        y=tx;
        return d;
    }
}
LL getinv(LL p)
{
    LL A=p,B=mod,x,y,K=1;
    LL d=exgcd(A,B,x,y);
    if(d!=1)return -1;
    x=(x*(K/d)%(B/d)+(B/d))%(B/d);
    return x;
}
matrix pow_mod(matrix u,LL b)
{
    matrix ret;
    for(int i=1;i<=3;i++)ret.m[i][i]=1;
    while(b)
    {
        if(b&1)ret=ret*u;
        u=u*u;b>>=1;
    }
    return ret;
}
int gg[1100000];
int vis[1100000];LL f[6100000],inv[1100000];
int main()
{
    scanf("%lld%lld%lld",&n,&K,&P);
    ans.m[1][2]=ans.m[1][3]=1;
    A.m[1][1]=A.m[1][2]=A.m[2][1]=A.m[3][3]=1;
    B.m[1][1]=B.m[2][2]=B.m[3][3]=1;B.m[3][1]=-1;
    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;
    }
    memset(inv,0,sizeof(inv));
    LL x=1;bool isring=false;//是否出现循环节 
    while(n)
    {
        if(!inv[x])inv[x]=getinv(x);
        if(inv[x]==-1){ans=ans*pow_mod(A,n);n=0;}
        else
        {
            if(!gg[x]||isring)//出现循环节或者之前没有访问过这个余数
            {
                gg[x]=1;
                if(!vis[inv[x]])//fib中没有这个数
                    ans=ans*pow_mod(A,n),n=0;
                else
                {
                    len[x]=vis[inv[x]];
                    if(n>=len[x])
                    {
                        n-=len[x];
                        tp[x]=pow_mod(A,len[x])*B;
                        ans=ans*tp[x];x=x*(LL)f[len[x]-1]%K;
                    }
                    else ans=ans*pow_mod(A,n),n=0;

                }
            }
            else//再次访问 出现循环节 
            {
                C.m[1][1]=C.m[2][2]=C.m[3][3]=1;
                LL sum=len[x];C=C*tp[x];
                for(LL i=x*(LL)f[len[x]-1]%K;i!=x;i=i*(LL)f[len[i]-1]%K)sum+=len[i],C=C*tp[i];
                ans=ans*pow_mod(C,n/sum);
                isring=true;n=n%sum;
            }
        }
    }
    printf("%lld\n",(ans.m[1][1]%P+P)%P);
    return 0;
}