bzoj2432: [Noi2011]兔农 快速幂+数论

不难发现，这个题就是求斐波那契数列改化，由于有一个很强的结论，斐波那契数列取模是一个周期数列，所以我们可以去找循环节，然后找到循环节后把这第一个循环节处理出来。 其实vfk说的很详细了，注意这里mod的数不一定是个质数，我们只能用拓展欧几里得求逆元。。。。

http://vfleaking.blog.163.com/blog/static/174807634201341721051604/

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <iostream>
using namespace std;
#define oo 2000000000000000000000LL
#define ll long long
#define maxn 1000100
ll n,len[maxn],next[maxn];
int inv[maxn],k,mod;
void init()
{
    inv[1]=1;
    for(int i=2;i<k;i++)
    {
        inv[i]=(k-k/i)*inv[i]%k;
    }
}
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 ni(ll x)
{
    ll p,q;
    if(gcd(x,k)!=1)return 0;
    gcd(p,q,x,k);
    return (p%k+k)%k;
}
int pos[maxn],fib[maxn*6];
void firstap()
{
    memset(pos,-1,sizeof(pos));
    fib[0]=1;fib[1]=1;fib[2]=2;
    for(int i=3;!( fib[i-1]==1&& fib[i-2]==1);i++)
    {
        fib[i]=fib[i-1]+fib[i-2];
        if(fib[i]>=k) fib[i]-=k;
        if(pos[fib[i]]==-1) pos[fib[i]]=i;
    }
}
struct matrix
{
    int a[3][3];
    void cl()
    {
        memset(a,0,sizeof(a));
    }
    void one()
    {
        memset(a,0,sizeof(a));
        for(int i=0;i<3;i++) a[i][i]=1;
    }
    void A()
    {
        a[0][0]=1;a[0][1]=1;a[0][2]=0;
        a[1][0]=1;a[1][1]=0;a[1][2]=0;
        a[2][0]=0;a[2][1]=0;a[2][2]=1;
    }
    void B()
    {
        a[0][0]=1;a[0][1]=1;a[0][2]=0;
        a[1][0]=1;a[1][1]=0;a[1][2]=0;
        a[2][0]=mod-1;a[2][1]=0;a[2][2]=1;
    }
};
matrix operator *(matrix aa,matrix bb)
{
    ll a[3][3];
    for(int i=0;i<3;i++)
    for(int j=0;j<3;j++)
    {
        a[i][j]=0;
        for(int k=0;k<3;k++)
        {
            a[i][j]+=(ll)aa.a[i][k]*bb.a[k][j]%mod;
            a[i][j]%=mod;
        }
    }
 
    matrix cc;
    for(int i=0;i<3;i++)
        for(int j=0;j<3;j++)
        {
            cc.a[i][j]=a[i][j];
        }
    return cc;
}
matrix matpow(matrix x,ll y)
{
    matrix res;res.one();
    while(y)
    {
        if(y&1){
            res=res*x;
        }
        x=x*x;
        y>>=1;
    }
    return res;
}
int vis[maxn];
int main()
{
    //freopen("rabbit20.in","r",stdin);
    scanf("%lld%d%d",&n,&k,&mod);
    for(int i=1;i<k;i++) inv[i]=ni(i);
 
    firstap();
    for(int i=0;i<k;i++)
    {
        if(!inv[i]||!pos[inv[i]])
        {
            next[i]=-1;len[i]=oo;
        }
        else
        {
            int kk=pos[inv[i]];
            next[i]=(ll)fib[kk-1]*i%k;
            len[i]=kk+1;
        }
    }
    int cirs;
    for (cirs=1;cirs!=-1 && !vis[cirs];cirs=next[cirs])
        vis[cirs]=1;
    matrix use;
    use.one();
    matrix AA,BB; AA.A();BB.B();
    for(int i=1;i!=cirs;i=next[i])
    {
        if(n<len[i])
        {
            use=use*matpow(AA,n);
            n=0;
            break;
        }
        else
        {
            use=use*matpow(AA,len[i]-1);
            use=use*BB;
            n-=len[i];
            if(n==0)
            {
                break;
            }
        }
    }
    if(n>0)
    {
        int cur = cirs;
        ll totL = 0;
        matrix matC;matC.one();
        do
        {
            matC = matC*matpow(AA,len[cur]-1);
            matC = matC*BB;
            totL += len[cur];
            cur = next[cur];
        }
        while (cur != cirs);
        use=use*matpow(matC, n / totL);
        n=n%totL;
        for(int i=cirs;;i=next[i])
        {
            if(n<len[i])
            {
                use=use*matpow(AA,n);
                break;
            }
            else
            {
                use=use*matpow(AA,len[i]-1);
                use=use*BB;
                n-=len[i];
                if(n==0)
                {
                    break;
                }
            }
        }
    }
    ll res=0;
    res=(ll)use.a[1][0]+(ll)use.a[2][0];
    res= res%mod;
    cout << res << endl;
    return 0;
}