【BZOJ2432】【NOI2011】兔农（数论，矩阵快速幂）

题面

BZOJ

题解

这题$75$分就是送的，我什么都不想写。

先手玩一下，发现每次每次出现$mod\ K=1$的数之后
把它减一，就变成了$0$。接着后面的数显然还是一个斐波那契数列
只是都乘了$0$之前的那个数作为倍数而已。
拿样例举个例子？以下数字都在模$7$意义下进行

1 1 2 3 5 0(1)
5 5 3 0(1)
3 3 6 2 0(1)
大概就是这样子。

当然，如果我们继续手玩下去，也许可以发现点什么？

1 1 2 3 5 0(1)
5 5 3 0(1)
3 3 6 2 0(1)
2 2 4 6 3 2 5 0
5 5 3 0(1)

似乎出现了循环？？？
那么，我们似乎可以按照找到末尾的$0$，找下一行的零，找到循环节这样的步骤来。

至于这个循环节的长度相关的问题，可以看看Vfk的博客。orz

考虑一下怎么计算这个$0$的位置？事实上是在找$1$的位置
而上面举例中的每一行都是一个斐波那契数列乘上$x$
其中$x$是上一行中倒数第二个数字
那么，$fib[len]*x=1$，而$x$对于我们来说是一个已知项
所以这个过程变成了一个求逆的过程。
而有根据$vfk$的博客，斐波那契数列在模$K$意义下的循环节长度不超过$6K$
所以我们可以暴力算一个循环节的斐波那契数列
这样子，我们的过程就变成了
找到当前行末尾的位置，对应的乘一下，得到下一行的倍数
如果下一行的开头这个数字已经被得到过，那么出现了全局的循环节，
直接暴力算就可以了。
如果发现此时逆元不存在，证明没有循环，直接矩阵快速幂即可。
否则的话继续矩阵快速幂找下一行即可。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define RG register
#define MAX 1001000
inline ll read()
{
    RG ll x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
ll MOD;
void add(ll &x,ll y){x+=y;if(x>=MOD)x-=MOD;}
struct Matrix
{
    ll s[4][4];
    ll* operator[](int x){return s[x];}
    void clear(){memset(s,0,sizeof(s));}
    void init(){clear();s[1][1]=s[2][2]=s[3][3]=1;}
}nt,lt,ans,ret[MAX];
Matrix operator*(Matrix a,Matrix b)
{
    Matrix ret;ret.clear();
    for(int i=1;i<=3;++i)
        for(int j=1;j<=3;++j)
            for(int k=1;k<=3;++k)
                add(ret[i][j],1ll*a[i][k]*b[k][j]%MOD);
    return ret;
}
Matrix fpow(Matrix a,ll b)
{
    Matrix s;s.init();
    while(b){if(b&1)s=s*a;a=a*a;b>>=1;}
    return s;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1;y=0;return;}
    exgcd(b,a%b,x,y);
    ll tmp=y;
    y=x-(a/b)*y;x=tmp;
}
ll f[MAX<<3],n,K,vis[MAX],inv[MAX],len[MAX];
bool book[MAX];
int main()
{
    n=read();K=read();MOD=read();
    f[1]=f[2]=1;bool fl=false;
    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]==1&&f[i-1]==1)break;
    }
    nt[1][2]=nt[2][1]=nt[2][2]=nt[3][3]=1;
    lt.init();lt[3][2]=-1;
    ans[1][1]=ans[1][3]=1;
    for(ll t=1;n;)
    {
        if(!inv[t])
        {
            if(__gcd(t,K)!=1)inv[t]=-1;
            else
            {
                ll x,y;exgcd(t,K,x,y);
                inv[t]=(x+K)%K;
            }
        }
        if(inv[t]==-1){ans=ans*fpow(nt,n);break;}
        if(!book[t]||fl)
        {
            book[t]=true;
            if(!vis[inv[t]]){ans=ans*fpow(nt,n);break;}
            len[t]=vis[inv[t]];
            if(n>=len[t])
            {
                n-=len[t];
                ret[t]=fpow(nt,len[t])*lt;
                ans=ans*ret[t];
                t=t*f[len[t]-1]%K;
            }
            else{ans=ans*fpow(nt,n);break;}
        }
        else
        {
            Matrix now;ll cnt=0;now.init();
            for(ll i=t*f[len[t]-1]%K;i!=t;i=i*f[len[i]-1]%K)
                now=now*ret[i],cnt+=len[i];
            now=ret[t]*now;cnt+=len[t];
            ans=ans*fpow(now,n/cnt);
            n%=cnt;fl=true;
        }
    }
    printf("%lld\n",(ans[1][2]%MOD+MOD)%MOD);
    return 0;
}