博客

$$ \text{对于矩阵}M,\text{构造其特征多项式}G\left( x \right) \\ \text{则有}M^n=G\left( M \right) Q\left( M \right) +R\left( M \right) \\ \text{即}M^n\equiv R\left( M \right) \left( mod\,\,G\left( M \right) \right) \\ \text{即}x^n\equiv R\left( x \right) \left( \text{mod }G\left( x \right) \right) \\ \text{考虑特征多项式如何计算。} \\ G\left( x \right) =\det \left( M-x\text{I} \right) ,\text{其中}I\text{为单位矩阵} \\ \text{由于}G\left( x \right) \text{为}n\text{次多项式}\left( n\text{为矩阵大小} \right) ,\text{故可以通过插值求出来对应的特征多项式。} \\
\\ \text{然后多项式快速幂}+\text{取模即可} $$

#include <algorithm>
#include <cstring>
#include <iostream>
#include <string>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;

const int_t mod = 1000000007;
const int_t LARGE = 103;
struct Matrix {
    int_t data[LARGE + 1][LARGE + 1];
    int_t n;
    Matrix(int_t n) {
        this->n = n;
        memset(data, 0, sizeof(data));
    }
    int_t at(int_t r, int_t c) const { return data[r][c]; }
    int_t* operator[](int_t r) { return data[r]; }
    Matrix operator*(const Matrix& rhs) const {
        Matrix result(n);
        for (int_t i = 1; i <= n; i++)
            for (int_t j = 1; j <= n; j++)
                for (int_t k = 1; k <= n; k++) {
                    result[i][j] =
                        (result[i][j] + data[i][k] * rhs.at(k, j) % mod) % mod;
                }
        return result;
    }
    Matrix operator+(const Matrix& rhs) const {
        Matrix result(n);
        for (int_t i = 1; i <= n; i++)
            for (int_t j = 1; j <= n; j++) {
                result[i][j] = (data[i][j] + rhs.at(i, j)) % mod;
            }
        return result;
    }
    Matrix operator-() const {
        Matrix result = *this;
        for (int_t i = 1; i <= n; i++)
            for (int_t j = 1; j <= n; j++)
                result[i][j] = (mod - result[i][j] + mod) % mod;
        return result;
    }
    Matrix operator*(int_t k) const {
        Matrix result = *this;
        for (int_t i = 1; i <= n; i++)
            for (int_t j = 1; j <= n; j++)
                result[i][j] = (result[i][j] * k) % mod;
        return result;
    }
    Matrix operator-(const Matrix& rhs) const { return (*this) + (-rhs); }
};
int_t power(int_t base, int_t index) {
    if (index < 0) {
        index = (index % (mod - 1) + mod - 1) % (mod - 1);
    }
    int_t result = 1;
    while (index) {
        if (index & 1) result = result * base % mod;
        base = base * base % mod;
        index >>= 1;
    }
    return result;
}

