51nod1195斐波那契数列的循环节

求 Fib 数模 n 的循环节：
1. 对 n 做因数分解：
n=p1^e1 * p2^e2 * … * pt^et;
2. 求出每个素数 pi 对应 Fib 数模 pi 的循环节mi0 ，则 pi^ei 对应的 Fib 数模 pi^ei 的 循环节 mi=mi0 * pi^(ei-1);
3. Fib 数模 n 的循环节就等于 lcm(mi)。
关键在于如何求Fib 数模素数 pi 的循环节mi0，有以下结论，如果 pi 是5的二次剩余， mi0 是 pi-1 的约数； 如果 pi 不是，则 mi0 是 2(pi+1) 的约数。

如果想了解上述方法证明过程可以看一下这篇文章
Charles W. Campbell II. The Period of the Fibonacci Sequence Modulo j. 2007.

以下是扩展
广义斐波那契数列

但这道题数据量非常大，很有可能会超时，所以需要一些技巧。
1. 将取模的操作尽量用加减法代替：
比如说 k = (a%p + b%p)%p 就可以写成
if((k = a%p + b%p)>=p) k -= p;
2. 筛出一些质数分解n的复杂度可以由O(sqrt(n))变成O(sqrt(n)/ln(n)):
3. 算模质数p意义下的循环节时，由于斐波那契数列的”循环节的倍数”次项在模意义下相等，所以可以用类似计算欧拉函数的形式将时间复杂度由O(sqrt(n))变成不到O(log(n)log(n))：
首先： 设 L 是模质数p意义下的循环节； 如果 p 是5的二次剩余，M = L 是 p-1 的约数，否则 M = L 是 2(p+1) 的约数。
可以设 L = p1^k1 * p2^k2 * … * pt^kt
M= p1^e1 * p2^e2 * … * pt^et * … * ps^es
由于如果 L|k ，则 k 也是循环节（不一定是最小，但求循环节默认是求最小） ，所以检查 M/(p1^j) 是否满足 (用 4 改进后的快速幂求) f(p1^j)=0 , f(p1^j +1) =1，如果 j 满足而 j - 1 不满足，则一定有 k1=j。依次类推可以将 k1,k2,…,kt求出。最差的情况是要试 (e1+e2+…+es+s) 次，为 log(M)级别，检查时也是 (log) 级别（或者更小）。
4. 斐波那契数列是线性递推数列，采用的公式Fib(n+m) = Fib(n-1) * Fib(m)+Fib(n) * Fib(m+1)：
A={{1,1},{1,0}}， {Fib(n+1),Fib(n)} = A^n * {Fib(1),Fib(0)}。
可以发现 其实 A^n = {{Fib(n+1),Fib(n)},{Fib(n),Fib(n-1)}} = {{Fib(n) + Fib(n-1),Fib(n) },{Fib(n),Fib(n-1)}} ，也就是说每次快速幂做矩阵乘法时没必要像普通矩阵乘法那样算， 比如算 A^(n+m) 本来是 由矩阵 A^n 和 A^m 相乘，现在就可以简化为 求 Fib(n+m) = Fib(n-1) * Fib(m)+Fib(n) * Fib(m+1) 和 Fib(n+m-1) = Fib(n-1) * Fib(m-1)+Fib(n) * Fib(m)，（这些值都在 矩阵 A^n 和 A^m 中） 然后 A^(n+m) ={{Fib(n+m) + Fib(n+m-1),Fib(n+m) },{Fib(n+m),Fib(n+m-1)}}
原先利用2 * 2的矩阵进行快速幂，每次的矩阵乘法是8次加法、8次乘法，但是使用这种改进后，3次加法、4次乘法，可以节省一半时间。
5. 打表（我打了1000以内的素数的情况）减少小质数情况的计算，其实最后这个技巧不加也已经能过了。

查阅博客 http://blog.csdn.net/acdreamers/article/details/10983813
http://blog.csdn.net/zhuangmezhuang/article/details/52627308
并感谢糖老师的讲解。

（这里打表不是必要的）

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#include <stack>
#include <vector>
#include <cstring>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;


int prime[3440];
bool notprime[32000];
void getPrime(void)
{
    ms(notprime);
    for(int i=2;i<32000;i++)
    {
        if(!notprime[i]) prime[++prime[0]]=i;
        for(int j=1;j<=prime[0]&&prime[j]*i<32000;j++)
        {
            notprime[prime[j]*i]=true;
            if(i%prime[j]==0) break;
        }
    }
}

