hdu 5895 Mathematician QSC 指数循环节+矩阵快速幂

Mathematician QSC

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 131072/131072 K (Java/Others)

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

思路：首先求A^B%C=A^(B%phi(C)+phi(C))%C  B>=phi(C)指数循环节;

　　   然后，求g函数，f(n)显然可以用矩阵快速幂写，g（n）=f（n）*f（n+1）/2；因为/2，模除法，首先想到逆元，然而模不一定是奇数，偶数的情况2无逆元；

　　　现在怎么处理2，f(n)与f(n+1)必然有一个是偶数，发现除2后的递推式更改为f(n)=6*f（n-1）-f(n-2);

　　　ps:一个小技巧处理2，mdzz，模的数*2，答案/2;

　　　详见代码；

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pi (4*atan(1.0))
const int N=1e5+10,M=1e6+1010,inf=1e9+10,mod=1e9+7;
const ll INF=1e18+10;
ll n,x,y,s;
ll m;
struct is
{
    ll a[10][10];
};
is juzhenmul(is a,is b,ll hang ,ll lie,ll mod)
{
    int i,t,j;
    is ans;
    memset(ans.a,0,sizeof(ans.a));
    for(i=1;i<=hang;i++)
    for(t=1;t<=lie;t++)
    for(j=1;j<=lie;j++)
    {
        ans.a[i][t]+=(a.a[i][j]*b.a[j][t]);
        ans.a[i][t]%=mod;
    }
    return ans;
}
is quickpow(is ans,is a,ll x,ll mod)
{
    while(x)
    {
        if(x&1)  ans=juzhenmul(ans,a,2,2,mod);
        a=juzhenmul(a,a,2,2,mod);
        x>>=1;
    }
    return ans;
}
void extend_Euclid(ll a, ll b, ll &x, ll &y)
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return;
    }
    extend_Euclid(b, a % b, x, y);
    ll tmp = x;
    x = y;
    y = tmp - (a / b) * y;
}
ll phi(ll n)
{
    ll i,rea=n;
    for(i=2;i*i<=n;i++)
    {
        if(n%i==0)
        {
            rea=rea-rea/i;
            while(n%i==0)  n/=i;
        }
    }
    if(n>1)
        rea=rea-rea/n;
    return rea;
}
ll Pow(ll a,ll n,ll mod)
{
    ll ans=1;
    while(n)
    {
        if(n&1)
        {
            ans=ans*a%mod;
        }
        a=a*a%mod;
        n>>=1;
    }
    if(ans==0) ans+=mod;
    return ans;
}
ll getans(ll x,ll mod)
{
    if(x==0)
    return 0;
    if(x==1)
    return 1;
    is ans,base;
    memset(ans.a,0,sizeof(ans.a));
    ans.a[1][1]=1;
    ans.a[2][2]=1;
    base.a[1][1]=0;
    base.a[1][2]=1;
    base.a[2][1]=1;
    base.a[2][2]=2;
    ans=quickpow(ans,base,x-2,mod);
    return (ans.a[2][1]+ans.a[2][2]*2)%mod;
}
ll getans2(ll x,ll mod)
{
    if(x==0)
    return 0;
    if(x==1)
    return 1;
    is ans,base;
    memset(ans.a,0,sizeof(ans.a));
    ans.a[1][1]=1;
    ans.a[2][2]=1;
    base.a[1][1]=6;
    base.a[1][2]=-1;
    base.a[2][1]=1;
    base.a[2][2]=0;
    ans=quickpow(ans,base,x-2,mod);
    return ((((ans.a[1][1]*6-ans.a[2][1])%mod)+mod)%mod);
}
int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);
        ll zhi=n*y;
        m=phi(s+1);
        ll k;
        if(zhi%2==0)
        k=(getans(zhi+1,m)%m)*(getans2(zhi/2,m)%m)%m;
        else
        k=(getans(zhi,m)%m)*(getans2((zhi+1)/2,m)%m)%m;
        ll out=Pow(x,k,s+1);
        printf("%lld
",out);
    }
    return 0;
}