void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
    static int_t Cx[LARGE + 1];
    std::fill(Cx, Cx + n + m + 1, 0);
    for (int_t i = 0; i <= n; i++)
        for (int_t j = 0; j <= m; j++)
            Cx[i + j] = (Cx[i + j] + A[i] * B[j] % mod) % mod;
    std::copy(Cx, Cx + n + m + 1, C);
}
void poly_div(const int_t* A, int_t n, const int_t* B, int_t m, int_t* R,
              int_t* Q) {
    while (B[m] % mod == 0) m--;
    if (n < m) {
        std::fill(Q, Q + n - m + 1, 0);
        std::copy(A, A + n + 1, R);
        return;
    }
    std::copy(A, A + n + 1, R);
    for (int_t i = n - m; i >= 0; i--) {
        Q[i] = R[i + m] * power(B[m], -1) % mod;
        for (int_t j = m; j >= 0; j--) {
            R[i + j] = (R[i + j] - Q[i] * B[j] % mod + mod) % mod;
        }
    }
}
void interpolate(int_t* X, int_t* Y, int_t n, int_t* result) {
    static int_t prod[LARGE + 1];
    prod[0] = 1;
    for (int_t i = 1; i <= n; i++) {
        int_t Z[] = {mod - X[i], 1};
        poly_mul(prod, i - 1, Z, 1, prod);
    }
#ifdef DEBUG
    for (int_t i = 0; i <= n; i++) cout << prod[i] << " ";
    cout << endl;

#endif
    for (int_t i = 1; i <= n; i++) {
        int_t coef = 1;
        for (int_t j = 1; j <= n; j++)
            if (i != j) {
                coef = coef * (X[i] - X[j] + mod) % mod;
            }
        coef = Y[i] * power(coef, -1) % mod;
        static int_t curr[LARGE + 1], R[LARGE + 1];
        int_t B[] = {mod - X[i], 1};
        poly_div(prod, n, B, 1, R, curr);
#ifdef DEBUG
        cout << "div x-" << X[i] << " = ";
        for (int_t i = 0; i <= n - 1; i++) cout << curr[i] << " ";
        cout << endl;
        cout << "coef = " << coef << endl;
        cout << endl;
#endif
        for (int_t i = 0; i < n; i++)
            result[i] = (result[i] + coef * curr[i] % mod) % mod;
    }
}
int_t det(Matrix mat) {
    int_t prod = 1;
    for (int_t i = 1; i <= mat.n; i++) {
        int_t next = -1;
        for (int_t j = i; j <= mat.n; j++)
            if (mat[i][j]) {
                next = j;
                break;
            }
        if (next == -1) return 0;
        if (next != i) {
            prod *= -1;
            std::swap(mat.data[i], mat.data[next]);
        }
        for (int_t j = i + 1; j <= mat.n; j++) {
            int_t f = mat[j][i] * power(mat[i][i], -1) % mod;
            for (int_t k = i; k <= mat.n; k++) {
                mat[j][k] = (mat[j][k] - f * mat[i][k] % mod + mod) % mod;
            }
        }
    }
    prod = mod + prod;
    for (int_t i = 1; i <= mat.n; i++) prod = prod * mat[i][i] % mod;
    return prod;
}
void poly_power(const int_t* A, int_t n, const int_t* B, int_t m,
                const std::string& str, int_t* result) {
    static int_t base[LARGE + 1], Q[LARGE + 1];
    result[0] = 1;
    std::copy(A, A + n + 1, base);
    for (int_t i = 0; i < str.length(); i++) {
        if (str[i] == '1') {
            poly_mul(result, m - 1, base, m - 1, result);
            poly_div(result, 2 * m - 1, B, m, result, Q);
        }
        poly_mul(base, m - 1, base, m - 1, base);
        poly_div(base, 2 * m - 1, B, m, base, Q);
    }
}
int main() {
    std::string index;
    cin >> index;
    std::reverse(index.begin(), index.end());

    int_t n;
    cin >> n;
    Matrix mat(n), I(n);
    static int_t X[LARGE + 1], Y[LARGE + 1], result[LARGE + 1];
    for (int_t i = 1; i <= n; i++) {
        for (int_t j = 1; j <= n; j++) cin >> mat[i][j];
        I[i][i] = 1;
    }

    for (int_t i = 1; i <= 1 + n; i++) {
        X[i] = i;
#ifdef DEBUG
        // auto curr=mat - I * i;
        // cout<<"sub "
#endif
        Y[i] = det(mat - I * i);
#ifdef DEBUG

        cout << "(" << X[i] << "," << Y[i] << ")" << endl;
#endif
    }
    interpolate(X, Y, n + 1, result);
    // for (int_t i = 0; i <= n; i++) cout << result[i] << " ";
    // cout << endl;
    static int_t A[LARGE + 1];
    static int_t result_poly[LARGE + 1];
    A[1] = 1;
    poly_power(A, 1, result, n, index, result_poly);
#ifdef DEBUG
    for (int_t i = 0; i < n; i++) cout << result_poly[i] << " ";
    cout << endl;
#endif
    Matrix result_mat(n), curr = I;
    for (int_t i = 0; i < n; i++) {
        result_mat = result_mat + curr * result_poly[i];
        curr = curr * mat;
    }
    for (int_t i = 1; i <= n; i++) {
        for (int_t j = 1; j <= n; j++) {
            cout << result_mat[i][j] << " ";
        }
        cout << endl;
    }
    return 0;
}

