矩阵 - HDOJ1575 - TrA

模板裸题
矩阵快速幂
借用了 zngg 的矩阵乘，写个快速幂就好了

#include <bits/stdc++.h>
using namespace std;

#define ll long long

const int MAX_MAT = 20;
const ll mod = 9973;
int T, n, k;

struct Mat{
    ll a[MAX_MAT][MAX_MAT];
    Mat(){
        for(int i = 0; i < MAX_MAT; i++){
            for(int j = 0; j < MAX_MAT; j++)
                a[i][j] = 0;
            a[i][i] = 1;
        }
    }
};

ll quickpow(ll x, ll y, ll MOD){
    ll ans = 1;
    while(y){
        if (y & 1)
            ans = (x * ans) % MOD;
        x = (x * x) % MOD;
        y >>= 1;
    }
    return ans;
}

ll A[MAX_MAT][MAX_MAT << 1];
ll get_inv(ll x){
    return quickpow(x, mod - 2, mod);
} 

void row_minus(int a, int b, ll k){
    for(int i = 0; i < 2 * MAX_MAT; i++)
        A[a][i] = (A[a][i] - A[b][i] * k % mod) % mod, A[a][i] = (A[a][i] + mod) % mod;
}

void row_multiplies(int a, ll k){
    for(int i = 0; i < 2 * MAX_MAT; i++)
        A[a][i] = (A[a][i] * k) % mod;
}

void row_swap(int a, int b){
    for(int i = 0; i < 2 * MAX_MAT; i++)
        swap(A[a][i], A[b][i]);
}

Mat getinv(Mat x){
    memset(A, 0, sizeof(A));
    for(int i = 0; i < MAX_MAT; i++)
        for(int j = 0; j < MAX_MAT; j++)
            A[i][j] = x.a[i][j], A[i][MAX_MAT + j] = (i == j);
    for(int i = 0; i < MAX_MAT; i++){
        if (!A[i][i]){
            for(int j = i + 1; j < MAX_MAT; j++)
                if (A[j][i]){
                    row_swap(i, j);
                    break;
                }
        }
        row_multiplies(i, get_inv(A[i][i]));
        for(int j = i + 1; j < MAX_MAT; j++)
            row_minus(j, i, A[j][i]);
    }
    for(int i = MAX_MAT - 1; i >= 0; i--)
        for(int j = i - 1; j >= 0; j--)
            row_minus(j, i, A[j][i]);
    Mat ret;
    for(int i = 0; i < MAX_MAT; i++)
        for(int j = 0; j < MAX_MAT; j++)
            ret.a[i][j] = A[i][MAX_MAT + j];
    return ret;
}

Mat operator * (Mat x, Mat y){
    Mat c;
    for(int i = 0; i < MAX_MAT; i++)
        for(int j = 0; j < MAX_MAT; j++)
            c.a[i][j] = 0;
    for(int i = 0; i < MAX_MAT; i++)
        for(int j = 0; j < MAX_MAT; j++)
            for(int k = 0; k < MAX_MAT; k++)
                c.a[i][j] = (c.a[i][j] + x.a[i][k] * y.a[k][j] % mod) % mod;
    return c;
}

Mat Mat_quickpow(Mat A, ll k){
    Mat ret;
    while(k){
        if (k & 1)
            ret = ret * A;
        A = A * A;
        k >>= 1;
    }
    return ret;
}

int main(){
    //freopen("input.txt", "r", stdin);
    scanf("%d", &T);
    Mat A, Ak;
    ll sum;
    while(T--){
        scanf("%d%d", &n, &k);
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                scanf("%lld", &A.a[i][j]);
        Ak = Mat_quickpow(A, k);
        sum = 0;
        for(int i = 0; i < n; i++)
            sum = (sum + Ak.a[i][i]) % mod;
        printf("%lld\n", sum);
    }
    return 0;
}