HDU 5895&&2016 ACM/ICPC Asia Regional Shenyang Online1004 Mathematician QSC [矩阵加速+欧拉降幂]【数论】

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5895
————————.
Mathematician QSC

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 131072/131072 K (Java/Others)
Total Submission(s): 222 Accepted Submission(s): 109

Problem Description
QSC dream of becoming a mathematician, he believes that everything in this world has a mathematical law.

Through unremitting efforts, one day he finally found the QSC sequence, it is a very magical sequence, can be calculated by a series of calculations to predict the results of a course of a semester of a student.

This sequence is such like that, first of all,f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2)Then the definition of the QSC sequence is g(n)=∑ni=0f(i)^2. If we know the birthday of the student is n, the year at the beginning of the semester is y, the course number x and the course total score s, then the forecast mark is x^g(n∗y)%(s+1).
QSC sequence published caused a sensation, after a number of students to find out the results of the prediction is very accurate, the shortcoming is the complex calculation. As clever as you are, can you write a program to predict the mark?

Input
First line is an integer T(1≤T≤1000).

The next T lines were given n, y, x, s, respectively.

n、x is 8 bits decimal integer, for example, 00001234.

y is 4 bits decimal integer, for example, 1234.
n、x、y are not negetive.

1≤s≤100000000

Output
For each test case the output is only one integer number ans in a line.

Sample Input
2
20160830 2016 12345678 666
20101010 2014 03030303 333

Sample Output
1
317

Source
2016 ACM/ICPC Asia Regional Shenyang Online

————————.

题目大意：
就是给你四个数n,y,x,s,
让你求x^g(n∗y)%(s+1).
其中g(n)=∑(i->n)f(i)^2;
f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2)

解题思路：
对于x的指数g(n) 是一个很大的数 所以需要想办法把它改成我们能计算的 就是欧拉降幂

然后只要求解g(n)就行了
很容易的想到g(n)=f(n)*f(n+1)/2 其中f(n)只用一个矩阵加速就能很快地求解
但是随后发现这样并不行 因为f(n)已经是对Phi(s+1)取模之后的数了 在/2之后之就会出错 然后想到求逆元的办法解决 但是随后发现虽然2是一个质数但是并不能满足gcd(2,s+1)==1 因为s+1%2可能等于0 于是这个思路就GG了。。。

最后还是看了ICPCCamp的题解 才知道解法
最终还是一个矩阵快速幂
（f[n]^2,f[n+1]^2,f[n]*f[n+1],g[n]）
↑这是左矩阵
[0,1,0,0]
[1,4,2,1]
[0,4,1,0]
[0,0,0,1] ←这是右矩阵

是这么解释的

相信你已经看懂了

/****************/
就是思维太局限了 首先欧拉降幂不知道
再后来矩阵不会构造 导致这场的GG。。
/****************/

附本题代码
————————————–.

#include<bits/stdc++.h>using namespace std;typedef long long LL;const int maxn = 505;const double  Pi =  acos(-1);#define pb push_back#define lalal puts("****");const int M = 4;int MOD ;struct Matrix{    LL m[M][M];    void clearO()    {        for(int i=0; i<M; i++) //初始化矩阵            for(int j=0; j<M; j++)                m[i][j]= 0;    }    void clearE()    {        for(int i=0; i<M; i++) //初始化矩阵            for(int j=0; j<M; j++)                m[i][j]= (i==j);    }    void display()    {        for(int i=0; i<M; i++)            {                for(int j=0; j<M; j++)                printf("%d ",m[i][j]);                puts("");            }    }};Matrix operator * (Matrix a,Matrix b){    Matrix c;    c.clearO();    for(int k=0; k<M; k++)        for(int i=0; i<M; i++) //实现矩阵乘法        {            if(a.m[i][k] <= 0)  continue;            for(int j=0; j<M; j++)            {                if(b.m[k][j] <= 0)    continue;                c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]+MOD)%MOD;            }        }    return c;}Matrix operator ^ (Matrix a,LL b){    Matrix c;    c.clearE();    while(b)    {        if(b&1) c= c * a ;        b >>= 1;        a = a * a ;    }    return c;}int Is_or[101001];int prime[13000],kpri;void Prime(){    int n=100001;    kpri=0;    memset(Is_or,1,sizeof(Is_or));    Is_or[0]=Is_or[1]=0;    for(int i=2; i<n; i++)    {        if(Is_or[i])        {            prime[kpri++]=i;            for(int j=i+i; j<n; j+=i)            {                Is_or[j]=0;            }        }    }    //prime[kpri]=10007;    //printf("%d\n",kpri);    return ;}LL Phi(LL n){    LL rea=n;    for(int i=0; prime[i]*prime[i]<=n; i++)    {        if(n%prime[i]==0)        {            rea=rea-rea/prime[i];            while(n%prime[i]==0)                n/=prime[i];        }    }    if(n>1)  rea=rea-rea/n;    return rea;}LL qmod(LL a,LL b){    LL res= 1;    while(b)    {        if(b&1) res=(res*a)%MOD;        b>>=1;        a=(a*a)%MOD;    }    return res;}int main(){    Prime();    int _;    while(~scanf("%d",&_))    {        while(_--)        {            LL n,y,x,s;            scanf("%I64d%I64d%I64d%I64d",&n,&y,&x,&s);            MOD = Phi(s+1);            Matrix a,b;            a.clearO(),b.clearO();            a.m[0][1]=1;            b.m[0][1]=1;            b.m[1][0]=1,b.m[1][1]=4,b.m[1][2]=2,b.m[1][3]=1;            b.m[2][1]=4,b.m[2][2]=1;            b.m[3][3]=1;            b=b^(n*y);            a=a*b;            LL zhi = a.m[0][3]%MOD+MOD;            MOD=s+1;            printf("%I64d\n",qmod(x%MOD,zhi));        }    }    return 0;}