HDU 3977 求斐波那契循环节

题意：求斐波那契数列模一个数的循环节的长度。

分析过程：首先我们知道fib数列模p如果出现了连续的1，0就意味这着开始循环了，因为接下来的项就是1 1 2 3 5等等。

那么很显然如果在第k位第一次出现了1，0，那么对于以后的1，0都可以表示为k*m。

 

那么，现在我们考虑如果fib数列模p在第pos位第一次出现了0，那么设0前面的那个数为a，则接下来的序列将是a，0，a，

a，2a，3a，5a，8a，....。可以看出a的系数就是一个fib数列，那么我们就可以得到fib(k+i)%p=a*fib(i)%p，其中i满

足0<i<k，所以进一步可以得到fib(i)=[a^j*fib(i-k*j)]%p。

那么我们现在先来说说如何求fib数模一个正整数n的循环节长度：

对于这个问题，我们先对n进行素因子分解，得到：

，然后先对每一个形如p^k的数计算循环节，然后它们

的最小公倍数就是n的循环节长度（当然这里面涉及到CRT等等方面的基础）。那么现在问题就是计算p^k的循环节，这个问题

可以进一步简化成求G(p)*p^(k-1)。其中G(p)表示fib数列模素数p的循环节长度，所以现在的问题是如何求fib数列模一个

小于10^6的素数p的循环节长度。

求fib数列模p(p是素数)的最小循环节方法：

暴力枚举fib[i]%p==0的最小的i，然后记录pos=i+1，设a为fib[i]%p==0的前一位数，即a=fib[i-1]

那么我们知道fib数列模p的最小循环节长度一定是pos*x，那么也就是说现在要求一个最小的数x，满足

,

求出x后，那么问题就解决了，剩下的就是合并等等。

#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define maxn 1000005
bool isprime[maxn];
typedef long long ll;
ll prime[maxn],nprime;
ll gcd(ll a,ll b)
{
return b?gcd(b,a%b):a;
}
void getprime()
{
ll i,j;
memset(isprime,1,sizeof(isprime));
nprime=0;
for(i=2; i<maxn; i++)
if(isprime[i])
{
prime[nprime++]=i;
for(j=i*i; j<maxn; j+=i) isprime[j]=0;
}
}
ll factor[100][2],tol;
void findfac(ll n)
{
ll x=n,l=(ll)sqrt((double)n);
tol=0,memset(factor,0,sizeof(factor));
for(int i=0; prime[i]<=l; i++)
if(x%prime[i]==0)
{
factor[tol][0]=prime[i];
while(x%prime[i]==0) factor[tol][1]++,x/=prime[i];
tol++;
}
if(x>1) factor[tol][0]=x,factor[tol++][1]++;
}
ll exp_mod(ll a,ll b,ll c)
{
ll ret=1;
while(b)
{
if(b&1) ret=ret*a%c;
b>>=1,a=a*a%c;
}
return ret;
}
ll getPrimeLoop(ll p)//求一个素数的循环节
{
ll pos=3,f1=1,f2=1,f3=2%p,k=1e9,l=(ll)sqrt((double)p-1);
while(f3) f1=f2,f2=f3,f3=(f1+f2)%p,pos++;//找到第一个值是0的点
for(ll i=1; i<=l; i++)
if((p-1)%i==0)
{
if(exp_mod(f2,(p-1)/i,p)==1) k=min(k,(p-1)/i);
if(exp_mod(f2,i,p)==1) k=min(k,i);
}
return pos*k;
}
ll solve(ll p,ll k)//求一个素数的k次方的循环节
{
ll ans=getPrimeLoop(p);
for(int i=0; i<k-1; i++) ans*=p;
return ans;
}
int main()
{
int t,ca=0;
ll n,ans;
getprime();
scanf("%d",&t);
while(t--)
{
scanf("%I64d",&n);
findfac(n);
ans=1;
ll temp;
for(int i=0; i<tol; i++)
temp=solve(factor[i][0],factor[i][1]),ans=ans/gcd(ans,temp)*temp;
printf("Case #%d: %I64d\n",++ca,ans);
}
return 0;
}