[NOI2011]兔农(斐波那契数列+乘法逆+矩阵加速)

12 篇文章 0 订阅
6 篇文章 0 订阅

几乎是抄的,啥都不要说了

贴上策爷和VFK的链接好了

http://jcvb.is-programmer.com/posts/39528.html

http://vfleaking.blog.163.com/blog/static/174807634201341721051604/


还是写两句吧。关键是将新数列每个元素模k后写出,在经减1余0的元素后面换行,并找规律

发现每行前两个数相同,因此每行第i个数为 Fibonacc[i]*行首的数,行末数为0,行长len[i]为:第一个模k等于[行首的数在模k下的逆元]的Fibonacc数的位置,每一行的后继行行首数字next[i]行首的数*Fibonacc[len[i]-1]%k

把这些信息预处理出来会方便一些

还发现有几行会循环出现,考虑矩阵加速优化之:

从第1行开始,每行内部用矩阵加速处理 
若n很大,进入循环,先处理一个循环节,然后在循环节之间也用矩阵加速处理 
最后,不足一循环的,每行矩阵加速处理 

注意特殊情况:不出现循环,即某一行的行首在模k下不存在逆元,或其逆元不等于任何Fibonacc[x]%k,此时设这一行无限长,无后继行 


#include<stdio.h>
#include<stdlib.h>
#include<string.h>
typedef long long LL;
const LL INF=2000000000000000000ll;//结尾加"ll"表示long long型常量 
LL len[1000005]={0};
int fib[6000005]={0},first[1000005]={0},next[1000005]={0},vis[1000005]={0};
LL mod;
struct juzhen
{
	LL s[5][5];
	juzhen()
	{
		memset(s,0,sizeof(s));
	}
};
void exgcd(int a,int b,int& d,int& x,int& y)
{
	int t;
	if(b==0)
	{
		d=a;
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,d,x,y);
	t=x;
	x=y;
	y=t-(a/b)*y;
}
int ni(int a,int n)
{
	int d=0,x=0,y=0;
	exgcd(a,n,d,x,y);
	if(d!=1) return 0;
	return (x+n)%n;
}
juzhen cheng(juzhen a,juzhen b)
{
	juzhen res;
	int i,j,k;
	for(i=1;i<=3;i++)
		for(j=1;j<=3;j++)
		{
			for(k=1;k<=3;k++)
				res.s[i][j]+=a.s[i][k]*b.s[k][j];
			res.s[i][j]%=mod;
		}
	return res;
}
juzhen ksm(juzhen a,LL n)
{
	juzhen res;
	res.s[1][1]=res.s[2][2]=res.s[3][3]=1;
	if(n==0) return res;
	if(n==1) return a;
	res=ksm(a,n/2);
	res=cheng(res,res);
	if(n%2==1) res=cheng(res,a);
	return res;
}
int main()
{
	juzhen A,B,C,Z;
	LL n,len_cir=0;
	int k,i,cir;
	scanf("%lld%d%lld",&n,&k,&mod);
	fib[1]=fib[2]=1;
	for(i=3;i==3||fib[i-1]!=1||fib[i-2]!=1;i++)//预处理出Fibonacc数列每个元素模k的值,直到出现循环 
	{
		fib[i]=(fib[i-1]+fib[i-2])%k;
		if(first[fib[i]]==0) first[fib[i]]=i;//记录余数最早出现的位置 
	}
	for(i=1;i<k;i++)
	{
		len[i]=(LL)first[ni(i,k)];//以两个i开头,到i的逆元的数字个数 
		if(ni(i,k)!=0&&len[i]!=0) next[i]=(int)((LL)fib[len[i]-1]*(LL)i%(LL)k);
		else//i无逆元或Fibonacc数列中不存在i的逆元
		{
			next[i]=-1;//这一行无限休无止,因此不存在下一行 
			len[i]=INF;//这一行无限长 
		}
	}
	for(cir=1;vis[cir]==0&&cir!=-1;cir=next[cir])//找圈的起点cir
		vis[cir]=1;
	A.s[1][2]=A.s[2][1]=A.s[2][2]=A.s[3][3]=1;//求Fibonacc数列的下一项 
	B.s[1][1]=B.s[2][2]=B.s[3][3]=1;//当前项减1
	B.s[3][2]=-1;
	Z.s[1][1]=Z.s[1][3]=1;
	for(i=1;i!=cir;i=next[i])//在圈之前的部分 
	{
		if(n<len[i])
		{
			Z=cheng(Z,ksm(A,n));
			n=0;
			break;
		}
		Z=cheng(Z,ksm(A,len[i]));
		Z=cheng(Z,B);
		n-=len[i];
	}
	C.s[1][1]=C.s[2][2]=C.s[3][3]=1;
	if(n>0)
	{
		do
		{
			C=cheng(C,ksm(A,len[i]));
			C=cheng(C,B);
			len_cir+=len[i];//计算圈上总共有多少项 
			i=next[i];
		}
		while(i!=cir);
		Z=cheng(Z,ksm(C,n/len_cir));
		n%=len_cir;
		for(i=cir;;i=next[i])//剩下不完整圈的部分 
		{
			if(n<len[i])
			{
				Z=cheng(Z,ksm(A,n));
				n=0;
				break;
			}
			Z=cheng(Z,ksm(A,len[i]));
			Z=cheng(Z,B);
			n-=len[i];
		}
	}
	printf("%lld",Z.s[1][2]);
	return 0;
}


  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值