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

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<<"*************************"<<endlconst 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;}

 

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