传送门：HDU 5895 Mathematician QSC
这是一篇很好的题解，我想讲的他基本都讲了http://blog.****.net/queuelovestack/article/details/52577212
【分析】
一开始想简单了,对于a^x mod p这种形式的直接用欧拉定理的数论定理降幂了
结果可想而知,肯定错,因为题目并没有保证gcd(x,s+1)=1,而欧拉定理的数论定理是明确规定的
所以得另谋出路
那么网上提供了一种指数循环节降幂的方法
具体证明可以自行从网上找一找
有了这种降幂的方法之后,我们要分析一下如何求g(n)
由于f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2)
可得,g(n)=f(n)*f(n+1)/2
这个是很好发现的
如果你发现不了的话,可以直接丢到OEIS里搜一下
然后,要求出g(n*y),就需要先求出f(n*y)和f(n*y+1)
这时,我们可以考虑用矩阵乘法
构造矩阵
套一下矩阵快速幂的模板就可以求出f(n*y)和f(n*y+1)
然后要求g(n)还有个除以2的操作,显然除法取模要用逆元
但考虑到2与模数不一定互质,无法用乘法逆元,所以要采用一点小技巧转化一下
这样我们就可以得到简化好的最终的指数部分
这样我们用快速幂就可以求x的幂次对(s+1)取模了
【时间复杂度&&优化】
O(1ogn)
/**************************************************************
Problem:hdu 5895 Mathematician QSC
User: youmi
Language: C++
Result: Accepted
Time:31MS
Memory:1584K
****************************************************************/
//#pragma comment(linker, "/STACK:1024000000,1024000000")
//#include<bits/stdc++.h>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <stack>
#include <set>
#include <sstream>
#include <cmath>
#include <queue>
#include <deque>
#include <string>
#include <vector>
#define zeros(a) memset(a,0,sizeof(a))
#define ones(a) memset(a,-1,sizeof(a))
#define sc(a) scanf("%d",&a)
#define sc2(a,b) scanf("%d%d",&a,&b)
#define sc3(a,b,c) scanf("%d%d%d",&a,&b,&c)
#define scs(a) scanf("%s",a)
#define sclld(a) scanf("%I64d",&a)
#define pt(a) printf("%d\n",a)
#define ptlld(a) printf("%I64d\n",a)
#define rep(i,from,to) for(int i=from;i<=to;i++)
#define irep(i,to,from) for(int i=to;i>=from;i--)
#define Max(a,b) ((a)>(b)?(a):(b))
#define Min(a,b) ((a)<(b)?(a):(b))
#define lson (step<<1)
#define rson (lson+1)
#define eps 1e-6
#define oo 0x3fffffff
#define TEST cout<<"*************************"<<endl
const double pi=*atan(1.0); using namespace std;
typedef long long ll;
template <class T> inline void read(T &n)
{
char c; int flag = ;
for (c = getchar(); !(c >= '' && c <= '' || c == '-'); c = getchar()); if (c == '-') flag = -, n = ; else n = c - '';
for (c = getchar(); c >= '' && c <= ''; c = getchar()) n = n * + c - ''; n *= flag;
}
ll Pow(ll base, ll n, ll mo)
{
ll res=;
while(n)
{
if(n&)
res=res*base%mo;
n>>=;
base=base*base%mo;
}
return res;
}
//*************************** ll n,y,x,s;
const ll mod=;
ll modp,modq;
const int maxn=; ll euler(ll nn)
{
ll res=nn,a=nn;
for(ll i=;i*i<=a;i++){
if(a%i==){
res=res/i*(i-);//先进行除法是为了防止中间数据的溢出
while(a%i==) a/=i;
}
}
if(a>) res=res/a*(a-);
return res;
}
struct matrix
{
ll mat[maxn][maxn];
matrix operator*(const matrix & rhs)const
{
matrix ans;
rep(i,,maxn-)
rep(j,,maxn-)
ans.mat[i][j]=;
rep(i,,maxn-)
rep(j,,maxn-)
rep(k,,maxn-)
ans.mat[i][j]=(ans.mat[i][j]+mat[i][k]*rhs.mat[k][j])%modp;
return ans;
}
matrix operator^(ll k)const
{
matrix rhs=*this;
matrix res;
rep(i,,maxn-)
rep(j,,maxn-)
res.mat[i][j]=(i==j);
while(k)
{
if(k&)
res=res*rhs;
rhs=rhs*rhs;
k>>=;
}
return res;
}
}xx; int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
#endif
int T_T;
scanf("%d",&T_T);
for(int kase=;kase<=T_T;kase++)
{
read(n),read(y),read(x),read(s);
modp=euler(s+)*;
modq=s+;
xx.mat[][]=,xx.mat[][]=,xx.mat[][]=,xx.mat[][]=;
matrix temp=xx^(n*y);
ll fn1=temp.mat[][];
ll fn=temp.mat[][];
ll gn=fn*fn1%modp/;
ll ans=Pow(x,gn,modq);
ptlld(ans);
}
return ;
}