常系数线性齐次递推

//矩阵类，支持矩阵的加减乘和幂运算
//------------------------------------------------------------------------------------
const int MAXN = 105;
const int MAXM = 105;
struct Martix
{
    int n, m;                    //n:矩阵行数         m:矩阵列数
    int a[MAXN][MAXM];
    void clear() {            //清空矩阵
        n = m = 0;
        memset(a, 0, sizeof(a));
    }
    Martix operator+(const Martix& b) const {
        Martix tmp;
        tmp.n = n; tmp.m = m;
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                tmp.a[i][j] = a[i][j] + b.a[i][j];
        return tmp;
    }
    Martix operator-(const Martix& b) const {
        Martix tmp;
        tmp.n = n; tmp.m = m;
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                tmp.a[i][j] = a[i][j] - b.a[i][j];
        return tmp;
    }
    Martix operator*(const Martix& b) const {
        Martix tmp;
        tmp.clear();
        tmp.n = n; tmp.m = b.m;
        for (int i = 0; i < tmp.n; i++)
            for (int j = 0; j < tmp.m; j++) 
                for (int k = 0; k < m; k++) {
                    tmp.a[i][j] += a[i][k] * b.a[k][j];
                    tmp.a[i][j] %= mod;
                }
        return tmp;
    }
    Martix operator^(int x)const {
        Martix base, ans;
        base = *this; 
        ans.clear(); ans.m = ans.n = n;
        for (int i = 0; i < n; i++) ans.a[i][i] = 1;
        if (x == 0) return ans;
        if (x == 1) return base;
        while (x) {
            if (x & 1) ans = ans * base;
            base = base * base;
            x >>= 1;
        }
        return ans;
    }
};
//------------------------------------------------------------------------------------
//常系数线性齐次递推模板
//a[]为常系数数组，b为初值数组，n为数组大小，t为要求解的项数
int solve(int a[],int b[],int n,int t){
    Martix M,F,E;
    M.clear(),F.clear(),E.clear();
    M.n=M.m=n;
    E.n=E.m=n;
    F.n=n,F.m=1;
    for(int i=0;i<n-1;i++)
    M.a[i][i+1]=1;
    for(int i=0;i<n;i++){
        M.a[n-1][i]=a[i];
        F.a[i][0]=b[i];
        E.a[i][i]=1;
    }
    if(t<n) return F.a[t][0];
    for(t-=n-1;t;t/=2){
        if(t&1) E=M*E;
        M=M*M;
    }
    F=E*F;
    return F.a[n-1][0];
}