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