uvalive 4387 Tower - 矩阵

题目大意

给了一个数列的递推公式，a[n]=2*a[2]*a[n-1]-a[n-2]，现在让你求数列前n项平方和模M。

数据范围：1<=a[2],m<=10^9, M<=10^9, 2<=n<=10^9

思路

给出数列递推式让你求第n项或者前n项的关系时，一般都是通过构造出一个矩阵，得出第n项的矩阵递推式，然后通过矩阵快速幂来解决，这题也一样，只不过递推式稍微难想一点。

构造一个1*4的矩阵
(a[n]*a[n], a[n-1]*a[n-1], a[n-1]*a[n-2], s[n-1]), 其中s[n]代表前n项平方和
s[n] = s[n-1] + a[n]^2
那么通过递推式可以得出这个矩阵的递推式，再使用快速幂就可求出答案。

/******************************************************
推公式。。
矩阵快速幂
S[n] 表示前n项平方和
a[n] = 2*a2*a[n-1] - a[n-2]
p = 2*a2;
初始状态 (a[3]^2,a[2]^2,a[2]*a[1],S[2])

(a[n+1]^2,a[n]^2,a[n]*a[n-1],S[n]) = (a[n]^2,a[n-1]^2,a[n-1]*a[n-2],S[n-1])*
    ( p^2,     1,  0,  1)
    ( 1-2*p^2, 0,  p,  0)
    ( 2*p,     0,  -1, 0)
    ( 0,       0,  0,  1)
***********************************************/
#include <cstdio>
#include <iostream>
#include <cstring>
using namespace std;
typedef long long LL;

LL a2,a3,n,mod,p;

/******************
矩阵快速幂
*******************/
const int maxn = 10;
struct Matrix{
    int n,m;
    LL a[maxn][maxn];

    Matrix(){
        memset(a,0,sizeof(a));
    }

    void InitOne(int x,int y){
        n =x, m = y;
        for(int i=0;i<n;i++){
            for(int j=0;j<m;j++) if(i==j) a[i][j] = 1;
        }
    }

    void debug(){
        printf("n:%d m:%d\n",n,m);
        for(int i=0;i<n;i++){
            for(int j=0;j<m;j++){
                cout << a[i][j] <<" ";
            }
            cout << endl;
        }
    }
};

Matrix multi(Matrix a,Matrix b){
    Matrix ret;
    for(int i=0;i<4;i++){
        for(int j=0;j<4;j++){

            for(int k=0;k<4;k++){
                ret.a[i][j] += a.a[i][k]*b.a[k][j];
            }
            ret.a[i][j] %= mod;
        }
    }
    return ret;
}

Matrix pow_mod(Matrix x,LL n){
    Matrix res;
    res.InitOne(4,4);
    while(n){
        if(n&1) res = multi(res,x);
        n >>= 1;
        x = multi(x,x);
    }
    return res;
}
///

LL solve(){
    if(a2 == 1 || mod == 1) return n%mod;
    if(n == 2) return (1 + a2*a2)%mod;

    a3 = (2*a2*a2%mod - 1 + mod)%mod;
    p = 2*a2%mod;

    LL x[5][5] = {
        {p*p%mod,     1,  0,  1},
        {(1-2*p%mod*p%mod + mod)%mod, 0,  p,  0},
        {2*p%mod,     0,  mod-1, 0},
        {0,       0,  0,  1}
    };

    Matrix ans,start;
    start.a[0][0] = a3*a3%mod;
    start.a[0][1] = a2*a2%mod;
    start.a[0][2] = a2%mod;
    start.a[0][3] = 1+a2*a2%mod;

    for(int i=0;i<4;i++){
        for(int j=0;j<4;j++) ans.a[i][j] = x[i][j];
    }

    ans = pow_mod(ans,n-2);
    long long ret = 0;
    for(int i=0;i<4;i++){
        ret += start.a[0][i]*ans.a[i][3]%mod;
        ret %= mod;
    }
    return ret;
}

int main(){
    //freopen("tow1.in","r",stdin);
    //freopen("tow1.out","w",stdout);
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%lld%lld%lld",&a2,&n,&mod);
        printf("%lld\n",solve());
    }
    return 0;
}