$$ \sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^ia_{i\text{mod}4}}=\sum_{0\le j\le 3}{a_j\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\left[ i\equiv j\left( mod\,\,4 \right) \right]}} \\ =\sum_{0\le j\le 3}{\begin{array}{c} a_j\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\left[ i-j\equiv 0\left( mod\,\,4 \right) \right]}\\ \end{array}} \\ \text{令}\omega _4\text{为模998244353意义下的四次单位根，则有}\sum_{0\le i<3}{\omega _{4}^{i}}=0 \\ \text{显然对于任意非零的}k,\text{有}\sum_{0\le i<3}{\omega _{4}^{ik}}=\text{0,而对于}k=\text{0则有}\sum_{0\le i<3}{\omega _{4}^{ik}}=4 \\ \text{故}\sum_{0\le j\le 3}{a_j\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\left[ i-j\equiv 0\left( mod\,\,4 \right) \right]}}=\sum_{0\le j\le 3}{\begin{array}{c} a_j\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\frac{1}{4}\sum_{0\le k\le 3}{\omega _{4}^{k\left( i-j \right)}}}\\ \end{array}} \\ \text{当且仅当}i-j=0\left( \text{mod}4 \right) \text{时，有}\frac{1}{4}\sum_{0\le k\le 3}{\omega _{4}^{k\left( i-j \right)}}=1 \\ \sum_{0\le j\le 3}{\begin{array}{c} a_j\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\frac{1}{4}\sum_{0\le k\le 3}{\omega _{4}^{k\left( i-j \right)}}}\\ \end{array}}=\frac{1}{4}\sum_{0\le j\le 3}{a_j\sum_{0\le k\le 3}{\omega _{4}^{-kj}\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\omega _{4}^{ki}}}} \\ =\frac{1}{4}\sum_{0\le j\le 3}{\begin{array}{c} a_j\sum_{0\le k\le 3}{\omega _{4}^{-kj}\sum_{0\le i\le n}{\left( \begin{array}{c} n\\ i\\ \end{array} \right) s^i\left( \omega _{4}^{k} \right) ^i}}\\ \end{array}} \\ =\frac{1}{4}\sum_{0\le j\le 3}{\begin{array}{c} a_j\sum_{0\le k\le 3}{\omega _{4}^{-kj}\left( s\omega _{4}^{k}+1 \right) ^n}\\ \end{array}} \\ $$

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <vector>
#include <inttypes.h>
#define debug(x) std::cout << #x << " = " << x << std::endl;

typedef long long int int_t;

using std::cin;
using std::endl;
using std::cout;

const int_t mod = 998244353;

int_t power(int_t base,int_t index) {
	int_t result=1;
	index=(index%(mod-1)+mod-1)%(mod-1);
	while(index) {
		if(index&1) result=result*base%mod;
		index>>=1;
		base=base*base%mod;
	}
	return result;
}
int main() {
	int_t T;
	cin>>T;
	const int_t g4=power(3,(mod-1)/4);
	while(T--) {
		int_t n,s;
		static int_t a[4];
		scanf("%lld%lld",&n,&s);
		for(int_t i=0; i<=3; i++)scanf("%lld",&a[i]);
		int_t result=0;
		for(int_t j=0; j<=3; j++) {
			int_t sum=0;
			for(int_t k=0; k<=3; k++) {
				sum=(sum+power(g4,-k*j)*power(1+s*power(g4,k)%mod,n)%mod)%mod;
			}
			sum=sum*a[j]%mod;
			result=(result+sum)%mod;
		}
		result=result*power(4,-1)%mod;
		printf("%lld\n",result);
	}
	return 0;
}

