[NOI2011] 兔农 矩阵乘法，矩阵的等比数列求和

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;
ll n,p;
ll k;
#define N 4000000
ll a[N],fr[N];
struct sq{
    ll a[3][3];
    sq(){memset(a,0,sizeof(a));}
    sq operator*(sq x)const{
           sq ans;
           for(ll i=1;i<=2;i++)for(ll j=1;j<=2;j++)for(ll z=1;z<=2;z++)(ans.a[i][j]+=a[i][z]*x.a[z][j])%=p;
           return ans;
    }
    void out(){
        for(ll i=1;i<=2;i++){
            for(ll j=1;j<=2;j++){
                cout<<a[i][j]<<" ";
            }
            cout<<endl;
        }
        cout<<endl;
    }
}tt;
ll hash[N];
struct sq2{
    ll a[5][5];
    sq2(){memset(a,0,sizeof(a));}
    sq2 operator*(sq2 x){
        sq2 ans;
        for(ll i=1;i<=4;i++)for(ll j=1;j<=4;j++)for(ll z=1;z<=4;z++)(ans.a[i][j]+=a[i][z]*x.a[z][j])%=p;
        return ans;
    }
    void out(){
        for(ll i=1;i<=4;i++){
            for(ll j=1;j<=4;j++){
                cout<<a[i][j]<<" ";
            }
            cout<<endl;
        }
        cout<<endl;
    }
};
sq ff(sq x,ll y){
    sq ans,k=x;
    for(ll i=1;i<=2;i++)ans.a[i][i]=1;
    for(ll i=0;(1LL<<i)<=y;i++){
        if((1LL<<i)&y)ans=ans*k;
        k=k*k;
    }
    return ans;
}
sq cheng(sq x,ll y){
    sq2 tr,ans1;
    for(ll i=1;i<=2;i++)for(ll j=1;j<=2;j++)ans1.a[i][j]=x.a[i][j];
    for(ll i=1;i<=2;i++){
        for(ll j=1;j<=2;j++){
           tr.a[i][j]=x.a[i][j];
        }
        tr.a[i][i+2]=1;
        tr.a[i+2][i+2]=1;
        ans1.a[i][i+2]=1;
    }
    for(ll i=0;(1LL<<i)<=y;i++){
        if(y&(1LL<<i)){
            ans1=ans1*tr;    
        }
        tr=tr*tr;
    }
    sq kk;
    for(ll i=1;i<=2;i++)for(ll j=1;j<=2;j++){
        kk.a[i][j]=ans1.a[i][j+2];
    }
    return kk;
}
ll gcd(ll &p,ll &q,ll x,ll y){
    if(x==0){
        p=1,q=1;
        return y;
    }
    ll r=gcd(q,p,y%x,x);
    p-=q*(y/x);
    return r;
}
ll gcd(ll x,ll y){
    return x==0?y:gcd(y%x,x);
}
ll bn;
ll bb[N],cir;
ll ben[N];
ll res;
ll ni(ll x){
    ll p,q;
    if(gcd(x,k)!=1)return 0;
    gcd(p,q,x,k);
    return (p%k+k)%k;
}
int main(){
    tt.a[1][1]=1;
    tt.a[2][1]=1;
    tt.a[1][2]=1;
    sq u;
    u.a[1][1]=1;
    a[1]=1;a[2]=1;
    scanf("%lld%lld%lld",&n,&k,&p);
    for(ll i=3;i<=k*3;i++){
        a[i]=(a[i-1]+a[i-2])%k;
        if(!fr[a[i]])fr[a[i]]=i;
    }
    ll i=0;
    ll j=1;
    cir=n+1;
    ll ss=0;
    for(;;){
        ll u=ni(j);
        if(!u||!fr[u])break;
        ll g=a[fr[u]-1]*j%k;
        if(i+fr[u]>n)break;
        if(!ben[g]){
            ben[g]=1;
            bb[++bn]=i+fr[u];
            hash[g]=bn;
        }else{
            bb[++bn]=i+fr[u];
            ss=hash[g];
            cir=bb[bn]-bb[ss];
            break;
        }
        i=i+fr[u];
        j=g;
    }
    sq yy=ff(tt,cir);
    res=(u*ff(tt,n)).a[1][2];
    for(ll i=1;i<=bn;i++){
        if(i<=ss){
            (res-=(u*ff(tt,n-bb[i])).a[1][1])%=p;
        }else{ll uu=(n-bb[i])%cir;
        (res-=((u*(ff(tt,uu)*cheng(yy,(n-bb[i])/cir)))).a[1][1])%=p;
    }}
    res=(res%p+p)%p;
    cout<<res;
}