$$ \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;
}