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;
}