[NOI2011]兔农 - jcvb

[NOI2011]兔农

jcvb posted @ 2013年6月07日 21:09 with tags 数论 矩阵乘法 , 5693 阅读

题意:f[1]=f[2]=1;i>2时,f[i]=f[i-1]+f[i-2],  if(f[i] % k == 1) f[n]--

求f[n]%p。

白送75分:同时记录f[i]%k和f[i]%p的值,暴力算即可。

100分算法:

以k=7为例,考虑f[i]%k组成的序列:

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,

把减1得0的位置标出,并以这些0为界分段,可以发现:

①每段开头必为相同两数,它恰是上一段的最末一位非0数;由于总共只有k-1种余数,所以不超过k段就会出现循环(如果有的话),比如上面k=7时的第3,4段就是循环节。

②记斐波那契数列为fib[i]。假如某段段首数字为x,那么这一段内第i个数即为x*fib[i] % k。若记这一段长度为len,则有x*fib[len]==1 (mod k)。


现在我们试图找到整个数列的循环结构:根据上式,①求x的逆元得到fib[len],②由fib[len]得知len,③用x*fib[len-1]%k算出下一段的段首,重复操作直到发现循环(或者发现这一段永远不终止)。
至于具体实现:①扩欧或者欧拉定理②预处理indfib[y]数组,表示斐波那契数列中模k余y的数第一次出现的下标③预处理fib[i]模k的值。有一个结论:斐波那契数列模k后一定是0,1,1开头的纯循环,而且这个循环节的长度<=6k(不知具体怎么证。。),所以只需暴力算fib数组并同时记录indfib[],发现循环即停止。
注意,假如第①步不存在逆元,或者第②步不存在符合的len,那么这一段将永远不会终止(比如k=8时就是这样),那么整个数列就不存在循环了(可以视作最后一段的长度为无穷大)。

 

接下来考虑如何用矩阵乘法计算f[n]%p。

[tex]\begin{pmatrix}
a_{i-1} & a_i & 1
\end{pmatrix}\begin{pmatrix}
0 &1  &0 \\
 1&1  &0 \\
 0&0  &1
\end{pmatrix}=\begin{pmatrix}
a_i &a_{i-1}+a_i  &1
\end{pmatrix}[/tex]

[tex]\begin{pmatrix}
a_{i-1} & a_i & 1
\end{pmatrix}\begin{pmatrix}
1 &0  &0 \\
 0&1  &0 \\
 0&-1  &1
\end{pmatrix}=\begin{pmatrix}
a_{i-1} &a_i-1  &1
\end{pmatrix}[/tex]

分别记这两个3*3矩阵为A,B。令初始矩阵为[tex]\begin{pmatrix}
1 & 0& 1
\end{pmatrix}[/tex],通过对其不断右乘A和B便能实现累加、减1两种操作。对于分出的每一段算出一个矩阵A^len * B,表示这一段的“效果”。

接下来是喜闻乐见的分类讨论时间:假如整个数列是循环的,判断第n项是否在混循环的那部分里,若是则直接把前面几段乘起来,n所在这一段的零头直接用A的次幂算;若不是则先把混循环全部乘起来,然后把循环节全部乘起来,算出循环次数再快速幂,然后再像刚才一样算零头乘上去。若数列不循环倒方便些,也与上面类似,不多说了。

代码很丑陋

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;
ll n,k,p;
ll fib[10000005],indfib[1000005]={0},inv[1000005];
ll exgcd(ll a,ll b,ll &lx,ll &ly){
    ll x=0,y=1,q,t;lx=1,ly=0;
    while(b){
    	q=a/b;
        t=b;b=a%b;a=t;
        t=x;x=lx-q*x;lx=t;
        t=y;y=ly-q*y;ly=t;
    }
    return a;
}
void buildinv(){
	ll tmp1,tmp2;inv[1]=1;
	for (ll i=2;i<k;i++)
		if(exgcd(i,k,tmp1,tmp2)!=1)inv[i]=0;
		else inv[i]=tmp1<0?tmp1+k:tmp1;
}
ll ord[1000005];ll tot=0;
ll lensum[1000005];
ll vis[1000005]={0};
struct mat{ll a[3][3];void cl(){memset(a,0,sizeof(a));}void set1(){cl();a[0][0]=a[1][1]=a[2][2]=1;}};
mat mul(const mat &a,const mat &b){
	mat c;
	for (int i=0;i<3;i++)
		for (int j=0;j<3;j++){
			c.a[i][j]=0;
			for (int k=0;k<3;k++)c.a[i][j]+=a.a[i][k]*b.a[k][j];
			c.a[i][j]%=p;
		}
	return c;
}
mat po(mat a,ll t){
	mat c;c.set1();
	while(t){
		if(t&1)c=mul(c,a);
		a=mul(a,a);
		t>>=1;
	}
	return c;
}
mat s[1000005],A,B;
int main()
{	
	cin>>n>>k>>p;
	fib[0]=0,fib[1]=fib[2]=1;
	for (ll i=3;;i++){ 
		fib[i]=fib[i-1]+fib[i-2];if(fib[i]>=k)fib[i]-=k;
		if(!indfib[fib[i]])indfib[fib[i]]=i;
		if(fib[i]==1 && fib[i-1]==1)break;
	}
	buildinv();
	ll cur=1;lensum[0]=0;
	A.cl();A.a[0][1]=A.a[1][0]=A.a[1][1]=A.a[2][2]=1;
	B.cl();B.a[0][0]=B.a[1][1]=B.a[2][2]=1;B.a[2][1]=-1;
	while(1){
		if(vis[cur])break;
		if(!inv[cur] || !indfib[inv[cur]])break;
		tot++;vis[ord[tot]=cur]=tot;
		lensum[tot]=lensum[tot-1]+indfib[inv[cur]];
		s[tot]=po(A,indfib[inv[cur]]);s[tot]=mul(s[tot],B);
		cur=1LL*cur*fib[indfib[inv[cur]]-1]%k;
	}	
	mat N;N.set1();
	if(vis[cur]){
		int pre=vis[cur]-1,i;mat M;
		for (i=1;i<=pre && n>=lensum[i];i++)N=mul(N,s[i]);
		if(i<=pre){
			M=po(A,n-lensum[i-1]);
			N=mul(N,M);
		}else{
			ll cylen=lensum[tot]-lensum[pre];
			n-=lensum[pre];
			mat C;C.set1();
			for (i=pre+1;i<=tot;i++)C=mul(C,s[i]);
			C=po(C,n/cylen);N=mul(N,C);n%=cylen;
			for (i=pre+1;n>=lensum[i]-lensum[pre];i++)N=mul(N,s[i]);
			M=po(A,n-(lensum[i-1]-lensum[pre]));
			N=mul(N,M);
		}
	}else{
		mat M;int i;
		for (i=1;i<=tot && n>=lensum[i];i++)N=mul(N,s[i]);
		M=po(A,n-lensum[i-1]);
		N=mul(N,M);
	}
	ll ans=(N.a[0][1]+N.a[2][1])%p;if(ans<0)ans+=p;
	cout<<ans<<'\n';
	return 0;
}
Host by is-Programmer.com | Power by Chito 1.3.3 beta | © 2007 LinuxGem | Design by Matthew "Agent Spork" McGee