题解：P2020 [NOI2011] 兔农

link

对于一般的斐波那契数列，显然可以用矩阵快速幂求解。在这道题中，当 \(F_{i-2}+F_{i-1}\equiv 1\pmod k\) 的时候，\(F_i\) 就需要减 \(1\)。如果我们设答案矩阵如下：

\[A=\begin{bmatrix}F_{i-1}&F_i&1\end{bmatrix} \]

则可以列出普通的转移矩阵：

\[B=\begin{bmatrix}0&1&0\\1&1&0\\0&0&1\end{bmatrix} \]

同理可以列出减 \(1\) 时的转移矩阵：

\[C=\begin{bmatrix}0&1&0\\1&1&0\\0&-1&1\end{bmatrix} \]

接下来我们考虑该以什么顺序来乘这些转移矩阵，对 \(F_i\mod k\) 打表（\(k=7\)，每次减 \(1\) 就换行），有：

 1 1 2 3 5 0
 5 5 3 0
 3 3 6 2 0
 2 2 4 6 3 2 5 0 5 5 3 0
 3 3 6 2 0
 2 2 4 6 3 2 5 0 5 5 3 0
 3 3 6 2 0
 2 ...

考虑每一行，设第一项为 \(x\)，由于上一行的末尾一定为 \(0\)，则这一行的第 \(j\) 项的值即为 \(fib_jx\)。若该行的长度为 \(l\)，则有 \(fib_lx\equiv1\pmod k\)。由于 \(k\) 不大，我们可以通过枚举 \(fib_i\)（当 \(fib_i\equiv fib_{i-1}\equiv 1\pmod k\) 时就结束），利用逆元求出对于每个 \(x\) 开头的行的长度和倒数第二个数。

发现出现了循环节，由于在模 \(k\) 意义下，显然循环节的行数是 \(O(k)\) 的，于是我们可以求出循环节。考虑完整的行的矩阵为 \(B^{l-1}C\)，分讨 \(n\) 是否涉及循环节、包含几个循环节，然后将循环节中的散行依次乘起来即可。

但还有一种可能是当前的 \(x\) 没有逆元，那么直接乘上 \(B^n\) 就行了（当然 \(n\) 是此时剩下的 \(n\)）。

注意矩阵乘法是在模 \(p\) 意义下，其他地方是在模 \(k\) 意义下，时间复杂度为 \(O(9k\log n)\)

参考资料：题解 P2020 【[NOI2011]兔农】

Code:

//kon ha watashi no suki desu
#include<iostream>
#include<cstring>
#define rep(i,l,r) for(int i=(l);i<=(r);i++)
#define per(i,l,r) for(int i=(l);i>=(r);i--)
using namespace std;
typedef long long ll;
const int maxk=1e6+5;
const ll inf=0x3f3f3f3f3f3f3f3f;
ll af[maxk<<3],al[maxk],siz[maxk],in,ik,ip;
int vis[maxk];
struct matrix{
	ll a[4][4];
	inline ll* operator[](int x){return a[x];}
	inline matrix operator*(matrix b){
		matrix res;
		rep(v1,1,3)rep(v2,1,3)rep(v3,1,3)res[v1][v2]=(res[v1][v2]+a[v1][v3]*b.a[v3][v2])%ip;
		return res;
	}
	inline void build(){
		rep(v1,1,3)a[v1][v1]=1;
	}
	inline void print(){
		rep(v1,1,3){
			rep(v2,1,3)cout<<a[v1][v2]<<' ';
			cout<<endl;
		}
	}
	matrix(){
		memset(a,0,sizeof(a));
	}
};
inline matrix qpow(matrix a,ll b){
	matrix res;
	res.build();
	while(b){
		if(b&1)res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
inline void ex_gcd(ll a,ll b,ll &x,ll &y,ll &g){
	if(!b){
		x=1;
		y=0;
		g=a;
		return;
	}
	ex_gcd(b,a%b,x,y,g);
	ll t=x;
	x=y;
	y=t-a/b*y;
}
int main()
{
	memset(al,0x3f,sizeof(al));
	cin>>in>>ik>>ip;
	if(in<=2){
		cout<<1<<endl;
		return 0;
	}
	af[1]=af[2]=1;
	rep(v1,3,1e7){
		af[v1]=(af[v1-1]+af[v1-2])%ik;
		if(af[v1]==1)al[af[v1]]=min(al[af[v1]],(ll)v1);
		if(af[v1]==1&&af[v1-1]==1)break;
		ll x,y,g;
		ex_gcd(af[v1],ik,x,y,g);
		x=(x%ik+ik)%ik;
		if(g==1)al[x]=min(al[x],(ll)v1);
	}
	ll p=1,tot=0,pre=0;
	bool flag=false;
	while(1){
		if(al[p]==inf){
			rep(v1,1,tot)pre+=siz[v1];
			flag=true;
			vis[p]=tot+1;
			break;	
		}
		if(vis[p]){
			rep(v1,1,vis[p]-1)pre+=siz[v1];
			break;
		}
		siz[++tot]=al[p];
		vis[p]=tot;
		p=p*af[al[p]-1]%ik;
	}
	matrix ans,tr1,tr2;
	ans[1][2]=ans[1][3]=tr1[1][2]=tr1[2][1]=tr1[2][2]=tr1[3][3]=1;
	tr2=tr1;
	tr2[3][2]=ip-1;
	in--;
	pre--;
	siz[1]--;
	if(in<=pre){//不涉及循环节
		rep(v1,1,tot){
			if(siz[v1]<=in){
				ans=ans*qpow(tr1,siz[v1]-1)*tr2;
				in-=siz[v1];
			}
			else{
				ans=ans*qpow(tr1,in);
				break;
			}
		}
	}
	else{//涉及循环节
		rep(v1,1,vis[p]-1){//处理循环节之前的行
			ans=ans*qpow(tr1,siz[v1]-1)*tr2;
			in-=siz[v1];
		}
		if(flag){//最后一行无逆元
			ans=ans*qpow(tr1,in);
		}
		else{
			ll lp=0;
			matrix fac;
			fac.build();
			rep(v1,vis[p],tot){
				fac=fac*qpow(tr1,siz[v1]-1)*tr2;
				lp+=siz[v1];
			}
			ans=ans*qpow(fac,in/lp);//整体的循环节
			in%=lp;
			if(in){
				rep(v1,vis[p],tot){//循环节内的散行
					if(siz[v1]<=in){
						ans=ans*qpow(tr1,siz[v1]-1)*tr2;
						in-=siz[v1];
					}
					else{
						ans=ans*qpow(tr1,in);
						break;
					}
				}
			}
		}
	}
	cout<<ans[1][2]<<endl;
	return 0;
}