HDU 5895 Mathematician QSC（欧拉定理推广）

来源：互联网发布：淘宝卖品牌要授权吗编辑：程序博客网时间：2024/05/21 15:12

分析：

附上链接：icpccamp
当 x>φ(m)时，ax≡aφ(m)+xmodφ(m)(modφ(m))总是对的。不需要互质的条件。
然后高一个矩阵快速幂就搞定了（是不是矩阵快速幂可以搞定所有序列的通项？）

代码：

#include <bits/stdc++.h>using namespace std;typedef long long LL;#define FOR(i,x,y)  for(int i = x;i < y;++ i)const int maxn = 1000010;int prime[maxn];bool check[maxn];void Mobius(){    memset(check,false,sizeof(check));    prime[0] = 0;    FOR(i,2,maxn){        if(!check[i])   {prime[++prime[0]] = i;}        FOR(j,1,prime[0]+1){            if(i*prime[j] > maxn)  break;            check[i*prime[j]] = true;            if(i% prime[j] == 0)    break;        }    }}long long quickpow(long long a,int n,int m){    long long ans=1;    while(n){        if(n&1) ans = (ans*a)%m;        a = (a*a)%m;        n>>=1;    }    return ans;}int phi(int n){    int nn = (int)sqrt(n);    int ans = 1;    for(int i = 1;i <= prime[0];++ i){        if(prime[i] > n && prime[i] > nn)   break;        int cnt = 0;        while(n%prime[i] == 0){            cnt ++;            n /= prime[i];        }        if(cnt == 0)    continue;        ans = ans*quickpow(prime[i],cnt-1,1e9+7)*(prime[i]-1);    }    if(n != 1)  ans = ans*(n-1);    return ans;}int Mod;const int N = 4;struct Matrix{    LL mat[N][N];    ///两个矩阵乘法    Matrix operator *(const Matrix& a) const{        Matrix c;        memset(c.mat,0,sizeof(c.mat));        for(int i=0;i<N;i++){            for(int j=0;j<N;j++){                for(int k=0;k<N;k++){                    c.mat[i][j] += mat[i][k]*a.mat[k][j]%Mod;                    c.mat[i][j] %= Mod;                }            }        }        return c;    }    ///矩阵的幂，用二分的思想写，就是快速幂算法    Matrix operator ^(LL n) const{        Matrix c;        ///初始化为单位矩阵        for(int i=0;i<N;i++){            for(int j=0;j<N;j++){                c.mat[i][j]=(i==j);            }        }        ///就是快速幂        Matrix a = *this;///毕竟直接敲的模板，没运行就写上去了，发现这里少了一个*号        while(n)        {            if(n%2) c=c*a;            a=a*a;            n=n>>1;        }        return c;    }}rhs;int n,y,x,s,a[N];int main(){    //freopen("test.in","r",stdin);    Mobius();    rhs.mat[0][0] = 1; rhs.mat[0][1] = 4; rhs.mat[0][2] = 4; rhs.mat[0][3] = 1;    rhs.mat[1][0] = 0; rhs.mat[1][1] = 4; rhs.mat[1][2] = 4; rhs.mat[1][3] = 1;    rhs.mat[2][0] = 0; rhs.mat[2][1] = 2; rhs.mat[2][2] = 1; rhs.mat[2][3] = 0;    rhs.mat[3][0] = 0; rhs.mat[3][1] = 1; rhs.mat[3][2] = 0; rhs.mat[3][3] = 0;    a[0] = 1; a[1] = 1; a[2] = 0; a[3] = 0;    int T;  scanf("%d",&T);    while(T--){        scanf("%d%d%d%d",&n,&y,&x,&s);        long long ny = (LL)n*y;        if(ny <= 10){            long long res = 1;            LL p = 1,q = 0;            for(int i = 2;i <= ny;++ i){                LL t = 2*p+q;                res += t*t;                q = p;  p = t;            }            printf("%lld\n",quickpow(x,res,s+1));            continue;        }        Mod = phi(s+1);        Matrix res = rhs^(ny-1);        LL t = 0;        for(int j = 0;j < N;++ j){            t += res.mat[0][j]*a[j]%Mod;            t %= Mod;        }        t += Mod;        printf("%lld\n",quickpow(x,t,s+1));    }    return 0;}

0 0