HDU-4794：Arnold（斐波拉契循环节 二次剩余）

本题我只是个搬运工，主要是抢救补板子，所以自己就没写。https://blog.csdn.net/u013534123/article/details/78058997

题意： 大致题意是给你一个N*N的矩阵，然后告诉你阿诺德变换，即原来坐标为(x,y)的点变换一次后变成((x+y)%N,(x+2y)%N)。然后告诉你阿诺德变换一定能够通过有限次变换之后变换回原本的矩阵，然后让你求这个周期。

思路： 不难发现是个斐波拉契循环，题意就是让我们找fib循环节。 然后就开始套板子了。

    1、把n质因数分解，即n=p1^a1*p2^a2...
    2、分别计算Fibonacci数列模每个的循环节，假设循环节为x1,x2,x3...
那么，现在问题就是如何求数列模的循环节。这里我们有一个定理，
Fibonacci数列模p^a的循环节长度等于G(p)*p^(a-1)，在p的剩余系下，
5是否是p的二次剩余，如果是那么周期就是p-1的因子，否则是2*(p+1)的因子。而验证循环节，用矩阵乘法得到对应一项对比即可。

而虽然我的代码是抄的，但是现在的hdu本题最快的，主要是他们的矩阵乘法可能是裸的，取mod次数也偏多，求因子也是sqrt的，balbla....可以参考如下，对此算法做一些优化。    https://blog.csdn.net/abc13068938939/article/details/53836482 

1，少用mod。

2，矩阵乘法可以稍加优化，主要是 fib比较特殊，满足f(n+m)=f(n-1)*f(m)+f(n)*f(m+1)。

3，因为要多次求因子，所以没必要sqrt，而是预处理出1到sqrt的素数。

#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long LL;
const LL maxn=1e5+5;
LL dp[maxn*10];
LL prime[maxn],s=0;
bool vis[maxn];
void init_prime(){
    for(LL i=2;i<maxn;i++){
        if(!vis[i]){
            prime[s++]=i;
        }
        for(LL j=0;j<s&&i*prime[j]<maxn;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){
                break;
            }
        }
    }
    return ;
}
LL pow_mod(LL a1,LL b1){
    LL ans=1;
    while(b1){
        if(b1&1){
            ans=ans*a1;
        }
        b1>>=1;
        a1*=a1;
    }
    return ans;
}
LL pow_mod2(LL a,LL b,LL p){
    LL ans=1;
    while(b){
        if(b&1){
            ans=ans*a%p;
        }
        b>>=1;
        a=a*a%p;
    }
    return ans;
}
LL gcd(LL a,LL b){
    return b?gcd(b,a%b):a;
}
bool f(LL n,LL p){
    if(pow_mod2(n,(p-1)>>1,p)!=1){
        return false;
    }
    return true;
}
struct matrix{
    LL x1,x2,x3,x4;
};
matrix matrix_a,matrix_b;
matrix M2(matrix aa,matrix bb,LL mod){
    matrix tmp;
    tmp.x1=(aa.x1*bb.x1%mod+aa.x2*bb.x3%mod)%mod;
    tmp.x2=(aa.x1*bb.x2%mod+aa.x2*bb.x4%mod)%mod;
    tmp.x3=(aa.x3*bb.x1%mod+aa.x4*bb.x3%mod)%mod;
    tmp.x4=(aa.x3*bb.x2%mod+aa.x4*bb.x4%mod)%mod;
    return tmp;
}
matrix M(LL n,LL mod){
    matrix a,b;
    a=matrix_a;b=matrix_b;
    while(n){
        if(n&1){
            b=M2(b,a,mod);
        }
        n>>=1;
        a=M2(a,a,mod);
    }
    return b;
}
LL fac[100][2],l,x,fs[1000];
void dfs(LL count,LL step){
    if(step==l){
        fs[x++]=count;
        return ;
    }
    LL sum=1;
    for(LL i=0;i<fac[step][1];i++){
        sum*=fac[step][0];
        dfs(count*sum,step+1);
    }
    dfs(count,step+1);
}
LL solve2(LL p){
    if(p<1e6&&dp[p]){
        return dp[p];
    }
    bool ok=f(5,p);//判断5是否为p的二次剩余
    LL t;
    if(ok){
        t=p-1;
    }
    else{
        t=2*p+2;
    }
    l=0;
    for(LL i=0;i<s;i++){
        if(prime[i]>t/prime[i]){
            break;
        }
        if(t%prime[i]==0){
            LL count=0;
            fac[l][0]=prime[i];
            while(t%prime[i]==0){
                count++;t/=prime[i];
            }
            fac[l++][1]=count;
        }
    }
    if(t>1){
        fac[l][0]=t;
        fac[l++][1]=1;
    }
    x=0;
    dfs(1,0);//求t的因子
    sort(fs,fs+x);
    for(LL i=0;i<x;i++){
        matrix m1=M(fs[i],p);
        if(m1.x1==m1.x4&&m1.x1==1&&m1.x2==m1.x3&&m1.x2==0){
            if(p<1e6){
                dp[p]=fs[i];
            }
            return fs[i];
        }
    }
}
LL solve(LL n){
    LL ans=1,cnt;
    for(LL i=0;i<s;i++){
        if(prime[i]>n/prime[i]){
            break;
        }
        if(n%prime[i]==0){
            LL count=0;
            while(n%prime[i]==0){
                count++;n/=prime[i];
            }
            cnt=pow_mod(prime[i],count-1);
            cnt*=solve2(prime[i]);
            ans=(ans/gcd(ans,cnt))*cnt;
        }
    }
    if(n>1){
        cnt=1;
        cnt*=solve2(n);
        ans=ans/gcd(ans,cnt)*cnt;
    }
    return ans;
}
int main(){
    LL t,n;
    init_prime();
    matrix_a.x1=matrix_a.x2=matrix_a.x3=1;
    matrix_a.x4=0;
    matrix_b.x1=matrix_b.x4=1;
    matrix_b.x2=matrix_b.x3=0;
    dp[2]=3;dp[3]=8;dp[5]=20;
    while(~scanf("%I64u",&n)){
        LL ans=solve(n);
        if(ans%2==0) ans/=2;
        printf("%I64u\n",ans);
    }
    return 0;
}