int factor[30][2];
int cnt;
void getFactors(int n)
{
    cnt=0;
    int tmp=n;
    for(int i=1;prime[i]<=tmp/prime[i];i++)
    {
        factor[cnt][1]=0;
        if(tmp%prime[i]==0){
            factor[cnt][0]=prime[i];
            while(tmp%prime[i]==0){
                factor[cnt][1]++;
                tmp/=prime[i];
            }
            cnt++;
        }
    }
    if(tmp!=1){
        factor[cnt][0]=tmp;
        factor[cnt++][1]=1;
    }
}

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

inline void multiply(LL a[2][2],LL b[2][2],int p)
{
    LL tmp01,tmp11;
    tmp01=a[0][0]*b[0][1]%p+a[0][1]*b[1][1]%p;
    if(tmp01>=p) tmp01-=p;
    tmp11=a[1][0]*b[0][1]%p+a[1][1]*b[1][1]%p;
    if(tmp11>=p) tmp11-=p;  
    if((a[0][0]=tmp01+tmp11)>=p) a[0][0]-=p;
    a[0][1]=a[1][0]=tmp01;
    a[1][1]=tmp11;
}
bool check(int n,int p)
{
    LL a[2][2]={{1,1},{1,0}},res[2][2]={{1,0},{0,1}};
    while(n){
        if(n&1) multiply(res,a,p);
        n>>=1;
        multiply(a,a,p);
    }
    return res[0][0]==1&&res[1][0]==0;
}

