HDU3977(斐波那契据数列模n的循环节长度)

HDU3977(斐波那契数列模n的循环节长度)

题目:http://acm.hdu.edu.cn/showproblem.php?pid=3977

 

题意:求斐波那契数列模p的循环节长度,注意p最大是2*10^9,但是它的素因子小于10^6。

 

分析过程:首先我们知道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进行素因子分解,得到:HDU3977(斐波那契据数列模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,满足HDU3977(斐波那契据数列模n的循环节长度),

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

#include <iostream>
#include <string.h>
#include <algorithm>
#include <stdio.h>
#include <math.h>

using namespace std;
typedef long long LL;
const int N=1000005;

bool prime[N];
int p[N];
int pri[N];
int num[N];
int f[N];
int fac[N];
int arr[N];
int k,cnt;

void isprime()
{
    k=0;
    int i,j;
    memset(prime,true,sizeof(prime));
    for(i=2; i<N; i++)
    {
        if(prime[i])
        {
            p[k++]=i;
            for(j=i+i; j<N; j+=i)
            {
                prime[j]=false;
            }
        }
    }
}

void Solve(LL n)
{
    cnt=0;
    int ct=0;
    int t=(int)sqrt(n*1.0);
    for(int i=0; p[i]<=t; i++)
    {
        if(n%p[i]==0)
        {
            ct=0;
            pri[cnt]=p[i];
            while(n%p[i]==0)
            {
                ct++;
                n/=p[i];
            }
            num[cnt]=ct;
            cnt++;
        }
    }
    if(n>1)
    {
        pri[cnt]=n;
        num[cnt]=1;
        cnt++;
    }
}

LL gcd(LL a,LL b)
{
    return b? gcd(b,a%b):a;
}

LL quick_mod(LL a,LL b,LL m)
{
    LL ans=1;
    a%=m;
    while(b)
    {
        if(b&1)
        {
            ans=ans*a%m;
            b--;
        }
        b>>=1;
        a=a*a%m;
    }
    return ans;
}

int main()
{
    isprime();
    int T,a,pos,tt=1;
    LL n;
    cin>>T;
    while(T--)
    {
        cin>>n;
        printf("Case #%d: ",tt++);
        Solve(n);
        pos=0;
        for(int k=0; k<cnt; k++)
        {
            f[0]=f[1]=1;
            for(int i=2;; i++)
            {
                f[i]=(f[i-1]+f[i-2])%pri[k];
                if(f[i]==0)
                {
                    a=f[i-1];
                    pos=i+1;
                    break;
                }
            }
            int cv=0;
            int tmp=pri[k]-1;
            int t=(int)sqrt(tmp*1.0);
            for(int i=1; i<=t; i++)
            {
                if(tmp%i==0)
                {
                    if(tmp/i==i)
                        fac[cv++]=i;
                    else
                    {
                        fac[cv++]=i;
                        fac[cv++]=tmp/i;
                    }
                }
            }
            int record=0;
            sort(fac,fac+cv);
            for(int i=0; i<cv; i++)
            {
                if(quick_mod(a,fac[i],pri[k])==1)
                {
                    record=fac[i];
                    break;
                }
            }
            LL ans=record*pos;
            for(int i=1; i<num[k]; i++)
                ans*=pri[k];
            arr[k]=ans;
        }
        LL ret=1;
        for(int i=0; i<cnt; i++)
            ret=ret/gcd(ret,arr[i])*arr[i];
        cout<<ret<<endl;
    }
    return 0;
}


 

其实本题还有一种方法,它是基于下面这两个定理的(证明暂时省略):

定理一:对于p>5,如果5是p的平方剩余,本原0位的位置出现在(p-1)及其因子位,如果是平方非剩余,则出现在2*(p+1)及

其因子位。

 

定理二:如果5是p的平方剩余,循环节的长度是p-1的因子,如果不是,循环节长度是2*p+2的因子。

 

那么我们可以这样做HDU3977了,找循环节的时候就先分两种情况,然后找出因子,再从小到大枚举因子,找最小的x符合

Fib[x-1]=1,Fib[x]=0,然后就是循环节了,当然这个最好用矩阵快速幂,不然数据大可能TLE。