HDU 5895 Mathematician QSC(矩阵乘法+循环节降幂+除法取模小技巧+快速幂)——2016 ACM/ICPC Asia Regional Shenyang Online

原文链接：http://blog.csdn.net/queuelovestack/article/details/52577212
解题思路:
【题意】
已知f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2)

给你n,y,x,s的值
求的值
【类型】
矩阵乘法+循环节降幂+除法取模小技巧+快速幂
【分析】
一开始想简单了,对于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)

∴f(0)=0,f(1)=1,f(2)=2,f(3)=5,f(4)=12,f(5)=29,f(6)=70……
∴g(0)=f(0)*f(0)=0=f(0)*f(1)/2
g(1)=g(0)+f(1)*f(1)=1=f(1)*f(2)/2
g(2)=g(1)+f(2)*f(2)=5=f(2)*f(3)/2
g(3)=g(2)+f(3)*f(3)=30=f(3)*f(4)/2
g(4)=g(3)+f(4)*f(4)=174=f(4)*f(5)/2
g(5)=g(4)+f(5)*f(5)=1015=f(5)*f(6)/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(logn)
题目链接：HDU 5895

/*Sherlock and Watson and Adler*/#pragma comment(linker, "/STACK:1024000000,1024000000")#include<stdio.h>#include<string.h>#include<stdlib.h>#include<queue>#include<stack>#include<math.h>#include<vector>#include<map>#include<set>#include<bitset>#include<cmath>#include<complex>#include<string>#include<algorithm>#include<iostream>#define eps 1e-9#define LL long long#define PI acos(-1.0)#define bitnum(a) __builtin_popcount(a)using namespace std;const int N = 2;const int M = 100005;const int inf = 1000000007;//const int mod = 1000000007;typedef struct node{    __int64 a[N][N];    void Init()    {        memset(a,0,sizeof(a));        for(int i=0;i<N;i++)            a[i][i]=1;    }}matrix;int mod;bool flag;matrix mul(matrix a,matrix b)//矩阵乘法{    matrix ans;    for(int i=0;i<N;i++)        for(int j=0;j<N;j++)        {            ans.a[i][j]=0;            for(int k=0;k<N;k++)            {                ans.a[i][j]+=a.a[i][k]*b.a[k][j];                if(ans.a[i][j]>mod)                    ans.a[i][j]=ans.a[i][j]%mod+mod;            }        }    return ans;}matrix pow(matrix a,__int64 n)//求a^n{    matrix ans;    ans.Init();    while(n)    {        if(n%2){            ans=mul(ans,a);        }        n/=2;        a=mul(a,a);    }    return ans;}int euler(int n){    int ans=1,i;    for(i=2;i*i<=n;i++)    {        if(n%i==0)        {            n/=i;            ans*=i-1;            while(n%i==0)            {                n/=i;                ans*=i;            }        }    }    if(n>1)        ans*=n-1;    return ans;}__int64 Quick_Mod(__int64 a, __int64 b, __int64 m){    __int64 res = 1,term = a % m;    while(b)    {        if(b & 1) res = (res * term) % m;        term = (term * term) % m;        b >>= 1;    }    return res%m;}int main(){    int t,n,y,x,s;    __int64 z;    matrix ans,k;    scanf("%d",&t);    while(t--)    {        k.a[0][0]=0;k.a[0][1]=1;        k.a[1][0]=1;k.a[1][1]=2;        ans.a[0][0]=0;ans.a[0][1]=0;        ans.a[1][0]=1;ans.a[1][1]=0;        scanf("%d%d%d%d",&n,&y,&x,&s);        mod=2*euler(s+1);        z=n;z*=y;        ans=mul(pow(k,z),ans);        z=(ans.a[0][0]*ans.a[1][0])%mod+mod;        z=(z%mod+mod)/2;        printf("%I64d\n",Quick_Mod(x,z,s+1));    }    return 0;}