【矩阵快速幂】[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

分析

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

FibiFibi+1Fibi+2Fibi+3Fibi+4Fibi+5=0=Fibi1=Fibi1=2Fibi1=3Fibi1=5Fibi1

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

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

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define MOD p
#define MAXK 1000000
typedef 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* a
Matrix 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);
}
posted @ 2016-07-13 09:26  outer_form  阅读(545)  评论(0编辑  收藏  举报