UVALive 4387 Tower

Source: Regionals 2008 :: Europe - Central
Problem: 已知a1,a2，且an=2a2an−1−an−2,sn=sn−1+a2n，求sn
Idea:

x=2a2,y=−1

⎡⎣⎢⎢⎢⎢a2na2n−1anan−1sn⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢x21xx2y200y22xy0y2xy0001⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a2n−1a2n−2an−1an−2sn−1⎤⎦⎥⎥⎥⎥

Code:

#include<bits/stdc++.h>using namespace std;#define fi first#define se second#define pb push_back#define lson o<<1#define rson o<<1|1#define CLR(A, X) memset(A, X, sizeof(A))#define bitcount(X) __builtin_popcountll(X)typedef long long LL;typedef pair<int, int> PII;typedef pair<string, string> PSS;const double eps = 1e-10;const double PI = acos(-1.0);//const LL MOD = 1e9+7;const auto INF = 0x3f3f3f3f;int dcmp(double x) { if(fabs(x) < eps) return 0; return x<0?-1:1; }const int MAXN = 2e3+5;int M;struct Matrix {    int c[4][4];    Matrix() { CLR(c, 0); }    void init() {        for(int i = 0; i < 4; i++) {            c[i][i] = 1;        }    }};Matrix mul(Matrix A, Matrix B) {    Matrix C;    for(int k = 0; k < 4; k++) {        for(int i = 0; i < 4; i++) if(A.c[i][k]) {            for(int j = 0; j < 4; j++) if(B.c[k][j]) {                C.c[i][j] = (C.c[i][j]+1LL*A.c[i][k]*B.c[k][j])%M;            }        }    }    return C;}Matrix qpow(Matrix a, int b) {    Matrix res;    res.init();    while(b) {        if(b & 1) res = mul(res, a);        a = mul(a, a);        b >>= 1;    }    return res;}int main() {    int X;    scanf("%d", &X);    while(X--) {        int y, n;        scanf("%d%d%d", &y, &n, &M);        int x = y<<1;        Matrix A, B, C;        A.c[0][0] = 1LL*x*x%M; A.c[0][1] = 1%M; A.c[0][2] = (M-2LL*x%M)%M;        A.c[1][0] = 1%M; A.c[2][0] = x%M; A.c[2][2] = (M-1)%M; A.c[3][3] = 1%M;        A.c[3][0] = A.c[0][0]; A.c[3][1] = 1%M; A.c[3][2] = A.c[0][2];        B.c[0][0] = 1LL*y*y%M; B.c[1][0] = 1%M; B.c[2][0] = y%M;        B.c[3][0] = (1LL*y*y%M+1)%M;        C = mul(qpow(A, n-2), B);        printf("%d\n", C.c[3][0]);    }    return 0;}