询问按照x从小到大分类。

然后按照x递增处理，每次把根到x的路径上深度为p的点的值加上$p^k-(p-1)^k$。

这样对于任何一条到根深度为p的路径，其权值为$p^k$。

使用树剖+线段树维护。

线段树的每个叶子节点有一个权值$p^k-(p-1)^k$，这个权值不会改变，同时需要维护区间$a_i\times(p_i^k-(p_i-1)^k)$的和。

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <vector>
#include <inttypes.h>
#define debug(x) std::cout << #x << " = " << x << std::endl;

typedef long long int int_t;

using std::cin;
using std::endl;
using std::cout;

const int_t mod = 998244353;
const int_t LARGE = 5e4;
int_t power(int_t base,int_t index){
    int_t result=1;
    while(index) {
        if(index&1) result=result*base%mod;
        index>>=1;
        base=base*base%mod;
    }
    return result;
}
struct Node{
    int_t begin,end;
    int_t value = 0;
    Node*left=nullptr,*right=nullptr;
    int_t val1=0,val2=0,mark=0;
    Node(int_t begin,int_t end):begin(begin),end(end){}
    
    void add(int_t x){
        mark+=x;
        val1+=x;
        value=(value+val2*x%mod)%mod;
    }
    void pushDown(){
        if(begin!=end){
            left->add(mark);
            right->add(mark);
            mark=0;
        }
    }
    void maintain(){
        if(begin!=end){
            value=(left->value+right->value)%mod;
            val1=left->val1+right->val1;
            val2=(left->val2+right->val2)%mod;
        }
    }
    void add(int_t begin,int_t end,int_t x){
        if(end<this->begin||begin>this->end) return;
        if(this->begin>=begin&&this->end<=end){
            this->add(x);
            return;
        }
        pushDown();
        left->add(begin,end,x);
        right->add(begin,end,x);
        maintain();
    }
    int_t query(int_t begin,int_t end){
        if(end<this->begin||begin>this->end) return 0;
        if(this->begin>=begin&&this->end<=end) return this->value;
        pushDown();
        return (left->query(begin,end)+right->query(begin,end))%mod;
    }
    static Node* build(int_t begin,int_t end,int_t* arr){
        Node* node=new Node(begin,end);
        if(begin!=end){
            int_t mid=(begin+end)/2;
            node->left=build(begin,mid,arr);
            node->right=build(mid+1,end,arr);
        }else{
            node->val2=arr[begin];
        }
        node->maintain();
        return node;
    }
};
Node* root;
int_t n,q,k;
std::vector<int_t> graph[LARGE+1];
int_t depth[LARGE+1],top[LARGE+1],maxChd[LARGE+1],size[LARGE+1];
int_t parent[LARGE+1];
int_t DFSN[LARGE+1],byDFSN[LARGE+1],target[LARGE+1];
void DFS1(int_t vtx,int_t from=0,int_t depth=1){
    parent[vtx]=from,::depth[vtx]=depth;
    size[vtx]=1;
    for(int_t to:graph[vtx]){
        DFS1(to,vtx,depth+1);
        if(size[to]>size[maxChd[vtx]]) maxChd[vtx]=to;
        size[vtx]+=size[to];
    }
}
void DFS2(int_t vtx){
    DFSN[vtx]=++DFSN[0];
    byDFSN[DFSN[vtx]]=vtx;
    target[DFSN[vtx]] = (power(depth[vtx],k) - power(depth[parent[vtx]],k) + mod)%mod;
    if(maxChd[vtx]){
        top[maxChd[vtx]]=top[vtx];
        DFS2(maxChd[vtx]);
    }
    for(int_t to:graph[vtx]) if(to!=maxChd[vtx]){
        top[to]=to;
        DFS2(to);
    }
}
int_t query(int_t v1){
    int_t result=0;
    while(top[v1]){
        result+=root->query(DFSN[top[v1]],DFSN[v1]);
        result%=mod;
        v1=parent[top[v1]];
    }
    return result;
}
void add(int_t v1,int_t val=1){
    while(top[v1]){
        root->add(DFSN[top[v1]],DFSN[v1],val);
        v1=parent[top[v1]];
    }
}
int_t result[LARGE+1];
struct Query{
    int_t prev, z;
    int_t* result;
};
std::vector<Query> querys[LARGE+1];
int main() {
    scanf("%lld%lld%lld",&n,&q,&k);
    for(int_t i=2;i<=n;i++){
        int_t x;
        scanf("%lld",&x);
        graph[x].push_back(i);
    }
    for(int_t i=1;i<=q;i++){
        int_t x,z;
        scanf("%lld%lld",&x,&z);
        querys[x].push_back(Query{x,z,&result[i]});
    }
    DFS1(1);
    top[1]=1;
    DFS2(1);
    #ifdef DEBUG
    for(int_t i=1;i<=n;i++) cout<<"top "<<i<<" = "<<top[i]<<endl;
    #endif
    
    root=Node::build(1,n,target);
    for(int_t i=1;i<=n;i++){
        add(i);
        for(const auto& query:querys[i]){
            *query.result=::query(query.z);
        }
    }
    for(int_t i=1;i<=q;i++){
        printf("%lld\n",result[i]);
    }
    return 0;
}

