【BZOJ】2432: [Noi2011]兔农 -矩乘&找规律&fibnacci循环节

传送门：bzoj2432

---

题解

良心题，暴力有$70pts$。。。

显然是分别维护斐波那契数列$\mathrm{m}\mathrm{o}\mathrm{d}k,\mathrm{m}\mathrm{o}\mathrm{d}p$意义下的值。
就题面中的例子($\mathrm{m}\mathrm{o}\mathrm{d}7$)来看(方便观察性质，每出现一个$0$就新开一行)：

$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$
$...$

观察到一些特殊的性质：

• 每行开头必然是两个相同的数(设为$s{t}_i$)，且每一行形式相同：$s{t}_i\ast f_1,s{t}_i\ast f_2,s{t}_i\ast f_3...$

• 如果有循环节的话，循环行数必然不超过$k$。

• 循环节中每行结尾都是以$s{t}_i\ast f_{len}\equiv 1\mathrm{m}\mathrm{o}\mathrm{d}k$结束的：
  $3,3,6,2,0$
  $2,2,4,6,3,2,5,0,5,5,3,0$
  这样表示比较直观。

只需要找到最小的$f_i\equiv s{t}_i^{-1}\mathrm{m}\mathrm{o}\mathrm{d}k$就得到了第$i$行的长度，且$s{t}_{i+1}=s{t}_i\ast f_{len-1}\mathrm{m}\mathrm{o}\mathrm{d}k$

具体来说，用$v{s}_i$表示最小的模$k$等于$i$的斐波那契数项。而可以证明模$k$意义下，斐波那契循环节长度不超过$6\ast k$，且循环节中前三个数必然为$0,1,1$(一个paper)

$exgcd$求解逆元，但$k$不保证是质数，无解的情况就等同于没有特殊的减一操作，直接矩乘。

所以我们只需要在不超过$k$次内找出循环节然后矩阵快速幂即可。

行内的转移矩阵：
$\left[\begin{array}{ccc}1& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]$
行之间的转移矩阵：
$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ -1& 0& 1\end{array}\right]$
初始矩阵:
$\left[\begin{array}{ccc}{f}_0=0& f_{-1}=1& 1\end{array}\right]$

---

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+10;

ll n;bool fwg[N];
int K,mod,pr,vs[N],f[N*6],inv[N];

inline int ad(int x,int y){x+=y;return x>=mod?x-mod:x;}
inline int dc(int x,int y){x-=y;return x<0?x+mod:x;}

struct Mat{
	int g[3][3];
	inline int *operator[](int x){return g[x];}
	Mat operator *(const Mat&ky)const{
	    Mat re;register int i,j,k,v;
	    for(i=0;i<3;++i)
	     for(j=0;j<3;++j){
	     	for(v=k=0;k<3;++k) v=ad(v,(ll)g[i][k]*ky.g[k][j]%mod);
	     	re.g[i][j]=v;
	     }
	    return re;
	}
	inline void itia()
	{
		register int i,j;
		for(i=0;i<3;++i)
		 for(j=0;j<3;++j)
		  g[i][j]=(i==j);
	}
}A,B,C,D,ori,mt[N];

inline Mat fp(Mat a,ll y)
{
	ori.itia();
	for(;y;y>>=1,a=a*a)
	 if(y&1) ori=ori*a;
	return ori;
}


int exgcd(int a,int b,ll &x,ll &y)
{
	if(!b) {x=1;y=0;return a;}
	int re=exgcd(b,a%b,y,x);y-=x*(a/b);
	return re;
}

inline int gtiv(int a,int b)
{
	ll x,y;
	return exgcd(a,b,x,y)==1?((x%b+b)%b):(-1);
}

int main(){
    int i,x,z,len;ll cot;
    scanf("%lld%d%d",&n,&K,&mod);
    f[1]=f[2]=1;
	for(i=3;;++i){
		f[i]=(f[i-1]+f[i-2])%K;
		if(!vs[f[i]]) vs[f[i]]=i;
		if(f[i]==1 && f[i]==f[i-1]) break;
    }
    A[0][0]=A[0][1]=A[1][0]=A[2][2]=1;
    B[0][0]=B[1][1]=B[2][2]=1;B[2][0]=mod-1;
    C[0][1]=C[0][2]=1;
    for(z=1,pr=0;n;){
    	if(!inv[z]) inv[z]=gtiv(z,K);
    	if(inv[z]==-1) {C=C*fp(A,n);break;}
    	if(pr || (!fwg[z])){
    		len=vs[inv[z]];
			if(!len || n<len){C=C*fp(A,n);break;}
    		if(!fwg[z]){fwg[z]=true;mt[z]=fp(A,len)*B;}
    		n-=len;C=C*mt[z];z=(ll)z*f[len-1]%K;
    	}else{
    		cot=0;cot=vs[inv[z]];D=mt[z];
    		for(x=(ll)z*f[vs[inv[z]]-1]%K;x!=z;x=(ll)x*f[vs[inv[x]]-1]%K)
    		  D=D*mt[x],cot+=vs[inv[x]];
    		C=C*fp(D,n/cot);n%=cot;pr=1;
    	}
    }
    printf("%d",C[0][0]);
    return 0;
}