hdu5895Mathematician QSC+矩阵快速幂+除法取余

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)=∑^n_{i=0}f(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

$f_n=2f_{n-1}+f_{n-2}$
$f_n*f_{n-1}=2f_{n-1}^2+f_{n-2}*f_{n-1}$
$2f_{n-1}^2=f_n*f_{n-1}-f_{n-2}*f_{n-1}$
/*****************************/
$2f_{n}^2=f_{n}*f_{n+1}-f_{n}*f_{n-1}$
$2f_{n-1}^2=f_{n-1}*f_{n}-f_{n-1}*f_{n-2}$
….
$2f_{2}^2=f_{2}*f_{3}-f_{2}*f_{1}$
$2f_{1}^2=f_{1}*f_{2}-f_{1}*f_{0}$
累加
$\sum_{i=1}^{n} 2f_{i}^2=f_{n}*f_{n+1}-f_{1}*f_{0}$
因为$f_{0}=0$
所以$\sum_{i=0}^{n} f_{i}^2=f_{n}*f_{n+1}/2$
所以直接构造fn的矩阵快速幂就可以求得了。。。
然后就是两个小结论要用到。。（orz，窝弱只知道其中一个啊，全场wa到死。。）
1.高次幂取膜
$a^b\%p=a^{φ(p)+b\%φ(p)}　　　b>=φ(p)$
φ(p)为p的欧拉函数
2.除法取膜
$\frac{a}{b}\%p=\frac{a\%pb}{b}$
所以小范围爆一下，大点>φ(p)就用矩阵优化。。。

#include<bits/stdc++.h>
using namespace std;

typedef long long LL;


//#define MOD 10000007
LL MOD;
struct Mat{
    int n,m;
    LL mat[9][9];
};
Mat operator *(Mat a,Mat b){
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    c.n = a.n,c.m = b.m;

    for(int i=1;i<=a.n;i++){
        for(int j=1;j<=b.m;j++){
            for(int k=1;k<=a.m;k++){
                c.mat[i][j] += (a.mat[i][k]*b.mat[k][j])%MOD;
                c.mat[i][j] %= MOD;
            }
        }
    }
    return c;
}
Mat operator +(Mat a,Mat b){
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    c.n = a.n,c.m = a.m;

    for(int i=1;i<=a.n;i++){
        for(int j=1;j<=a.m;j++){
            c.mat[i][j] = (a.mat[i][j]+b.mat[i][j])%MOD;
        }
    }
    return c;
}
Mat operator ^(Mat a,LL k){
    Mat c;
    memset(c.mat,0,sizeof(c.mat));
    c.n = a.n,c.m = a.n;
    for(int i=1;i<=a.n;i++)c.mat[i][i] = 1;

    while(k){
        if(k&1){
            c = c*a;
        }
        a = a*a;
        k>>=1;
    }
    return c;
}
void out(Mat a){
    for(int i=1;i<=a.n;i++){
        for(int j=1;j<=a.m;j++){
            printf(j==a.m? "%I64d\n":"%I64d ",a.mat[i][j]);
        }
    }
}
LL quick_pow(LL a,LL b,LL p){
    LL ret=1;
    a%=p;
    while(b){
        if(b&1) ret=(ret*a)%p;
        a=(a*a)%p;
        b/=2;
    }
    return ret;
}

LL euler_phi(LL n)
{
    LL m = (LL) sqrt(n+0.5);
    LL ans = n;
    for(LL i = 2; i <= m;i++)
        if(n%i == 0)
        {
            ans = ans/i *(i-1);
            while(n%i == 0) n/=i;
        }
    if(n>1) ans = ans/n * (n-1);
    return ans;
}
LL SUM2(int n){
    int a[13]={0,1,2,5,12,29,70,169,408,985,2378,5741};
    LL ans=0;
    for(int i=0;i<=n;i++) ans+=1LL*a[i]*a[i];
    return ans;
}

LL SUM(LL n){
    Mat pp;
    pp.n=pp.m=2;
    pp.mat[1][1]=2%MOD;
    pp.mat[1][2]=1%MOD;
    pp.mat[2][1]=1%MOD;
    pp.mat[2][2]=0%MOD;

    Mat A0;
    A0.n=1,A0.m=2;
    A0.mat[1][1]=1%MOD;
    A0.mat[1][2]=0%MOD;

    if(n==0) return 0;
    else if(n==1) return 1;

    Mat ans=pp^(n-1);
    Mat ans2=pp^(n);
    ans=A0*ans;
    ans2=A0*ans2;

    return ((ans.mat[1][1]*ans2.mat[1][1])%MOD)/2;
}

LL read()
{
    LL x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
int main()
{
    int T_T;
    scanf("%d",&T_T);
    LL n,y,x,s;
    while(T_T--){
        n = read();
        y = read();
        x = read();
        s = read();

        if(n*y<=11){
            cout<<quick_pow(x,SUM2(n*y),s+1)<<endl;
        }
        else{
            MOD=euler_phi(s+1);
            MOD=MOD*2;
            LL temp=SUM(n*y);
            cout<<quick_pow(x,temp+MOD,s+1)<<endl;
        }
    }
    return 0;
}