Mathematician QSC（欧拉定理+矩阵快速幂）

Mathematician QSC

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

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  xg(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

 

Recommend

wange2014

首先这题的关键是要降幂，由欧拉定理（，φ(n)为小于等于n的质数的个数，具体怎么求参考网上），可以写出原题

所要求的可以转化为x​g(n×y) mod φ(s+1)+φ(s+1)​​ mod (s+1).这样就可以用快速幂在解决。

接下来就是求g(n)的过程，一开始的反应是矩阵快速幂，所以就要构造矩阵，这里写下怎么构造矩阵的思路。

我们要求g(n)，g(n+1)=g(n)+f​2​​(n+1)，所以右矩阵中肯定有g(n)和f​2​​(n+1)。接下来，因为有f​2​​(n+1)，那么最后的矩阵会有

f​2​​(n+2)，我们把f​2​​(n+2)展开，f^2(n + 2) = (f(n) + 2f(n + 1))^2 = f^2(n) + 4f(n)f(n + 1) + 4f^2(n + 1)f​2​​(n+2)=(f(n)+2f(n+1))​2​​=f​2​​(n)+4f(n)f(n+1)+4f​2​​(n+1)，所以右边的矩阵还有

f​2​​(n)和f(n)f(n+1)，这样我们就可以写出矩阵了

接下来的是就是用模板了，刚好借这个机会存个模板。

f^2(n + 2) = (f(n) + 2f(n + 1))^2 = f^2(n) + 4f(n)f(n + 1) + 4f^2(n + 1的

#include<iostream>
#include<cstdio>
#include<string.h>
#include<algorithm>
using namespace std;
long long SMod;
bool valid[1000005];
int prime[1000005];
int Max=1000005;
struct Matrix
{
    long long m[5][5];
};
Matrix Mul(Matrix a,Matrix b)
{
    Matrix c;
    memset(c.m,0,sizeof(c.m));
    for(int i=0;i<4;i++)
        for(int j=0;j<4;j++)
            for(int k=0;k<4;k++)
                c.m[j][k] += ((a.m[j][i]*b.m[i][k])%SMod + SMod)%SMod;
    return c;
}
Matrix fastm(Matrix a,long long n) 
{
    Matrix res;
    memset(res.m,0,sizeof(res.m));
    res.m[0][0] = res.m[1][1] = res.m[2][2] =res.m[3][3]= 1;
    while(n)
    {
        if(n&1)
        {
            res = Mul(res,a);
            
        }
        n>>=1;
        a = Mul(a,a);
        
        //cout<<"yes"<<endl;
    }
    return res;
}

long long euler(long long n)  
{  
    long long ans=1;  
    int 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;  
}  
long long PowerMod(long long a, long long b, long long c)
{
    long long ans = 1;
    a = a % c;
    while(b>0) {
        if(b % 2 == 1)
        ans = (ans * a) % c;
        b = b/2;
        a = (a * a) % c;
    }
    return ans;
}
int main()
{
    int cas;
    scanf("%d",&cas);
    while(cas--)
    {
        long long n,y,x,s;
        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);
        SMod=euler(s+1);
        //cout<<SMod<<endl;
        Matrix mm;
        mm.m[0][0]=0;mm.m[0][1]=1;mm.m[0][2]=0;mm.m[0][3]=0;
        mm.m[1][0]=1;mm.m[1][1]=4;mm.m[1][2]=4;mm.m[1][3]=0;
        mm.m[2][0]=0;mm.m[2][1]=2;mm.m[2][2]=1;mm.m[2][3]=0;
        mm.m[3][0]=0;mm.m[3][1]=1;mm.m[3][2]=0;mm.m[3][3]=1;
        Matrix tem=fastm(mm,n*y);
        //cout<<tem.m[3][1]<<endl;
        long long ss=(tem.m[3][1]%SMod+SMod);
        long long aa=PowerMod(x,ss,s+1);
        //cout<<PowerMod(3,2,555)<<endl;
        cout<<aa<<endl;
    }
}