单调栈板子题..?

首先按位拆分，然后转变为统计一个01矩阵中所有全1的子矩阵（按位与）以及所有全0的子矩阵（用总矩阵数减掉即为按位或的答案）。

打表可知$n\times n$的矩阵有$(\frac{n(n+1)}2)^2$个子矩阵。

统计全1矩阵和全0矩阵本质一样，故我们先考虑统计全0矩阵。

考虑处理出来$f(i,j)$，表示如果第i行第j列的元素时0，那么往上到什么位置的元素也全都是0。对于是1的位置，$f(i,j)=-1$。

然后我们枚举一个下界$low$，表示现在我们要统计下边缘在$low$上的矩形个数。

然后从左到右扫，维护一个单调栈，单调栈内的元素是列，设$S_i$为栈中自底向上第$i$个元素，那么满足$f(low,i+1)>f(low,i)$。

先不考虑出栈，考虑入栈。

处理右下角为$(i,j)$时的答案时，设$x$表示合法的矩形数量

假设第i列的元素$i$入栈的时候，我们设$S_{-1}$表示栈内最后一个元素，$S_{-2}$以此类推。

如果加入元素i不会使得任何元素出栈，那么我们可以得到以$(i,j)$为右下角的所有合法的矩形数量为$x+f(low,i)$。

由于加入i不会导致元素出栈，故$f(low,i)$一定大于栈中其他元素的值，所以第i列较低的$S_{-1}$列可以和前面的拼起来，比$S_{-1}$高的部分可以自己成一个宽度为1的矩形。

如果导致了元素出栈呢？

左侧到$S_{-2}$，上边界到$f(low,S_{-1})$的矩形不可能再次被取到，减掉即可。

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <vector>
#include <inttypes.h>
#define debug(x) std::cout << #x << " = " << x << std::endl;
/*
°´Î»ÌÖÂÛ
ANDֵΪ1->¾ØÕóÈ«²¿Îª1
ORֵΪ1->¾ØÕóÖÁÉÙÓÐÒ»¸ö1
ͳ¼ÆÈ«0ºÍÈ«1µÄ¾ØÕó¸öÊý£¬ÖÁÉÙÓÐÒ»¸ö1µÈ¼ÛÓÚ×Ü×Ó¾ØÕóÊý-È«0¾ØÕó¡£
*/
typedef int int_t;

