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

版权声明:欢迎转载(请附带原链接)ヾ(๑╹◡╹)ノ" https://blog.csdn.net/corsica6/article/details/84134841

传送门:bzoj2432


题解

良心题,暴力有 70 p t s 70pts 。。。

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

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

观察到一些特殊的性质:

  • 每行开头必然是两个相同的数(设为 s t i st_i ),且每一行形式相同: s t i f 1 , s t i f 2 , s t i f 3 . . . st_i*f_1,st_i*f_2,st_i*f_3...

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

  • 循环节中每行结尾都是以 s t i f l e n 1 m o d    k st_i*f_{len}\equiv 1\mod k 结束的:
    3 , 3 , 6 , 2 , 0 3,3,6,2,0
    2 , 2 , 4 , 6 , 3 , 2 , 5 , 0 , 5 , 5 , 3 , 0 2,2,4,6,3,2,5,0,5,5,3,0
    这样表示比较直观。

只需要找到最小的 f i s t i 1 m o d    k f_i\equiv st_i^{-1}\mod k 就得到了第 i i 行的长度,且 s t i + 1 = s t i f l e n 1 m o d    k st_{i+1}=st_i*f_{len-1}\mod k

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

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

扫描二维码关注公众号,回复: 4246083 查看本文章

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

行内的转移矩阵:
[ 1 1 0 1 0 0 0 0 1 ] \left[ \begin{matrix} 1& 1& 0\\ 1&0&0\\ 0&0&1 \end{matrix} \right]
行之间的转移矩阵:
[ 1 0 0 0 1 0 1 0 1 ] \left[ \begin{matrix} 1& 0& 0\\ 0&1&0\\ -1&0&1 \end{matrix} \right]
初始矩阵:
[ f 0 = 0 f 1 = 1 1 ] \left[ \begin{matrix} f_0=0&f_{-1}=1 &1 \end{matrix} \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;
}

猜你喜欢

转载自blog.csdn.net/corsica6/article/details/84134841