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