传送门：HDU 5895 Mathematician QSC
这是一篇很好的题解，我想讲的他基本都讲了http://blog.csdn.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 ",a) #define ptlld(a) printf("%I64d ",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=4*atan(1.0); using namespace std; typedef long long ll; template <class T> inline void read(T &n) { char c; int flag = 1; for (c = getchar(); !(c >= '0' && c <= '9' || c == '-'); c = getchar()); if (c == '-') flag = -1, n = 0; else n = c - '0'; for (c = getchar(); c >= '0' && c <= '9'; c = getchar()) n = n * 10 + c - '0'; n *= flag; } ll Pow(ll base, ll n, ll mo) { ll res=1; while(n) { if(n&1) res=res*base%mo; n>>=1; base=base*base%mo; } return res; } //*************************** ll n,y,x,s; const ll mod=1000000007; ll modp,modq; const int maxn=2; ll euler(ll nn) { ll res=nn,a=nn; for(ll i=2;i*i<=a;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; } struct matrix { ll mat[maxn][maxn]; matrix operator*(const matrix & rhs)const { matrix ans; rep(i,0,maxn-1) rep(j,0,maxn-1) ans.mat[i][j]=0; rep(i,0,maxn-1) rep(j,0,maxn-1) rep(k,0,maxn-1) 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,0,maxn-1) rep(j,0,maxn-1) res.mat[i][j]=(i==j); while(k) { if(k&1) res=res*rhs; rhs=rhs*rhs; k>>=1; } return res; } }xx; int main() { #ifndef ONLINE_JUDGE freopen("in.txt","r",stdin); #endif int T_T; scanf("%d",&T_T); for(int kase=1;kase<=T_T;kase++) { read(n),read(y),read(x),read(s); modp=euler(s+1)*2; modq=s+1; xx.mat[0][0]=2,xx.mat[0][1]=1,xx.mat[1][0]=1,xx.mat[1][1]=0; matrix temp=xx^(n*y); ll fn1=temp.mat[0][0]; ll fn=temp.mat[1][0]; ll gn=fn*fn1%modp/2; ll ans=Pow(x,gn,modq); ptlld(ans); } return 0; }