Xmas Contest 2021 选做
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 中国国家集训队论文
