【矩阵快速幂】[NOI2011]兔农

题目描述

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

Source

Day1

分析

有一个结论，模P意义下Fibonacci数列的循环节≤6P
这个结论让我们可以暴力求一次模k意义下的循环节。
再来看题目，当第一次遇到Fibi≡1(modk)，Fibi−−。
那我们来看看减一之后在模k意义下的这个数列。

FibiFibi+1Fibi+2Fibi+3Fibi+4Fibi+5=0=Fibi−1=Fibi−1=2∗Fibi−1=3∗Fibi−1=5∗Fibi−1⋮

我们先这又是一个Fibonacci数列，只不过乘了Fibi−1，那么，当 Fibonacci中第一次出现Fibi−1模k意义下的乘法逆元的时候，将其乘 Fibi−1模k就等于1。
我们可以先预处理出乘法逆元，然后在求循环节的时候，假设在Fibi遇到b的乘法逆元，也就是Fibi的乘法逆元是b，lenb=i,nextb=b∗Fibi−1modk，最后再看next数组中存不存在循环，如果存在，就先求一个循环（循环不一定包括1），再快速幂，如果没有，就直接快速幂。

矩阵不能像一般的Fibonacci数列那样构造，因为矩阵加减和乘法放在一起是没有结合律的，而如果有循环而且不包括1就会交换乘法与减法的运算顺序，所以我们只能有乘法。将矩阵再增加一阶就好了。

代码

#include<cstdio>#include<cstring>#include<algorithm>using namespace std;#define MOD p#define MAXK 1000000typedef long long LL;int k,p,a[MAXK*10+10],cir,c[MAXK+10],nxt[MAXK+10],ans,inv[MAXK+10],bgc;bool vis[MAXK+10];LL n,dep[MAXK+10];template<class T>void Read(T &x){    char c;    while(c=getchar(),c!=EOF)        if(c>='0'&&c<='9'){            x=c-'0';            while(c=getchar(),c>='0'&&c<='9')                x=x*10+c-'0';            ungetc(c,stdin);            return;        }}void exgcd(int a,int b,int &d,int &x,int &y){    if(!b){        d=a;        x=1;        y=0;        return;    }    exgcd(b,a%b,d,y,x);    y-=a/b*x;}// gcd (a,b)=gcd(b,a-a/b*b)//ax+by=bx'+(a-a/b*b)y'//ax+by=bx'+ay'-a/b*by'=ay'+b(x'-a/b*y')//x=y' y=x'-a/b*y'void prepare(){    int i,d,x,y;    inv[1]=1;    for(i=2;i<k;i++){        exgcd(i,k,d,x,y);        if(d==1)            inv[i]=((x%k)+k)%k;    }    a[0]=a[1]=1;    for(i=2;;i++){        a[i]=(a[i-1]+a[i-2])%k;        if(a[i]&&inv[a[i]]&&!c[inv[a[i]]])            c[inv[a[i]]]=i+1,nxt[inv[a[i]]]=1ll*a[i-1]*inv[a[i]]%k;        if(!a[i-1]&&a[i]==1)            break;    }    cir=i;}struct Matrix{    int a[3][3];    inline Matrix(){        memset(a,0,sizeof a);    }    inline Matrix(int){        memset(a,0,sizeof a);        a[0][0]=a[1][1]=a[2][2]=1;    }    inline Matrix operator * (const Matrix &b)const{        Matrix c;        static int i,j,k;        for(i=0;i<3;i++)            for(j=0;j<3;j++)                if(a[i][j])                    for(k=0;k<3;k++)                        c.a[i][k]=(c.a[i][k]+1ll*a[i][j]*b.a[j][k])%MOD;        return c;    }    inline Matrix& operator *=(const Matrix &b){        return *this=*this*b;    }}x,y(1),z(1),xx;//f[0] f[-1]  1* aMatrix quick_pow(Matrix a,LL b){    Matrix ret(1);    while(b){        if(b&1)            ret*=a;        a*=a;        b>>=1;    }    return ret;}void solve(){    int b;    LL tot=0;    b=1;    x.a[0][0]=x.a[1][0]=x.a[0][1]=x.a[2][2]=1;    xx=x,xx.a[2][0]=p-1;    while(!vis[b]&&nxt[b]&&tot+c[b]<=n){        vis[b]=1;        dep[b]=tot;        tot+=c[b];        y*=quick_pow(x,c[b]-1);        y*=xx;        b=nxt[b];    }    if(nxt[b]&&tot+c[b]<=n){        bgc=b;        do{            z*=quick_pow(x,c[b]-1);            z*=xx;            b=nxt[b];        }while(b!=bgc);        y*=quick_pow(z,(n-tot)/(tot-dep[b]));        tot+=(n-tot)/(tot-dep[b])*(tot-dep[b]);        while(tot+c[b]<=n){            tot+=c[b];            y*=quick_pow(x,c[b]-1);            y*=xx;            b=nxt[b];        }    }    y*=quick_pow(x,n-tot);    ans=(y.a[1][0]+y.a[2][0])%MOD;}int main(){    Read(n),Read(k),Read(p);    prepare();    solve();    printf("%d\n",ans);}