Xmas Contest 2021
Xmas Contest 其實已經在 AtCoder 舉辦了幾年了，不過不能直接通過首頁進入，而且題目沒有英文。\(\newcommand{\me}{\mathrm e}\)
C - Count Me
題目大意
給定一個 0、1、? 構成的字符串。將每個 ? 枚舉其填 0 或 1 后得到一個 01 串 \(S\)，對所有 \(S\) 的 \(f(S)\) 求和，\(f(S)\) 為如下序列的數量：
- \(S_0,S_1,\dots,S_N=S\) 皆為
01序列，其中第 \(k\) 項的長度為 \(k\)。 - \(S_k\) 是 \(S_{k+1}\) 的子序列。
對 \(998244353\) 取模。
- \(N\le 2.5\times 10^5\)
經過一些特殊情況的打表我們猜想，答案等價於如下 \(N+1\) 階排列計數：將 0 替換為 <、將 1 替換為 > 之后滿足相鄰不等關系的排列。
證明是很簡單的，首先我們發現操作序列等價於每次把某個全 0 段或者全 1 段的長度減去 \(1\)，考慮每次去掉最大的那個數，這能和刪除一個 0 或 1 的操作序列進行一一映射。
那么我們直接把 LOJ575 不等關系 的代碼蒯過來就好了。
復雜度 \(O(N\log ^2N)\) 或 \(O\left(\frac{N\log^2 N}{\log\log N}\right)\) 或更低。
點擊查看代碼
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int L = 20, N = 1 << L, P = 998244353, R = 3;
int n;
int inv[N], fac[N], ifac[N], f[N];
bool key[N];
char s[N], t[N];
int norm(int x) { return x >= P ? (x - P) : x; }
int reduce(int x) { return x < 0 ? x + P : x; }
void add(int& x, int y) {
if ((x += y) >= P)
x -= P;
}
int mpow(int x, int k) {
int ret = 1;
while (k) {
if (k & 1)
ret = ret * (ll)x % P;
k >>= 1;
x = x * (ll)x % P;
}
return ret;
}
void prepare(int n) {
inv[1] = 1;
for (int i = 2; i <= n; ++i) inv[i] = -(P / i) * (ll)inv[P % i] % P + P;
fac[0] = 1;
for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * (ll)i % P;
ifac[0] = 1;
for (int i = 1; i <= n; ++i) ifac[i] = ifac[i - 1] * (ll)inv[i] % P;
}
struct NTT {
int L;
vector<int> root;
NTT() : L(-1) {}
void prepRoot(int l) {
L = l;
root.resize((1 << L) + 1);
int i, n = 1 << L;
int *w2 = root.data();
*w2 = 1;
w2[1 << l] = mpow(31, 1 << (21 - l));
for (i = l; i; --i)
w2[1 << (i - 1)] = (ull) w2[1 << i] * w2[1 << i] % P;
for (i = 1; i < n; ++i)
w2[i] = (ull) w2[i & (i - 1)] * w2[i & -i] % P;
}
void DIF(int *a, int l) {
int *j, *k, n = 1 << l, len = n >> 1, r, *o;
for (; len; len >>= 1)
for (j = a, o = root.data(); j != a + n; j += len << 1, ++o)
for (k = j; k != j + len; ++k) {
r = (ull) *o * k[len] % P;
k[len] = reduce(*k - r);
add(*k, r);
}
}
void DIT(int *a, int l) {
int *j, *k, n = 1 << l, len = 1, r, *o;
for (; len != n; len <<= 1)
for (j = a, o = root.data(); j != a + n; j += len << 1, ++o)
for (k = j; k != j + len; ++k) {
r = norm(*k + k[len]);
k[len] = ull(*k - k[len] + P) * *o % P;
*k = r;
}
}
void fft(int *a, int lgn, int d = 1) {
if (L < lgn) prepRoot(lgn);
int n = 1 << lgn;
if (d == 1) DIF(a, lgn);
else {
DIT(a, lgn);
reverse(a + 1, a + n);
ull nv = P - (P - 1) / n;
for (int i = 0; i < n; ++i) a[i] = a[i] * nv % P;
}
}
} ntt;
const int BRUTE_N2_LIMIT = 50;
int pool[N << 2];
void dc(int l, int r, int* top) {
if (l == r) {
if (l == 0) return;
if (key[l])
f[l] = norm(P - f[l]);
else
f[l] = 0;
return;
}
if (r - l <= BRUTE_N2_LIMIT) {
for (int i = l; i <= r; ++i) {
for (int j = l; j < i; ++j)
f[i] = (f[i] + f[j] * (ll)ifac[i - j]) % P;
dc(i, i, top);
}
return;
}
int lg = 31 - __builtin_clz(r - l);
int d = (r - l) / lg + 1;
int lgd = 0;
while ((1 << lgd) < d) ++lgd;
++lgd;
for (int i = 0; i < lg; ++i) {
copy(ifac + i * d, ifac + (i + 2) * d, top + (i << lgd));
fill(top + (i << lgd) + 2 * d, top + ((i + 1) << lgd), 0);
ntt.fft(top + (i << lgd), lgd, 1);
}
fill(top + (lg << lgd), top + (lg << (lgd + 1)), 0);
for (int i = 0; i < lg; ++i) {
if (i) ntt.fft(top + ((lg + i) << lgd), lgd, -1);
for (int j = 0; j < min(d, r - l - i * d + 1); ++j)
f[l + i * d + j] = norm(f[l + i * d + j] + top[((lg + i) << lgd) + d + j]);
dc(l + i * d, min(l + (i + 1) * d - 1, r), top + (lg << (lgd + 1)));
if (i + 1 < lg) {
copy(f + l + i * d, f + l + (i + 1) * d, top + ((lg + i) << lgd));
fill(top + ((lg + i) << lgd) + d, top + ((lg + i + 1) << lgd), 0);
ntt.fft(top + ((lg + i) << lgd), lgd, 1);
}
for (int j = i + 1; j < lg; ++j) {
for (int k = 0; k < (1 << lgd); ++k)
top[((lg + j) << lgd) + k] = (top[((lg + j) << lgd) + k] + top[((j - i - 1) << lgd) + k] * (ll)top[((lg + i) << lgd) + k]) % P;
}
}
}
int solve(int n) {
for (int i = 1; i <= n; ++i)
key[i] = s[i] == '>';
key[0] = true;
key[n + 1] = true;
fill(f, f + n + 2, 0);
f[0] = 1;
dc(0, n + 1, pool);
int ans = f[n + 1];
bool par = false;
for (int i = 1; i <= n + 1; ++i)
if (key[i])
par ^= true;
if (par)
ans = norm(P - ans);
return ans;
}
int main() {
prepare((1 << L) - 1);
int n; scanf("%d%s", &n, t + 1);
int lst = 0;
int ans = fac[n + 1];
t[n + 1] = '?';
for (int i = 1; i <= n + 1; ++i) {
if (t[i] == '?') {
if (i - lst > 1) {
for (int j = lst + 1; j < i; ++j) s[j - lst] = t[j] == '0' ? '>' : '<';
ans = ans * (ll)solve(i - lst - 1) % P;
}
lst = i;
}
}
printf("%d\n", ans);
return 0;
}
D - Determinant?
題目大意
給 \(N\) 個 \(K\times K\) 的矩陣 \(A_i\)，計算 \(\sum_{\sigma\in \mathfrak S_N}\operatorname{sgn}(\sigma)A_{\sigma(1)}\cdots A_{\sigma(N)}\)，答案對 \(998244353\) 取模。
- \(1\le N\le 32\)
- \(1\le K\le 8\)
我們先從簡單的問題想起。假設 \(K=1\)，\(A\) 退化為標量，這時，當 \(N\ge 2\) 時，奇偶排列配對得知答案必為 \(0\)。
考慮計算最終矩陣的 \(M_{i_0i_N}\) 項，將其貢獻拆解為 \(A_{\sigma(1)i_0i_1}\cdots A_{\sigma(N)i_{N-1}i_N}\)。此時，對於一條具體的路徑 \(i_0 \to i_1 \to \cdots \to i_N\)，設 \(B_{jk} = A_{j i_{k-1}i_k}\)，那么這部分貢獻的求和就是
於答案 \(M_{i_0i_N}\)，我們應該可以將其寫作
那么問題來了，對於一個給定的邊集，后面這個值有什么樣的性質呢？我們注意到，如果一個點有兩次從它出發又回到自身，那么這兩個回路的邊順序交換之后，\(\sigma\)>那么問題來了，對於一個給定的邊集，后面這個值有什么樣的性質呢？我們注意到，如果一個點有兩次從它出發又回到自身，那么這兩個回路的邊順序交換之后，\(\sigma\) 的奇偶性改變。因此，起點的入度 \(\ge2\) 的話就會貢獻抵消，其它點的入度 \(\ge 3\) 的話就會貢獻抵消，因此邊數不能超過 \(2K-1\)。
上面這個證明好像漏了環長是偶數的情況, 我有時間再修修, 但是好像有些代數證明比如說 這個.
所以，當 \(N\ge 2K\) 的時候答案一定是 \(O\)，否則我們直接狀壓 DP 就好了，復雜度為 \(O(2^N NK^3) = O(4^KK^4)\)。
C T T 2 0 1 9 情 景 再 現
點擊查看代碼
#include <bits/stdc++.h>
using namespace std;
using ull = unsigned long long;
using ll = long long;
const int P = 998244353;
int norm(int x) { return x >= P ? x - P : x; }
int reduce(int x) { return x < 0 ? x + P : x; }
int neg(int x) { return x ? P - x : 0; }
void add(int& x, int y) { if ((x += y - P) < 0) x += P; }
void sub(int& x, int y) { if ((x -= y) < 0) x += P; }
void fam(int& x, int y, int z) { x = (x + y * (ull)z) % P; }
int mpow(int x, unsigned k) {
if (k == 0) return 1;
int ret = mpow(x * (ull)x % P, k >> 1);
if (k & 1) ret = ret * (ull)x % P;
return ret;
}
int quo2(int x) { return ((x & 1) ? x + P : x) >> 1; }
const int _ = 8;
int dp[1 << 15][_][_];
int A[15][_][_];
int main() {
ios::sync_with_stdio(false); cin.tie(NULL);
int N, K; cin >> N >> K;
if (N >= K * 2) {
for (int i = 0; i != K; ++i)
for (int j = 0; j != K; ++j)
cout << 0 << " \n"[j == K - 1];
return 0;
}
for (int i = 0; i != N; ++i)
for (int j = 0; j != K; ++j)
for (int k = 0; k != K; ++k)
cin >> A[i][k][j];
for (int i = 0; i != K; ++i) dp[0][i][i] = 1;
for (int s = 1; s != 1 << N; ++s)
for (int i = 0; i != N; ++i) if ((s >> i) & 1) {
if (__builtin_parity(s >> (i + 1))) {
for (int j = 0; j != K; ++j)
for (int k = 0; k != K; ++k)
for (int l = 0; l != K; ++l)
fam(dp[s][j][k], dp[s ^ (1 << i)][j][l], P - A[i][k][l]);
} else {
for (int j = 0; j != K; ++j)
for (int k = 0; k != K; ++k)
for (int l = 0; l != K; ++l)
fam(dp[s][j][k], dp[s ^ (1 << i)][j][l], A[i][k][l]);
}
}
for (int i = 0; i != K; ++i)
for (int j = 0; j != K; ++j)
cout << dp[(1 << N) - 1][i][j] << " \n"[j == K - 1];
return 0;
}
E - E and PI
題目大意
對於正整數 \(j\)，數軸上點 \(P_j\) 的坐標為 \(j\me \bmod 1\)，其中 \(\me\) 是自然對數的底。對於正整數 \(n\)，考慮點 \(P_1,\dots,P_{n-1}\) 將線段 \([0,1]\) 分割成 \(n\) 段，記 \(f(n)\) 為這 \(n\) 段中的最大值。
記 \(n_K\) 是第 \(K\) 個滿足 \(f(n)>f(n+1)\) 的正整數，給定 \(K\)，求出 \(n_K\) 和 \(f(n_K)\)。顯然 \(f(n_K)\) 可以表示成 \(a + b\me\) 的形式，只需要輸出 \(n_K,a,b\) 同余 \(998244353\) 的結果。
- \(1\le K\le 10^{11}\)
由於本人不太懂具體道理，所以這里主要講講大概找到什么規律。
經過打表和 OEIS 搜索之后我們可以發現，\(n_K = \texttt{A006259}(K+3)-1\)，這個數列由 \(\me\) 的最佳有理逼近給出。也就是，將 \(\me\) 在 Stern-Brocot 樹上走的字符串 \(\me = \mathrm{RL^0RLR^2LRL^4RLR^6}\dots\) 截斷到前 \(K+2\) 個字符，得到的數 \(-1\) 就是答案：
顯然我們可以在 \(O(1)\) 的時間內完成 \(\mathrm{L}^n\) 或 \(\mathrm{R}^n\) 的計算，因此我們就可以在 \(O(\sqrt K)\) 的時間內暴力計算出 \(n_K\) 了。
然后經過打表，我們又發現對應的 \(f(n_K)\) 是這樣的：
我們發現除了第一個，對於剩下的 \(a+b\me\)，\(\dfrac {-a}b\) 的值依然都是出現在有理逼近里的！接着我們更具體地看，容易想到 \(a,b\) 何者為負只取決於 \(\dfrac ab\) 和 \(\me\) 的大小關系，也就是說，如果 \(\dfrac ab\) 的下一個字符是 \(\mathrm L\)，那么說明 \(\dfrac ab > \me\)，出現在答案里的應該是 \(a - b\me\)，反之則是 \(-a+b\me\)。再仔細看，我們發現有些數出現的位置相比於原序列是變得靠后了：
看起來規律還是不夠顯然，我們再把 \(\mathrm{LR}\) 標出來，並且用紅色表示位置不對的數，灰色表示這些數在逼近序列里的位置。
如果用每個數接下來的字符代表它，那么我們發現這個情況就是它恰好把前面提到的 \(\mathrm{LR}^{2n}/\mathrm{RL}^{2n}\)，替換成了 \(\mathrm{{\color{gray}L}R}^{2n-1}\mathrm{\color{red}L}\mathrm{R}/\mathrm{{\color{gray}R}L}^{2n-1}\mathrm{\color{red}R}\mathrm{L}\)。
這也是好處理的，於是我們就可以 \(O(\sqrt K)\) 解決這道題了。
（那么如果我們再用類似整式遞推的方法加個速，可以做到 \(O(K^{1/4}\log K)\)，這題就能出到 \(10^{18}\) 了）
點擊查看代碼
#include <bits/stdc++.h>
using namespace std;
using ull = unsigned long long;
const int P = 998244353;
int norm(int x) { return x >= P ? x - P : x; }
int reduce(int x) { return x < 0 ? x + P : x; }
int neg(int x) { return x ? P - x : 0; }
void add(int& x, int y) { if ((x += y - P) < 0) x += P; }
void sub(int& x, int y) { if ((x -= y) < 0) x += P; }
void fam(int& x, int y, int z) { x = (x + y * (ull)z) % P; }
int mpow(int x, unsigned k) {
if (k == 0) return 1;
int ret = mpow(x * (ull)x % P, k >> 1);
if (k & 1) ret = ret * (ull)x % P;
return ret;
}
int quo2(int x) { return ((x & 1) ? x + P : x) >> 1; }
void cal1(ull K) {
int a = 1, b = 0;
auto apply = [&](bool f, int N) {
N = min((ull)N, K);
K -= N;
if (!f) fam(a, b, N);
else fam(b, a, N);
return K == 0;
};
for (int i = 0; ; ++i) {
if (apply(i & 1, 1)) break;
if (apply(~i & 1, i * 2)) break;
if (apply(i & 1, 1)) break;
}
cout << norm(a + b - 1);
}
int main() {
ios::sync_with_stdio(false);
cin.tie(NULL); cout.tie(NULL);
ull K; cin >> K;
cal1(K + 2);
if (K == 1) { cout << " 1 0\n"; return 0; }
if (K == 2) { cout << " 998244351 1\n"; return 0; }
--K;
int x = 0, y = 1, a = 1, b = 0;
auto apply = [&](bool f, int N) {
N = min((ull)N, K);
K -= N;
if (!f) { fam(a, b, N); fam(x, y, N); }
else { fam(b, a, N); fam(y, x, N); }
return K == 0;
};
for (int i = 0; ; ++i) {
if (K >= (i + 1) * 2) {
apply(i & 1, 1); apply(~i & 1, i * 2), apply(i & 1, 1);
continue;
}
if (K < i * 2 - 1) ++K;
else if (K == i * 2 - 1) K = 0;
bool flag = (K == 0 || K == i * 2 + 1) ^ (i & 1);
int a1, a2;
if (!apply(i & 1, 1))
if (!apply(~i & 1, i * 2))
apply(i & 1, 1);
a1 = norm(x + y); a2 = norm(a + b);
if (flag) a1 = reduce(-a1);
else a2 = reduce(-a2);
cout << ' ' << a1 << ' ' << a2 << '\n';
break;
}
return 0;
}
H - Homework from Zhejiang
題目大意
給定大小分別為 \(MN\) 的有向帶權圖 \(G_1,G_2\)，權值矩陣分別為 \(A,B\)。令 \(G = G_1 \square G_2\)，也即點集由 \(V_1 \times V_2\) 構成，其中 \((i,j)\) 到 \((k,j)\) 有邊，權值為 \(A_{ik}\)，另外 \((i,j)\) 到 \((i,l)\) 有邊，權值為 \(B_{jl}\)。
對於每個 \((i,j)\)，求以 \((i,j)\in G_1 \square G_2\) 為根的外向樹權值和，權值定義為邊權之積。答案對 \(998244353\) 取模。
- \(2\le M,N\le 500\)
這就是經典的圖的 Cartesian 積。令 \(L^{(1)},L^{(2)}\) 分別為 \(G_1,G_2\) 的 Laplacian 矩陣，易得 \(G_1 \mathop \square G_2\) 的 Laplacian 矩陣為 \(L=L^{(1)} \otimes I_N + I_M \otimes L^{(2)}\)。我們知道，計算以一個點 \(u\) 為根的生成樹數量就是在 \(L_{uu}\) 處的主余子式，只需計算 Laplacian 矩陣的伴隨矩陣。由於 Laplacian 矩陣不滿秩，伴隨矩陣能寫作 \(x y^{\mathsf T}\) 的形式，其中 \(Lx=0\)，\(y^{\mathsf T}L=0\)，根據 Laplacian 矩陣的性質，顯然 \(y\) 只需取全 \(1\) 向量，而我們可以在 \(O(N^3)\) 時間內消元得到 \(x\)。得到了 \(x^{(1)},x^{(2)}\) 之后，易見 \(L(x^{(1)}\otimes x^{(2)})=0\)，可知 \(x^{(1)}\otimes x^{(2)}\) 就是答案了。但其實還有點問題，求出來的向量其實整體差一個常數倍，我們還需要定出這個常數倍。
設 \(L^{(1)},L^{(2)}\) 的特征值分別為 \(\lambda_1,\dots,\lambda_M\) 和 \(\mu_1,\dots,\mu_N\)，我們有 \(L\) 的特征值為 \(\lambda_i + \mu_j\)（證明見 [1]），由於一定有 \(0\) 這個特征值，不妨設 \(\lambda_1 = \mu_1 = 0\)，那么 \(G\) 以每個點為根的生成樹數量之和，也就是特征多項式的 \(1\) 次項系數，也即求出除 \(\lambda_1 + \mu_1\) 以外的特征值乘積。這可以用求出 \(L^{(1)}\) 和 \(L^{(2)}\) 的特征多項式之后利用結式計算，答案記為 \(T\)。這樣，我們乘以常數 \(c\) 使得 \(c (x^{(1)}\otimes x^{(2)})^{\mathsf T} \cdot 1 = T\) 之后，答案就對了。
由於答案之和有可能剛好是 \(998244353\) 的倍數，上述做法會寄，我們需要稍微仔細一點處理。
設 \(x^{(1)}, x^{(2)}\) 分別為 \(G_1,G_2\) 的答案向量，這是容易算的，首先得到他們的常數倍，然后對一個非零位置求一次行列式就能得到。那么注意到他們的和分別是 \(\lambda_2 \cdots \lambda_M\) 和 \(\mu_2\cdots \mu_N\)，所以此時 \(G\) 的答案向量一定是
而此時的乘積就是 \(\operatorname{res}\left(\frac{|tI - A|}t, \frac{|(-t)I-B|}{-t}\right)\)。
復雜度 \(O(M^3 + N^3)\)。
點擊查看代碼
#include <bits/stdc++.h>
using namespace std;
using ull = unsigned long long;
using ll = long long;
const int P = 998244353;
int norm(int x) { return x >= P ? x - P : x; }
int reduce(int x) { return x < 0 ? x + P : x; }
int neg(int x) { return x ? P - x : 0; }
void add(int& x, int y) { if ((x += y - P) < 0) x += P; }
void sub(int& x, int y) { if ((x -= y) < 0) x += P; }
void fam(int& x, int y, int z) { x = (x + y * (ull)z) % P; }
int mpow(int x, unsigned k) {
if (k == 0) return 1;
int ret = mpow(x * (ull)x % P, k >> 1);
if (k & 1) ret = ret * (ull)x % P;
return ret;
}
int quo2(int x) { return ((x & 1) ? x + P : x) >> 1; }
namespace PinkRabbit {
const int MN = 510, Mod = P;
typedef long long LL;
int N, A[MN][MN], H[MN][MN], Ans[MN];
void Hessenberg() {
for (int i = 1; i <= N; ++i)
for (int j = 1; j <= N; ++j) {
H[i][j] = A[i][j];
}
for (int i = 1; i <= N - 2; ++i) {
if (!H[i + 1][i]) for (int j = i + 2; j <= N; ++j) if (H[j][i]) {
for (int k = i; k <= N; ++k) std::swap(H[i + 1][k], H[j][k]);
for (int k = 1; k <= N; ++k) std::swap(H[k][i + 1], H[k][j]);
break;
}
if (!H[i + 1][i]) continue;
int val = mpow(H[i + 1][i], P - 2);
for (int j = i + 2; j <= N; ++j) {
int coef = (LL)val * H[j][i] % Mod;
for (int k = i; k <= N; ++k) H[j][k] = (H[j][k] + (LL)H[i + 1][k] * (Mod - coef)) % Mod;
for (int k = 1; k <= N; ++k) H[k][i + 1] = (H[k][i + 1] + (LL)H[k][j] * coef) % Mod;
}
}
}
void Characteristic() {
Hessenberg();
for (int i = 1; i <= N; ++i)
for (int j = 1; j <= N; ++j)
H[i][j] = Mod - H[i][j];
static int P[MN][MN];
P[0][0] = 1;
for (int i = 1; i <= N; ++i) {
P[i][0] = 0;
for (int j = 1; j <= i; ++j) P[i][j] = P[i - 1][j - 1];
int val = 1;
for (int j = i - 1; j >= 0; --j) {
int coef = (LL)val * H[j + 1][i] % Mod;
for (int k = 0; k <= j; ++k) P[i][k] = (P[i][k] + (LL)P[j][k] * coef) % Mod;
if (j) val = (LL)val * (Mod - H[j + 1][j]) % Mod;
}
}
for (int i = 0; i <= N; ++i) Ans[i] = P[N][i];
}
}
const int _ = 510;
namespace Gauss {
int bas[_][_], mat[_][_];
void solve(int N, int adj[]) {
memset(bas, 0, sizeof(bas));
for (int i = 1; i <= N; ++i) bas[i][i] = 1;
for (int i = 1; i <= N; ++i) {
int piv = 1;
while (piv <= N && !mat[i][piv]) ++piv;
if (piv > N) {
memcpy(adj, bas[i], sizeof(bas[i]));
break;
}
int c = mpow(mat[i][piv], P - 2);
for (int j = piv; j <= N; ++j)
mat[i][j] = mat[i][j] * (ull)c % P;
for (int j = 1; j <= i; ++j)
bas[i][j] = bas[i][j] * (ull)c % P;
for (int j = i + 1; j <= N; ++j) {
c = P - mat[j][piv];
for (int k = piv; k <= N; ++k)
fam(mat[j][k], mat[i][k], c);
for (int k = 1; k <= i; ++k)
fam(bas[j][k], bas[i][k], c);
}
}
}
}
struct G {
int N;
int mat[_][_];
int adj[_];
vector<int> charp;
int det() {
int ret = 1;
for (int i = 1; i <= N; ++i) {
int pivot = i;
while (!mat[pivot][i] && pivot <= N) ++pivot;
if (pivot > N) return 0;
if (pivot != i) {
ret = P - ret;
for (int j = i; j <= N; ++j) std::swap(mat[i][j], mat[pivot][j]);
}
ret = ret * (ull)mat[i][i] % P;
int nv = mpow(mat[i][i], P - 2);
for (int j = i; j <= N; ++j) mat[i][j] = mat[i][j] * (ull)nv % P;
for (int j = i + 1; j <= N; ++j) for (int k = N; k >= i; --k)
fam(mat[j][k], P - mat[j][i], mat[i][k]);
}
return ret;
}
void solve() {
for (int i = 1; i <= N; ++i)
for (int j = 1; j <= N; ++j) {
int x; cin >> x;
add(mat[j][j], x);
sub(mat[i][j], x);
}
PinkRabbit::N = N;
memcpy(PinkRabbit::A, mat, sizeof(mat));
PinkRabbit::Characteristic();
charp.resize(N + 1);
memcpy(charp.data(), PinkRabbit::Ans, sizeof(int) * (N + 1));
for (int i = 1; i <= N; ++i)
for (int j = 1; j != i; ++j)
swap(mat[i][j], mat[j][i]);
memcpy(Gauss::mat, mat, sizeof(mat));
Gauss::solve(N, adj);
int piv = 1; while (!adj[piv]) ++piv;
for (int i = 1; i <= N; ++i) mat[piv][i] = mat[i][piv] = 0;
mat[piv][piv] = 1;
int r = det() * (ull)mpow(adj[piv], P - 2) % P;
for (int i = 1; i <= N; ++i) adj[i] = adj[i] * (ull)r % P;
}
} A, B;
int res(vector<int> p, vector<int> q) {
int ret = 1, N = p.size() - 1, M = q.size() - 1;
if (N < M) {
swap(N, M); swap(p, q);
if ((N * M) % 2) ret = reduce(-ret);
}
while (M) {
for (int i = N - M; i >= 0; --i) {
int c = reduce(-p[M + i]);
for (int j = 0; j <= M; ++j)
fam(p[j + i], q[j], c);
}
while (N && p[N] == 0) --N; p.resize(N + 1);
if (p[N] == 0) return 0;
ret = ret * (ull)mpow(p[N], M) % P;
int c = mpow(p[N], P - 2);
for (int i = 0; i <= N; ++i) p[i] = p[i] * (ull)c % P;
swap(N, M); swap(p, q);
if ((N * M) % 2) ret = reduce(-ret);
}
return ret;
}
int main() {
ios::sync_with_stdio(false); cin.tie(NULL);
cin >> A.N >> B.N;
A.solve(); B.solve();
vector<int> p = A.charp, q = B.charp;
for (int i = p.size() - 2; i >= 0; i -= 2)
p[i] = reduce(-p[i]);
p.erase(p.begin()); q.erase(q.begin());
int r = res(p, q);
for (int i = 1; i <= A.N; ++i)
for (int j = 1; j <= B.N; ++j)
cout << (A.adj[i] * (ull)B.adj[j] % P * r % P) << " \n"[j == B.N];
return 0;
}
[1] 潘佳奇，淺談線性代數與圖論的關系，IOI 2021 中國國家集訓隊論文