矩阵乘+坑 [NOI2011] 兔农

     

问题 C: 兔农

时间限制: 1 Sec  内存限制: 256 MB

题目描述

农夫栋栋近年收入不景气，正在他发愁如何能多赚点钱时，他听到隔壁的小 朋友在讨论兔子繁殖的问题。 问题是这样的：第一个月初有一对刚出生的小兔子，经过两个月长大后，这 对兔子从第三个月开始，每个月初生一对小兔子。新出生的小兔子生长两个月后 又能每个月生出一对小兔子。问第 n 个月有多少只兔子？ 聪明的你可能已经发现，第 n 个月的兔子数正好是第 n 个 Fibonacci(斐波那 契)数。栋栋不懂什么是 Fibonacci 数，但他也发现了规律：第 i+2 个月的兔子数 等于第 i 个月的兔子数加上第 i+1 个月的兔子数。前几个月的兔子数依次为： 1 1 2 3 5 8 13 21 34 … 栋栋发现越到后面兔子数增长的越快，期待养兔子一定能赚大钱，于是栋栋 在第一个月初买了一对小兔子开始饲养。 每天，栋栋都要给兔子们喂食，兔子们吃食时非常特别，总是每 k 对兔子围 成一圈，最后剩下的不足 k 对的围成一圈，由于兔子特别害怕孤独，从第三个月 开始，如果吃食时围成某一个圈的只有一对兔子，这对兔子就会很快死掉。 我们假设死去的总是刚出生的兔子，那么每个月的兔子数仍然是可以计算的。 例如，当 k=7 时，前几个月的兔子数依次为： 1 1 2 3 5 7 12 19 31 49 80 … 给定 n，你能帮助栋栋计算第 n 个月他有多少对兔子么？由于答案可能非常 大，你只需要告诉栋栋第 n 个月的兔子对数除 p 的余数即可。

输入

输入一行，包含三个正整数 n, k, p。

输出

输出一行，包含一个整数，表示栋栋第 n 个月的兔子对数除 p 的余数。

样例输入

6 7 100

样例输出

7

提示

    好一道国赛的题，觉得巨坑很好。。

    菲波数%k所得到的数排成排，在需要-1（也就是%k==1）时，重启一行，就会发现能搞出一个循环，也不一定，我下面说。。

    有两个结论，我不会证，但对这道题很有帮助，

     一：如果有循环节，最多找6*k个数，就一定能循环一次，

    二：设x为这一行行首的值，x*fib[i]%k==这一行第i位的值。

    然后发现，x*fib[i]%k==1就是行尾，而fib[i]即为x的乘法逆元，所以能求出fib[i]的值，

    并且如结论一所提，在这里找时，记录某个fib[i]第一次出现的位置，因此求出逆元后，就变相得知了这行的长度。当把循环节求出后，就好说了。。

    注：有可能出现每个fib[i]%k都无逆元的情况，就是说这一行无限延伸下去了。。直接求就行了。

    当然乘的过程中用矩阵乘优化，其实那个矩阵怎么构造的我也不太懂

    

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<iostream>
#define ll long long
using namespace std;
ll read()
{
	ll sum=0,f=1;char x=getchar();
	while(x<'0'||x>'9'){if(x=='-')f=-1;x=getchar();}
	while(x>='0'&&x<='9'){sum=sum*10+x-'0';x=getchar();}
	return sum*f;
}
ll f[6000405],vis[1000405],ine[1000405],don[1000405],len[1000405];
ll n,p,k;
bool flag=0;
struct hh
{
	ll n,m,f[3][3];
	hh(){}
	hh(int x,int y)
	{
		n=x,m=y;
		memset(f,0,sizeof(f));
	}
	
} mul,d,ans,res[1000505];
inline hh operator * (hh a,hh b)
{
	hh c(a.n,b.m);
	for(int i=0;i<a.n;i++)
	   for(int j=0;j<b.m;j++)
	      for(int k=0;k<a.m;k++)
	          (c.f[i][j]+=(a.f[i][k]*b.f[k][j])%p)%=p;
	return c;
}
inline ll ny(ll a,ll b,ll c)
{
	if(a==0)return -1;
	else if((c%a)==0)return c/a;
	ll t=ny(b%a,a,((-c%a)+a)%a);
	if(t==-1)return -1;
	return (t*b+c)/a;
}
inline hh quick(hh a,ll x)
{
	hh sum(a.n,a.m);
	for(int i=0;i<a.n;i++)
	   sum.f[i][i]=1;
	for(;x;x>>=1,a=a*a)
	    if(x&1)sum=sum*a;
	return sum;
}
int yjn()
{
//	freopen("noi2011_rabbit.in","r",stdin);
//	freopen("noi2011_rabbit.out","w",stdout);
	n=read();k=read();p=read();
	f[1]=f[2]=1;
	for(int i=3;i<=6*k;i++)
	{
		f[i]=(f[i-1]+f[i-2])%k;
		if(!vis[f[i]])vis[f[i]]=i;
		if(f[i]==1&&f[i-1]==1)break;
	}
	mul.n=mul.m=d.n=d.m=3;
	ans.n=1;ans.m=3;
	ans.f[0][0]=ans.f[0][2]=1;
	mul.f[0][1]=mul.f[1][0]=mul.f[1][1]=mul.f[2][2]=1;
	d.f[0][0]=d.f[1][1]=d.f[2][2]=1;
	d.f[2][1]=-1;
	for(ll i=1;n;)
	{
		if(!ine[i])ine[i]=ny(i,k,1);
		if(ine[i]==-1)
		{
		    ans=ans*quick(mul,n);
			n=0;
		}
		else
		{
			if(!don[i]||flag)
			{
				don[i]=1;
				if(!vis[ine[i]])
				{
					ans=ans*quick(mul,n);
					n=0;
				}
				else
				{
					len[i]=vis[ine[i]];
					if(n>=len[i])
					{
						n-=len[i];
						res[i]=quick(mul,len[i])*d;
						ans=ans*res[i];
						(i*=f[len[i]-1])%=k;
					}
					else
					{
						ans=ans*quick(mul,n);
						n=0;
					}
				}
			}
			else
			{
				ll cnt=0;
				hh l(3,3);
				l.f[0][0]=l.f[1][1]=l.f[2][2]=1;
				for(ll j=i*f[len[i]-1]%k;j!=i;(j*=f[len[j]-1])%=k)
				    cnt+=len[j],l=l*res[j];
				cnt+=len[i],l=res[i]*l;
				ans=ans*quick(l,n/cnt);
				n%=cnt;
				flag=1;
			}
		}
	}
	cout<<(ans.f[0][1]+p)%p;
}
int qty=yjn();
int main(){;}