using std::cin;
using std::endl;
using std::cout;
const int_t mod = 1e9 + 7;
const int_t LARGE = 1e3+3;
int_t mat[LARGE+1][LARGE+1];

//last:ij¸öλÖÃÉÏ·½µÄµÚÒ»¸ö1µÄλÖÃ
int_t count[LARGE+1][LARGE+1],last[2][LARGE+1][LARGE+1];
int_t resultAnd,resultOr;
int_t n;
int_t C2(int_t n) {
	return ((int64_t)n*(n-1)/2)%mod;
}
void countFor(int_t bit) {
	//ö¾ÙÁÐ
	for(int_t i=1; i<=n; i++) {
		last[0][0][i]=last[1][0][i]=1;
		for(int_t j=1; j<=n; j++) {
			for(int_t k=0; k<=1; k++) {
				if(count[j-1][i]!=k) {
					last[k][j][i]=j;
				} else {
					last[k][j][i]=last[k][j-1][i];
				}
				if(count[j][i]!=k) last[k][j][i]=-1;
			}
		}
	}
	int_t zero=0,one=0;
	//ö¾Ùµ×
	for(int_t low=1; low<=n; low++) {
		static int_t val[2][LARGE+1];
		for(int_t i=1; i<=n; i++) {
			for(int_t j=0; j<=1; j++) {
				if(last[j][low][i]==-1) val[j][i]=0;
				else val[j][i]=low-last[j][low][i]+1;
			}
		}
		const auto count=[&](int_t bit) {
			int_t result=0;
			int_t count = 0;
			//´Ó×óÍùÓÒɨ
			std::vector<int_t> stack;
			stack.push_back(0);
			val[bit][0]=-1;
			for(int_t i=1; i<=n; i++) {
				//ºÍÇ°ÃæÁ¬µ½Ò»ÆðµÄ & ×Ô¼ºµÄ
				count += val[bit][i];
				while(stack.size() >= 1 && val[bit][stack.back()]>=val[bit][i]) {
					count -= (stack.back() - stack[stack.size() - 2]) * (val[bit][stack.back()] - val[bit][i]);
					stack.pop_back();
				}
				stack.push_back(i);
				result += count;
			}
			return result;
		};
		zero=(zero+count(0))%mod;
		one=(one+count(1))%mod;
	}
	resultAnd=(resultAnd+(int64_t)one*(1ll<<bit)%mod)%mod;
	const auto pow2=[](int64_t x) {
		x%=mod;
		return x*x%mod;
	};
	int_t least = (pow2((int64_t)n*(n+1)/2) - zero + mod) % mod;
	resultOr=((int64_t)resultOr+least*(1ll<<bit)%mod)%mod;
}
namespace SimpleIO{
    const int_t LARGE = 5e7;
    char buf[LARGE + 1];
    char* used = buf;
    void init(){
        fread(buf,1,LARGE,stdin);
    }
    template<class T>
    void nextInt(T& x){
        x = 0;
        while(isdigit(*used) == false) used++;
        while(isdigit(*used)){
            x = x * 10 + (*used) - '0';
            used++;
        }
    }
}
int main() {
//	scanf("%d",&n);
    SimpleIO::init();
    SimpleIO::nextInt(n);
	int_t maxval=0;
	for(int_t i=1; i<=n; i++) {
		for(int_t j=1; j<=n; j++) {
//			scanf("%d",&mat[i][j]);
            SimpleIO::nextInt(mat[i][j]);
			maxval=std::max(maxval,mat[i][j]);
		}
	}
	for(int_t i=0; (1ll<<i)<=maxval; i++) {
		for(int_t j=1; j<=n; j++) {
			for(int_t k=1; k<=n; k++) {
				count[j][k]=(mat[j][k]&(1ll<<i))!=0;
			}
		}
		countFor(i);
	}
	cout << resultAnd << " " << resultOr << endl;
	return 0;
}

