hdu5895Mathematician QSC

链接：http://acm.hdu.edu.cn/showproblem.php?pid=5895

题意：给定f(0)=0,f(1)=1,f(n)=f(n-2)+2*f(n-1)，g(n)=sigma(f(i)^2){0<=i<=n}。给定多组n,y,x,s，求x^g(n*y)%(s+1)。

分析：写出矩阵递推式，然后就是矩阵快速幂啦。矩阵递推式见代码。

代码：

#include<map>#include<set>#include<cmath>#include<queue>#include<bitset>#include<math.h>#include<vector>#include<string>#include<stdio.h>#include<cstring>#include<iostream>#include<algorithm>#pragma comment(linker, "/STACK:102400000,102400000")using namespace std;const int N=3010;const int M=50010;const int mod=1000000007;const int MOD1=1000000007;const int MOD2=1000000009;const double EPS=0.00000001;typedef long long ll;const ll MOD=1000000007;const int INF=1000000010;const ll MAX=1ll<<55;const double eps=1e-5;const double inf=~0u>>1;const double pi=acos(-1.0);typedef double db;typedef unsigned int uint;typedef unsigned long long ull;struct martrix{    ll x[5][5];}p;ll phi(ll n) {    ll i,ans=n,m=(ll)sqrt(n+0.5);    for (i=2;i<=m;i++)    if (n%i==0) {        ans=ans/i*(i-1);        while (n%i==0) n/=i;    }    if (n>1) ans=ans/n*(n-1);    return ans;}martrix mul(martrix a,martrix b,ll mo) {    martrix ret;    memset(ret.x,0,sizeof(ret.x));    for (int i=1;i<5;i++)        for (int j=1;j<5;j++)            for (int k=1;k<5;k++) (ret.x[i][j]+=a.x[i][k]*b.x[k][j])%=mo;    return ret;}martrix qpowm(ll x,ll mo) {    martrix ret,now=p;    memset(ret.x,0,sizeof(ret.x));    for (int i=1;i<5;i++) ret.x[i][i]=1;    while (x) {        if (x&1) ret=mul(ret,now,mo);        now=mul(now,now,mo);x>>=1;    }    return ret;}ll getmi(ll n,ll mo) {    if (n<=1) return n;    martrix ans=qpowm(n-1,mo);    return ((ans.x[2][4]+ans.x[4][4])%mo+mo)%mo;}ll qpow(ll x,ll y,ll mo) {    ll ret=1;    while (y) {        if (y&1) (ret*=x)%=mo;        (x*=x)%=mo;y>>=1;    }    return (ret+mo)%mo;}int main(){    /*    矩阵递推式    [f(n-2)^2,f(n-1)^2,f(n-2)*f(n-1),g(n-1)]*[0,1,0,1]                                             [1,4,2,4]                                             [0,4,1,4]                                             [0,0,0,1]    =[f(n-1)^2,f(n)^2,f(n-1)*f(n),g(n)]    */    int i,t;    ll n,y,x,s,z;    p.x[1][1]=0;p.x[1][2]=1;p.x[1][3]=0;p.x[1][4]=1;    p.x[2][1]=1;p.x[2][2]=4;p.x[2][3]=2;p.x[2][4]=4;    p.x[3][1]=0;p.x[3][2]=4;p.x[3][3]=1;p.x[3][4]=4;    p.x[4][1]=0;p.x[4][2]=0;p.x[4][3]=0;p.x[4][4]=1;    scanf("%d", &t);    while (t--) {        scanf("%lld%lld%lld%lld", &n, &y, &x, &s);        z=phi(s+1);        printf("%lld\n", qpow(x,getmi(n*y,z)+z,s+1));    }    return 0;}