首页 > 代码库 > HDU 5895 矩阵快速幂+高次幂取模

HDU 5895 矩阵快速幂+高次幂取模

HDU 5895 Mathematician QSC

题意:已知f(n)=2*f(n-1)+f(n-2), g(n)=∑f(i)²(0<=i<=n), 给出n,x,y,s, 求x^(g(n*y))%(s+1);

思路:OEIS查到了g(n)=f(n)*f(n+1)/2, f(n)可以用矩阵快速幂求得, 有一个定理可以用于高次幂取模 x^n %k=x^(n%phi(k)+phi(k)) %k, 此处phi(x)为欧拉函数,但是在对幂次取模时存在一个除2,

又因为(a/b)%k=(a%bk)/b,所以这个问题得以解决(这个方法和逆元有点分不清, 还得好好看看).

技术分享
#include<cstdio>#include<cstdlib>#include<cstring>#include<cmath>#include<iostream>#include<algorithm>using namespace std;typedef long long LL;LL n,x,y,s1,S;LL euler(LL n){    LL res=n, a=n;    for(LL i=2;i*i<=n;i++){        if(a%i==0){            res=res/i*(i-1);            while(a%i==0) a/=i;        }    }    if(a>1) res=res/a*(a-1);    return res;}LL pow_mod(LL a,LL n){    LL t=a, res=1;    while(n){        if(n&1) res=(res*t)%(S+1);        n/=2;        t=(t*t)%(S+1);    }    return res;}struct Mat{    LL a[2][2];    void init(){        a[0][0]=2;        a[1][0]=a[0][1]=1;        a[1][1]=0;    }};Mat operator *(Mat a,Mat b){    Mat c;    for(LL i=0;i<2;i++)    for(LL j=0;j<2;j++){        c.a[i][j]=0;        for(LL k=0;k<2;k++)            c.a[i][j]+=a.a[i][k]*b.a[k][j];        c.a[i][j]=c.a[i][j]%s1;    }    return c;}Mat operator ^(Mat p,LL k){    Mat ans; ans.init();    while(k){        if(k&1)            ans=ans*p;        k/=2;        p=p*p;    }    return ans;}int main(){    LL t; cin>>t;    while(t--){        cin>>n>>y>>x>>S;        s1=2*euler(S+1);        if(n==-1) break;        if(n==0){            cout<<1<<endl;            continue;        }        Mat s; s.init();        s=s^(n*y-1);        LL u=(s.a[0][0]*s.a[1][0]);        u=(u%s1+s1)/2;        LL ans=pow_mod(x,u);        cout<<ans<<endl;    }    return 0;}
Psong

 

HDU 5895 矩阵快速幂+高次幂取模