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+,M=1e6+,inf=1e9+,mod=1e9+;
const ll INF=1e18+;
ll n,x,y,s;
ll m;
struct is
{
    ll a[][];
};
is juzhenmul(is a,is b,ll hang ,ll lie,ll mod)
{
    int i,t,j;
    is ans;
    memset(ans.a,,sizeof(ans.a));
    for(i=;i<=hang;i++)
    for(t=;t<=lie;t++)
    for(j=;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&)  ans=juzhenmul(ans,a,,,mod);
        a=juzhenmul(a,a,,,mod);
        x>>=;
    }
    return ans;
}
void extend_Euclid(ll a, ll b, ll &x, ll &y)
{
    if(b == )
    {
        x = ;
        y = ;
        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=;i*i<=n;i++)
    {
        if(n%i==)
        {
            rea=rea-rea/i;
            while(n%i==)  n/=i;
        }
    }
    if(n>)
        rea=rea-rea/n;
    return rea;
}
ll Pow(ll a,ll n,ll mod)
{
    ll ans=;
    while(n)
    {
        if(n&)
        {
            ans=ans*a%mod;
        }
        a=a*a%mod;
        n>>=;
    }
    if(ans==) ans+=mod;
    return ans;
}
ll getans(ll x,ll mod)
{
    if(x==)
    return ;
    if(x==)
    return ;
    is ans,base;
    memset(ans.a,,sizeof(ans.a));
    ans.a[][]=;
    ans.a[][]=;
    base.a[][]=;
    base.a[][]=;
    base.a[][]=;
    base.a[][]=;
    ans=quickpow(ans,base,x-,mod);
    return (ans.a[][]+ans.a[][]*)%mod;
}
ll getans2(ll x,ll mod)
{
    if(x==)
    return ;
    if(x==)
    return ;
    is ans,base;
    memset(ans.a,,sizeof(ans.a));
    ans.a[][]=;
    ans.a[][]=;
    base.a[][]=;
    base.a[][]=-;
    base.a[][]=;
    base.a[][]=;
    ans=quickpow(ans,base,x-,mod);
    return ((((ans.a[][]*-ans.a[][])%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+);
        ll k;
        if(zhi%==)
        k=(getans(zhi+,m)%m)*(getans2(zhi/,m)%m)%m;
        else
        k=(getans(zhi,m)%m)*(getans2((zhi+)/,m)%m)%m;
        ll out=Pow(x,k,s+);
        printf("%lld\n",out);
    }
    return ;
}