ACdream 1124 喵喵的遗憾

题目http://acdream.info/problem?pid=1124

这道题首先要明白

F(F(F(N)))%P=F(        F (     F ( N) % L ( L ( P ) )        ) %L(P)        ) %P   其中 L（P）是f（n）%p的循环节

之后这道题就转化成了求斐波那契数列的循环节问题

参考这里

http://blog.csdn.net/acdreamers/article/details/10983813

对于一个正整数n，我们求Fib数模n的循环节的长度的方法如下：

    （1）把n素因子分解，即

    （2）分别计算Fib数模每个的循环节长度，假设长度分别是

    （3）那么Fib模n的循环节长度

从上面三个步骤看来，貌似最困难的是第二步，那么我们如何求Fib模的循环节长度呢？

 

     这里有一个优美的定理：Fib数模的最小循环节长度等于，其中表示Fib数模素数的最小循环节长度。可以看出我们现在最重要的就是求

 

 

对于求我们利用如下定理：

 

   如果5是模的二次剩余，那么循环节的的长度是的因子，否则，循环节的长度是的因子。

顺便说一句，对于小于等于5的素数，我们直接特殊判断，loop(2)=3,loop(3)=8,loop(5)=20。

那么我们可以先求出所有的因子，然后用矩阵快速幂来一个一个判断，这样时间复杂度不会很大。

代码如下：

/**

吉林大学
Jilin U

Author:     sinianluoye (JLU_LiChuang)
Date:       2014-07-24
Usage:      ACdream 1124

**/

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>

#define ll long long
#define eps 1e-8
#define ms(x,y) (memset(x,y,sizeof(x)))
#define fr(i,x,y) for(ll i=x;i<=y;i++)

using namespace std;

/****************************fib数模N的循环节*********************************/


struct matrix
{
    ll m[2][2];
    matrix(ll a,ll b,ll c,ll d)
    {
        m[0][0]=a;
        m[0][1]=b;
        m[1][0]=c;
        m[1][1]=d;
    }
    matrix(ll x=1,ll y=1)
    {
        m[0][0]=x;
        m[0][1]=y;
        m[1][0]=1;
        m[1][1]=0;
    }
    matrix operator %= (ll M)
    {
        for(ll i=0;i<2;i++)
            for(ll j=0;j<2;j++)
            m[i][j]%=M;
        return *this;
    }
};

matrix mul(matrix A,matrix B,ll mod)
{
    matrix ret(0,0,0,0);
    fr(i,0,1)fr(j,0,1)fr(k,0,1)
        ret.m[i][j]=(ret.m[i][j]+(A.m[i][k]*B.m[k][j])%mod)%mod;
    return ret;
}

matrix fast_mod(matrix T,ll n,ll mod)
{
    matrix E(1,0,0,1);
    while(n)
    {
        if(n&1)
        {
            E=mul(E,T,mod);
            E%=mod;
        }
        T=mul(T,T,mod);
        T%=mod;
        n>>=1;
    }
    return E;
}

ll fib(ll t1,ll t2,ll a0,ll a1,ll n,ll mod)///t1 t2 递推式 a0 a1 初值
{
    if(n==0)return a0;
    if(n==1)return a1;
    matrix T(t1,t2);
    T=fast_mod(T,n-1,mod);
    return ((T.m[0][0]*a1)%mod+(T.m[0][1]*a0)%mod)%mod;
}

ll fib(ll n,ll mod)
{
    return fib(1,1,1,1,n,mod);
}

ll gcd(ll a,ll b)
{
    while(b)swap(a%=b,b);return a;
}

const ll N = 400005;
const ll NN = 5005;

ll num[NN],pri[NN];
ll fac[NN];
ll cnt,c;

bool prime[N];
ll p[N];
ll k;

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

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;
}

ll legendre(ll a,ll p)
{
    if(quick_mod(a,(p-1)>>1,p)==1) return 1;
    else                           return -1;
}

void Solve(ll n,ll pri[],ll num[])
{
    cnt = 0;
    ll t = (ll)(sqrt(1.0*n)+0.5);

    for(ll i=0; p[i]<=t; i++)
    {
        if(n%p[i]==0)
        {
            ll a = 0;
            pri[cnt] = p[i];
            while(n%p[i]==0)
            {
                a++;
                n /= p[i];
            }
            num[cnt++] = a;
        }
    }
    if(n > 1)
    {
        pri[cnt] = n;
        num[cnt] = 1;
        cnt++;
    }
}

void Work(ll n)
{
    c = 0;
    ll t = (ll)sqrt(1.0*n);
    for(ll i=1; i<=t; i++)
    {
        if(n % i == 0)
        {
            if(i * i == n) fac[c++] = i;
            else
            {
                fac[c++] = i;
                fac[c++] = n / i;
            }
        }
    }
}

ll find_loop(ll n)
{
    Solve(n,pri,num);
    ll ans=1;
    for(ll i=0; i<cnt; i++)
    {
        ll record=1;
        if(pri[i]==2)
            record=3;
        else if(pri[i]==3)
            record=8;
        else if(pri[i]==5)
            record=20;
        else
        {
            if(legendre(5,pri[i])==1)
                Work(pri[i]-1);
            else
                Work(2*(pri[i]+1));
            sort(fac,fac+c);
            for(ll k=0; k<c; k++)
            {
                matrix A;
                matrix a = fast_mod(A,fac[k]-1,pri[i]);
                ll x = (a.m[0][0]%pri[i]+a.m[0][1]%pri[i])%pri[i];
                ll y = (a.m[1][0]%pri[i]+a.m[1][1]%pri[i])%pri[i];
                if(x==1 && y==0)
                {
                    record = fac[k];
                    break;
                }
            }
        }
        for(ll k=1; k<num[i]; k++)
            record *= pri[i];
        ans = ans/gcd(ans,record)*record;
    }
    return ans;
}
/********************************--------------end-------------**************************************/

ll L(ll n){return find_loop(n);}
ll F(ll n,ll m){return fib(n,m);}

int main()
{
    ll T;
    cin>>T;
    isprime();
    while(T--)
    {
        ll n,p;
        scanf("%lld%lld",&n,&p);
        if(p!=1)
        {
            ll mod1=p;
            ll mod2=L(p);
            ll mod3=L(mod2);
            ll f1=F(n,mod3);
            ll f2=F(f1,mod2);
            cout<<F(f2,mod1)<<endl;
        }
        else
            cout<<0<<endl;
    }
}

/*************copyright by sinianluoye (JLU_LiChuang)***********/

注意  不要直接写成=F(        F (     F ( N) % L ( L ( P ) )        ) %L(P)        ) %P

否则会因为重复计算而TLE