HDU 5895 Mathematician QSC(矩阵乘法+循环节降幂+除法取模小技巧+快速幂)

传送门：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\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=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;
}