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
^2,f[n+1]^2,f
*f[n+1],g
）

↑这是左矩阵

[0,1,0,0]

[1,4,2,1]

[0,4,1,0]

[0,0,0,1] ←这是右矩阵

是这么解释的



相信你已经看懂了

/[b]****************[/b]/

就是思维太局限了 首先欧拉降幂不知道

再后来矩阵不会构造 导致这场的GG。。

/[b]****************[/b]/

附本题代码

————————————–.

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