$$ \text{设}f\left( n \right) \text{表示}2\times n\text{的矩形的方案数。} \\ \text{设}g\left( n \right) \text{为}2\times n,\text{全部使用}1\times \text{2的砖块填充的方案数} \\ \text{转移则有}f\left( n \right) =f\left( n-1 \right) +f\left( n-2 \right) +2\sum_{1\le j\le n-2}{g\left( j-1 \right)}\left( \text{枚举在第}j\text{列放置第一个}1\times \text{1的块，第}n\text{列放置第二个} \right) \\ =f\left( n-1 \right) +f\left( n-2 \right) +\sum_{0\le j\le n-3}{g\left( j \right)} \\ \text{由}g\left( 0 \right) =g\left( 1 \right) =\text{1,}g\left( n \right) =g\left( n-1 \right) +g\left( n-2 \right) \text{知}\sum_{0\le j\le n-3}{g\left( j \right)}=g\left( n-1 \right) -1 \\ \text{故}f\left( n \right) =f\left( n-1 \right) +f\left( n-2 \right) +2g\left( n-1 \right) -2 \\ \,\,\text{构造矩阵} \\ M_1=\left[ \begin{matrix}{l} f\left( 1 \right)& f\left( 2 \right)& g\left( 1 \right)& g\left( 2 \right)& -2\\ \end{matrix} \right] \text{和转移矩阵}T=\left[ \begin{matrix}{l} 0& 1& 0& 0& 0\\ 1& 1& 0& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 2& 1& 1& 0\\ 0& 1& 0& 0& 1\\ \end{matrix} \right] \\ \text{则}M_1T^{n-1}\text{即为结果矩阵。} \\ f\left( 1 \right) =f\left( 2 \right) =\text{0,}g\left( 1 \right) =\text{1,}g\left( 2 \right) =2 \\ $$

#include <cstring>
#include <iostream>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;

const int_t mod = 1e9 + 7;

struct Matrix {
    int_t data[6][6];
    int_t n;
    Matrix(int_t n) : n(n) { memset(data, 0, sizeof(data)); }
    Matrix operator*(const Matrix& rhs) const {
        Matrix result(n);
        for (int_t i = 1; i <= n; i++) {
            for (int_t j = 1; j <= n; j++) {
                for (int_t k = 1; k <= n; k++) {
                    result.data[i][j] = (result.data[i][j] +
                                         data[i][k] * rhs.data[k][j] % mod) %
                                        mod;
                }
            }
        }
        return result;
    }
    Matrix pow(int_t index) {
        Matrix result(n), base = *this;
        for (int_t i = 1; i <= n; i++) result.data[i][i] = 1;
        while (index) {
            if (index & 1) result = result * base;
            base = base * base;
            index >>= 1;
        }
        return result;
    }
};

int main() {
    int_t T;
    cin >> T;
    Matrix trans(5);
    trans.data[1][2] = 1;
    trans.data[2][1] = 1;
    trans.data[2][2] = 1;
    trans.data[3][4] = 1;
    trans.data[4][2] = 2;
    trans.data[4][3] = 1;
    trans.data[4][4] = 1;
    trans.data[5][2] = 1;
    trans.data[5][5] = 1;
    Matrix M0(5);
    M0.data[1][3] = 1;
    M0.data[1][4] = 2;
    M0.data[1][5] = mod - 2;
    while (T--) {
        int_t n;
        cin >> n;
        if (n == 0)
            cout << 0 << endl;
        else {
            auto res = M0 * trans.pow(n - 1);
            cout << res.data[1][1] << endl;
        }
    }
    return 0;
}