HDU 5451——递推式&&循环节

题意

设 $y = (5+2sqrt 6)^{1+2^x}$，给出 $x, M$（$0leq x leq 2^{32}, M leq 46337$），求 $[y]\%M$.

分析

由通项推递推式？？

设 $A_n = (5 + 2sqrt 6)^n, B_n = (5 - 2sqrt 6)^n,C_n = A_n + B_n$,

显然 $C_n$ 是整数，且 $B_n$ 是小于1的，所以答案就是 $C_n - 1$.

通过推导：

$C_n = A_n + B_n = (5+2sqrt6)^n +  (5-2sqrt6)^n$

$C_{n+1} = A_{n+1} + B_{n+1} = (5+2sqrt6)(5+2sqrt6)^n +  (5-2sqrt6)(5-2sqrt6)^n$

$C_{n+2} = A_{n+2} + B_{n+2} = (49+20sqrt6)(5+2sqrt6)^n +  (49-20sqrt6)(5-2sqrt6)^n$，

观察得 $C_{n+2} = 10C_{n+1} - C_n$

写成矩阵快速幂的形式，即

$$egin{bmatrix} C_n\  C_{n-1} end{bmatrix} = {egin{bmatrix} 10 & -1\  1 & 0 end{bmatrix}}^{n-1}egin{bmatrix} C_1\  C_0 end{bmatrix}$$

幂太大，直接用快速幂肯定TLE。

我们可以找循环节，

由于模和转移矩阵是确定的，可以暴力打表找规律。

也可以用结论，因为 $M$ 为素数，存在循环节 $(M+1)(M-1)$。这不一定是最小的循环节，可以枚举其因子找到最小的循环节。

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N=1000000+10;
ll x0,x1,a,b,n,mod;
char s[N];

const int maxn = 100000 + 10;       //p最为2e9，不会有两个超过1e5的质因数
int prime[maxn], pcnt;            //prime[i]表示第i个素数
bool is_prime[maxn + 1];    //is_prime[i]为true表示i是素数
int sieve(int n)
{
    int cnt = 0;
    for (int i = 0; i <= n; i++)  is_prime[i] = true;
    is_prime[0] = is_prime[1] = false;
    for (ll i = 2; i <= n; i++)
    {
        if (is_prime[i])
        {
            prime[cnt++] = i;
            for (ll j = i * i; j <= n; j += i)  is_prime[j] = false;  //i * i可能爆int
        }
    }
    return cnt;
}

ll solve(ll x){
    ll ans1=1,ans2=1,xx=x;
    for(int i=0;i<pcnt;i++){
        if(1ll*prime[i]*prime[i]>x) break;
        if(x%prime[i]==0){
            ans1*=(prime[i]-1)*(prime[i]+1);
            ans2*=prime[i];
            while(x%prime[i]==0) x/=prime[i];
        }
    }
    if(x>1){
        ans1*=(x-1)*(x+1);
        ans2*=x;
    }
    return xx/ans2*ans1;
}
ll qmul(ll x,ll y,ll p){        //快速乘
    x%=p;
    y%=p;
    ll ans=0;
    while(y){
        if(y&1){
            ans+=x;
            if(ans>=p) ans-=p;  //这样写不能有负数
        }
        x<<=1;
        if(x>=p) x-=p;
        y>>=1;
    }
    return ans;
}

struct Mat{
    int r,c;
    ll m[3][3];
    Mat(){
        memset(m,0,sizeof(m));
    }
};

Mat mmul(Mat x,Mat y,ll p){
    Mat ans;
    ans.r=x.r;
    ans.c=y.c;
    for(int i=0;i<x.r;i++)
        for(int k=0;k<x.c;k++)
            for(int j=0;j<y.c;j++){
                ans.m[i][j]+=qmul(x.m[i][k],y.m[k][j],p);
                if(ans.m[i][j]>=p) ans.m[i][j]-=p;
            }
    return ans;
}
Mat mpow(Mat x,ll y,ll p){
    Mat ans;
    ans.r=x.r;
    ans.c=x.c;
    for(int i=0;i<ans.c;i++) ans.m[i][i]=1;
    while(y){
        if(y&1) ans=mmul(ans,x,p);
        x=mmul(x,x,p);
        y>>=1;
    }
    return ans;
}
int main(){
    pcnt = sieve(100000);
    while(scanf("%lld%lld%lld%lld",&x0,&x1,&a,&b) == 4){
        scanf("%s%lld",s,&mod);
        ll lop=solve(mod);  //循环节长度
        n=0;
        int lens=strlen(s);
        for(int i=0;i<lens;i++){
            n=qmul(n,10,lop)+s[i]-'0';
            if(n>=lop) n-=lop;
        }
        Mat A,T;
        A.r=2; A.c=1;
        A.m[0][0]=x1; A.m[1][0]=x0;
        T.r=2; T.c=2;
        T.m[0][0]=a; T.m[0][1]=b; T.m[1][0]=1;
        if(n>1){
            T=mpow(T,n-1,mod);
            A=mmul(T,A,mod);
        }
        printf("%lld
",A.m[0][0]);
    }
    return 0;
}

参考链接：https://www.cnblogs.com/addf/p/4834108.html