HDU5895 Mathematician QSC(构造矩阵+矩阵快速幂+幂次循环节)

Mathematician QSC

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


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

题面的意思就是
求x^g(n*y)%(s+1),其实就是A^BmodC的问题。
这里的g(n*y)怎么求呢，其实一眼看上去就知道要构造矩阵然后进行矩阵快速幂。
由于f(n)=2*f(n-1)+f(n-2),这里有g(n) = ∑f(i)^2.
那么我们可以平方展开得到 f(n)^2=4*f(n-1)^2+4*f(n-1)*f(n-2)+f(n-2)^2;
这里我们就可以来构造矩阵啦：
我下面的代码是六维的，其实只要四维就好了，f(n-1)和f(n-2)没有也没影响。
我就写个四维的吧。
| 4 4 1 0 |  * | f(n-1)^2   |  =  | f(n)^2   |
| 2 1 0 0 |    |f(n-1)f(n-2)|     |f(n)f(n-1)|
| 1 0 0 0 |    | f(n-2)^2   |     | f(n-1)^2 |
| 4 4 1 1 |    |    g(n)    |     |  g(n)    |
f(n-1)f(n-2)->f(n)(n-1)其实就是f(n-1)f(n-2)+2*f(n-1)^2
这样矩阵就构造出来了。
那问题来了，这里A^BmodC的B的值很大，那我们就需要降幂。
这里需要用到指数循环节的知识。
A^BmodC = A^(B%phi(C)+phi(C)) (mod C) (B>=phi(C))
(phi(C)这里是欧拉函数)。
矩阵快速幂和指数循环节结合在一起就可以求g(n*y)了。
接下来就用快速幂随便搞一下就好了。 

#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cctype>
#include<cmath>
#include<ctime>
#include<string>
#include<stack>
#include<deque>
#include<queue>
#include<list>
#include<set>
#include<map>
#include<cstdio>
#include<limits.h>
#define fir first
#define sec second
#define fin freopen("/home/ostreambaba/文档/input.txt", "r", stdin)
#define fout freopen("/home/ostreambaba/文档/output.txt", "w", stdout)
#define mes(x, m) memset(x, m, sizeof(x))
#define pii pair<int, int>
#define Pll pair<ll, ll>
#define INF 1e9+7
#define Pi 4.0*atan(1.0)
#define MOD 1000000007

#define lowbit(x) (x&(-x))
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
#define ls rt<<1
#define rs rt<<1|1


typedef long long ll;
typedef unsigned long long ull;
const double eps = 1e-12;
const int maxn = 6;
using namespace std;
//#define TIME

inline int read(){
    int 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;
}
ull mod;
int s;
struct matrix{
    ull mat[maxn][maxn];
    void init(){
        mes(mat, 0);
        for(int i = 0; i < 6; ++i){
            mat[i][i] = 1;
        }
    }
    void clear(){
        mes(mat, 0);
    }
    void output(){
        for(int i = 0; i < 6; ++i){
            for(int j = 0; j < 6; ++j){
                //  printf("%lld ", mat[i][j]);
                cout << mat[i][j] << " ";
            }
            printf("\n");
        }
    }
    matrix operator *(const matrix &base){
        matrix tmp;
        tmp.clear();
        for(int i = 0; i < 6; ++i){
            for(int j = 0; j < 6; ++j){
                for(int k = 0; k < 6; ++k){
                    tmp.mat[i][j] = (tmp.mat[i][j] + mat[i][k]*base.mat[k][j]);
                    if(tmp.mat[i][j] > mod){ //指数循环节运用
                        tmp.mat[i][j] %= mod;
                    }
                }
            }
        }
        return tmp;
    }
};
int eular(int n)  //欧拉函数
{
    int res = 1;
    if(1 == n){
        return res;
    }
    int end = sqrt(n+1);
    for(int i = 2; i <= end; ++i){
        if(n%i == 0){
            n/=i;
            res *= (i-1);
            while(n%i == 0){
                n/=i;
                res*=i;
            }
        }
    }
    if(n > 1){
        res *= (n - 1);
    }
    return res;
}

matrix matrix_fast_mod(ull m, matrix base) //矩阵快速幂
{
    matrix res;
    //  res.output();
    res.init();
    while(m){
        if(m&1){
            res = res*base;
        }
        base = base*base;
        m >>= 1;
    }
    return res;
}

ull fast_mod(ull m, ull x){  //快速幂
    ull res = 1;
    while(m){
        if(m&1){
            res = (res*x)%(s+1);
        }
        x=(x*x)%(s+1);
        m >>= 1;
    }
    return res;
}


int main()
{
    fin;
    int Case, n, x, y;
    ull ret, cx;
    Case = read();
    while(Case--){
        scanf("%d%d%d%d", &n, &y, &x, &s);
        matrix base = {  //构造矩阵
            4, 4, 1, 0, 0, 0,
            2, 1, 0, 0, 0, 0,
            1, 0, 0, 0, 0, 0,
            0, 0, 0, 2, 1, 0,  //
            0, 0, 0, 1, 0, 0,   //  这一行和上面一行去掉也行，f(n-1), f(n-2)对结果无影响
            4, 4, 1, 0, 0, 1
        };
        matrix p;
        p.clear();
        p.mat[0][0] = 1;
        p.mat[1][0] = 0;
        p.mat[2][0] = 0;
        p.mat[3][0] = 1;
        p.mat[4][0] = 0;
        p.mat[5][0] = 1;
        cx = (ull)n*y;
        mod = eular(s+1);
        if(1 == cx){
            ret = 1;
        }
        else{
            base = matrix_fast_mod(cx-1, base);
            p = base*p;
            ret = p.mat[5][0];
        }
        ret =  fast_mod(ret, (ull)x)%(s+1);
        cout << ret << endl;
    }
    return 0;
}