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,