luoguP4000 斐波那契数列

Posted 2020-12-27 sssy

题目链接

luoguP4000 斐波那契数列

题解

根据这个东西
https://www.cnblogs.com/sssy/p/9418732.html
我们可以找出%p意义下的循环节
然后就可以做了
人傻，自带，大，常数

代码

#include<bits/stdc++.h> 
using namespace std; 
#define LL long long 
const LL maxn = 1000007; 
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;
        } 
    } 
} 
LL pow_mod(LL a1,LL b1){
    LL ret = 1;
    for(;b1;b1 >>= 1,a1 *= a1) 
        if(b1 & 1) ret = ret * a1; 
    return ret; 
} 
LL pow_mod2(LL a,LL b,LL p) { 
    LL ret = 1; 
    for(;b;b >>= 1,a = a * a % p)  
        if(b & 1) ret = ret * a % p; 
    return ret; 
} 
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;
    for(;n;n >>= 1,a = M2(a,a,mod))  
        if(n & 1) b = M2(b,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 get_M(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; 
} 
char Ss[30000007]; 
struct bign { 
    int z[30000007],l; 
    void init() { 
        memset(z,0,sizeof(z)); 
        scanf("%s",Ss + 1); 
        l = strlen(Ss + 1); 
        for(int i = 1;i <= l;i ++) 
            z[i] = Ss[l - i + 1] - '0'; 
    } 
    LL operator % (const long long & a) const { 
        LL b = 0; 
        for (int i = l;i >= 1;i --) 
            b = (b * 10 + z[i]) % a; 
        return b; 
    } 
}z;  
LL m1; 
struct Matrix { 
    LL a[3][3]; 
    Matrix () { memset(a,0,sizeof a); }     
    Matrix  operator * (const Matrix & p) const { 
        Matrix ret; 
        for(int i = 0;i <= 1;++ i) 
            for(int j = 0;j <= 1;++ j) 
                for(int k = 0;k <= 1;++ k) 
                    ret.a[i][j] = (ret.a[i][j] + ( a[i][k] * p.a[k][j] ) % m1 ) % m1; 
        return ret; 
    } 
}; 
LL solve(LL x) { 
    Matrix p,q; 
    p.a[0][0] = 1; p.a[0][1] = 1; p.a[1][0] = 1; 
    q.a[0][1] = 1;  
    for(;x;x >>= 1,p = p * p) 
        if(x & 1) q = q * p; 
    return q.a[0][0]; 
} 
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; 
    z.init(); scanf("%lld",&n); m1 = n; 
    n = get_M(n); 
    //printf("%lld
",n % m1); 
    n = z % n; 
    printf("%lld
",solve(n)); 
    return 0; 
}