【BZOJ2432】【NOI2011】兔农（数论，矩阵快速幂）
【BZOJ2432】【NOI2011】兔农（数论，矩阵快速幂）
题面
题解
这题\(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;
}
