2016 ACM/ICPC Asia Regional Shenyang Online 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;
}