int biao[1001]={0,3,8,20,16,10,28,36,18,48,14,30,76,40,88,32,108,58,60,136,70,148,78,168,44,196,50,208,72,108,76,256,130,276,46,148,50,316,328,336,348,178,90,190,388,396,22,42,448,456,114,52,238,240,250,516,176,268,270,556,56,568,588,88,310,628,636,110,676,232,174,236,358,736,748,378,768,388,796,200,408,418,84,430,868,438,888,448,916,46,928,936,478,976,490,498,1008,254,26,1048,90,1096,124,376,568,570,1156,1176,1188,598,600,1216,1228,1236,206,630,640,1288,1296,1308,658,220,1348,452,1368,138,700,118,718,1456,1468,738,496,750,1516,380,192,1548,1576,228,202,270,820,1648,1656,276,838,1708,1716,78,1728,1756,176,1768,1776,1816,70,102,928,1876,470,1896,212,176,970,652,1968,198,1996,126,2028,1018,510,206,2068,1038,262,1050,530,2128,356,128,1090,2188,732,96,554,2236,2248,1128,230,2308,2328,1170,1180,2376,2388,1200,2428,812,816,614,410,2476,624,1258,852,426,2568,322,430,2596,1300,2608,872,1318,1320,2656,680,2736,2748,460,1398,704,2848,168,1428,2868,1438,2896,1450,2908,1458,490,1480,424,2976,1488,2988,1498,302,1016,1530,3088,774,1036,1558,3136,1570,526,3168,68,160,3216,1608,3228,1618,810,3256,3276,3316,3328,3336,834,3388,3396,566,854,430,3448,1156,870,3496,3508,1758,3556,3568,3576,894,1800,1810,1216,1830,3696,930,3736,374,3748,1252,1878,1888,1900,3816,1276,1930,3868,1948,390,1316,1978,3976,3988,3996,666,4008,2010,4036,1352,1014,2038,4108,4128,1034,130,4168,4176,1044,2098,2110,4228,2128,2130,4276,2140,4288,4308,80,198,4408,64,4428,148,1492,746,1496,750,1512,324,4548,760,4576,4588,4596,2308,2310,1556,2338,2340,4696,2350,4716,790,4756,2380,4768,398,4788,2398,2410,124,4848,4876,1220,1632,2458,4936,4948,4956,5008,120,2530,2538,5088,2548,2550,5116,2578,518,5188,2608,5236,1310,5268,5296,5316,886,1776,2670,5356,5368,1792,896,5388,2698,5416,2710,5428,2718,682,390,2740,916,1836,5536,1852,164,2790,5596,1400,5608,2818,5668,5676,5688,2850,5716,1430,2878,5776,5796,5808,2908,5836,5856,2938,5908,5916,5928,424,2970,2998,100,3010,3018,864,6076,160,762,3060,6136,162,6168,1544,3108,3118,520,6276,6328,192,1584,3180,6376,3190,6408,3208,6436,3220,3228,650,6508,6516,3258,1090,194,1650,6616,6628,3318,2216,416,1110,6688,2232,3358,3360,3370,6748,1694,1130,2272,6828,6868,3448,532,1730,6928,2312,102,3490,1166,3510,7036,7056,588,2356,3538,3540,7096,2372,3558,34,1790,7168,7188,7216,7228,2412,2416,1210,7276,7288,3658,3670,7348,7356,3690,7396,3700,3708,3718,7456,7468,534,3760,2512,3768,3778,7588,7596,7608,1910,7648,7668,7696,770,7708,7728,7756,1940,3888,7816,3910,7836,1306,7848,982,3930,7888,2632,7936,1994,2000,8008,8016,892,4018,134,8056,2024,1350,8116,8148,4078,4090,8188,4098,4110,8256,4128,2756,4138,8308,1188,1386,8356,2100,4210,8436,4218,2114,1410,4240,8488,8508,4258,1420,4270,8548,1224,1072,8596,8656,8676,4338,4348,8716,8728,324,878,8796,4408,4420,8848,222,8896,4450,8916,8928,128,8968,2996,9016,244,9036,4518,9048,3032,2274,760,9136,3056,4590,9196,9208,462,3092,4638,3096,1162,930,9316,9328,9348,4678,938,3136,4720,9448,2364,1052,950,1586,9568,9576,2394,9588,4798,4800,9628,9636,690,810,974,9756,4888,9808,1636,4918,4930,9868,9876,9888,4950,9916,9936,184,9948,9976,9988,4998,10008,1252,5010,2510,10048,5038,5050,1686,10156,254,10176,5098,2550,10216,10228,5118,1144,10308,10336,5170,5178,5188,10396,2604,10456,5230,10468,3492,2630,10548,5278,528,3532,10608,5308,10648,3556,10696,5350,1076,10776,10788,5398,10816,10828,516,5418,5430,10876,1360,10888,5448,5470,10956,5478,10968,5500,11008,11016,5518,5520,11056,1106,11116,11128,5568,11148,1860,1118,11248,5638,376,11296,5650,1028,11316,5658,2834,11368,1422,3796,1900,5710,11436,604,410,11488,5748,54,11568,5790,1160,3872,11628,1940,11656,1946,3896,2924,5850,11716,2930,11736,978,5878,392,3932,11808,11848,11856,5938,11908,230,11976,12016,1202,6028,12076,12088,12096,12108,12136,12148,2026,6088,870,1220,12228,6120,6130,12268,12288,6150,536,12348,12396,6198,12408,1242,12436,3110,2076,12496,12516,432,6268,1254,12556,12576,6298,6300,6310,12636,12648,904,12676,12688,12708,6358,1590,12736,12748,6378,3194,1828,3210,12856,6448,6450,2156,12948,3240,1298,3260,3264,13096,6550,13108,4376,6568,6570,13156,470,6598,13216,2206,13276,13308,6658,3330,284,742,6688,2230,6700,13408,86,6718,1924,4492,3380,13528,6778,2260,6790,13588,13608,13648,13656,6828,13668,684,13716,13728,3434,458,13768,6898,13816,6910,4612,13896,3474,6958,1740,13936,6970,13956,1552,2330,13996,7000,4676,7018,2008,7038,14088,14116,1178,7078,4736,7108,7120,14256,7128,7150,7158,14356,4792,14388,14416,1442,14428,2406,3614,14476,14488,4832,14508,14568,14596,1624,3654,1830,7330,14668,3674,7350,3684,14788,7410,14836,2124,1490,14916,7458,14956,7480,14976,832,7498,15016,5012,15048,7528,15076,3770,15096,3774,7558,7560,15148,15156,15168,3794,7590,15208,15216,7620,2546,15288,478,2556,15348,7680,15376,7690,7698,15408,15436,15448,15456,7740,15508,5172,7758,7788,15588,5212,5216,7828,1960,5236,2248,15748,15756,202,1752,7900,15816,7918};
int fnd(int n)
{
    int l=0,m,r=1001;
    while(l<r){
        m=(l+r)>>1;
        if(prime[m]<n) l=m+1;
        else if(prime[m]>n) r=m;
        else return biao[m];
    }
}
int main(int argc, char const *argv[])
{
    int _,_i=0;
    scanf("%d",&_);
    getPrime();
    while(++_i<=_){
        int p;LL ans=1ll;
        scanf("%d",&p);
        getFactors(p);
        for(int i=0;i<cnt;i++)
        {
            LL rd=1ll;
            int &pr=factor[i][0];
            if(pr<=prime[1000]) rd=fnd(pr);
            else{
                int m;
                if(pr%5==1||pr%5==4) m=pr-1;
                else m=(pr+1)<<1;
                int tmp=m;
                for(int i=1;prime[i]<=tmp/prime[i];i++)
                    if(m%prime[i]==0){
                        while(check(m/prime[i],pr))
                            m/=prime[i];
                        while(tmp%prime[i]==0) tmp/=prime[i];
                    }
                if(tmp!=1&&check(m/tmp,pr)) m/=tmp;
                rd=m;
            }
            for(int k=1;k<factor[i][1];k++)
                rd*=pr;
            ans=ans/gcd(ans,rd)*rd;
        }
        printf("%I64d\n",ans );
    }
    return 0;
}