Goodbye Renyin E 新年的找对题解

前排提醒：如果数学公式看不清楚，可以把鼠标悬停上去，会放大。

简要题意

给定一张 $n$ 个点的简单无向图，每个点可能有一个颜色 $c_{i} > 0$ , 也可能没有颜色，记作 $c_{i} = 0$ .

你需要找出尽可能多的简单路径，使得：

这些路径间两两不交（不经过重复点）；
路径的起点与终点需要是有颜色但是颜色不同的点；
除了起点终点外，路径只能经过没有颜色的点。

需要构造出答案。

$n ⩽ 300$ .

背景与参考文献

这个问题是经典问题，叫做 S-path problem。

本文将陈述此问题的 $O (n^{3})$ 算法。若使用更快的 $O (n^{w})$ 的矩阵乘法，则可以优化到 $O (n^{w})$ 。若将来 $w = 2$ ，则复杂度为 $O (n^{2} \log n)$ 。

参考文献为

Algebraic Algorithms for Linear Matroid Parity Problems, Ho Yee Cheung, Lap Chi Lau, Kai Man Leung

理论与算法

1. 问题转化

我们先把原问题转化为：

问题一：选出原图中尽量多的边，使得这些边构成森林，且森林中每棵树里最多有两个有颜色的点，并且不能有两个相同颜色的点。

这样的话，最优解中每棵树里肯定要么有一个有颜色的点，要么有两个有颜色且颜色不同的点（除非原图中某个连通块里一个有颜色的点都没有，这样的连通块显然可以直接扔掉）。

所以原问题的答案一定就是总的有颜色的点的个数，减去新问题中连通块个数。而新问题中连通块个数又等于 $n$ 减去选出的边数。

因此选的边越多，得出的简单路径就越多。这些简单路径也可以由新问题里选出的边轻易构造出来。

这一步转化是容易的，我们接下来会将这个新问题再转化为一个线性代数问题。

2. 线性拟阵配对问题

这个线性代数问题如下。

问题二(线性拟阵配对)：给定 $N$ 维空间里的 $M$ 对向量，即 $2 M$ 个向量 $(b_{1}, c_{1}), (b_{2}, c_{2}), \dots, (b_{M}, c_{M})$ . 你需要找出其中尽可能多对，使得这些对向量线性无关。

如果你并不知道“线性无关”的定义，建议参考OI-Wiki: 线性基。 此外，OI-Wiki 中线性代数相关的其他概念也可能对你理解本文有帮助。

如何将前面的问题一转化为该问题呢？有一个非常精妙的构造如下。

我们先取一个大空间 $R^{2 n}$ ，其中每个向量 $x$ 的坐标记作 $x_{1}, x_{2}, \dots, x_{2 n}$ 。此外我们记 $e_{i}$ 为第 $i$ 个分量是 $1$ , 别的分量是 $0$ 的向量。

为了方便，我们记 $x (u) = (x_{2 u - 1}, x_{2 u})$ 为一个二维向量。

接下来对原图中每条边 $(u, v)$ ，我们将其对应一个二维子空间 $L_{u, v} = {x \in R^{2 n} ∣ x (u) + x (v) = 0 且 \forall w \notin {u, v}, x (w) = 0}$ 也可以说是对应了两个向量 $b_{u, v} = e_{2 u - 1} - e_{2 v - 1}, c_{u, v} = e_{2 u} - e_{2 v}$ 。

如果只是这样的话，那么只有当选出的边构成一个环的时候，才会线性相关：

比如如果选出了 $x - y - z - x$ 这样三条边，那么 $b_{x, y} + b_{y, z} + b_{z, x} = 0, c_{x, y} + c_{y, z} + c_{z, x} = 0$ 。

为了满足我们的额外条件（每棵树最多有两个有颜色的点，并且颜色不能相同），还需要引入一个额外的东西。

对于每个颜色，我们选取 $R^{2}$ 里的一个非零向量 $f_{c}$ ，并且他们两两不平行（即不存在 $c \neq c^{'}$ 使得 $f_{c}$ 是 $f_{c^{'}}$ 的若干倍）。此外记 $f_{0} = 0$ 。

对每个点 $u$ ，我们取一个向量 $a_{u}$ ，使得 $a_{u} (u) = f_{c_{u}}$ ，而 $a_{u} (v) = 0$ .

也就是说 $a_{u}$ 的 $u$ 分量是二维向量 $f_{u}$ ，其他分量是 $0$ 。如果 $c_{u} = 0$ 那么也有 $a_{u} = 0$ .

可以发现，如果我们选出了一些边，这些边连通了某两个点 $s, t$ ，那么就可以用这些边对应的向量加出 $L_{s, t}$ 的元素，即使 $(s, t)$ 这条边不存在或者没被选择。

那么如果我们选出的边连接了两个颜色相同的点，比如 $1$ 和 $2$ ，就有 $a_{1} - a_{2} \in L_{1, 2}$ 。

而如果我们选出的边连接了三个有颜色的点，比如 $1, 2, 3$ 且其颜色也是 $1, 2, 3$ ，我们可以找出三个实数 $x_{1}, x_{2}, x_{3}$ 使得 $x_{1} f_{1} + x_{2} f_{2} + x_{3} f_{3} = 0$ 。那么总可以找出 $L_{1, 2}$ 和 $L_{1, 3}$ 里的元素，加起来恰好得到 $x_{1} f_{1} + x_{2} f_{2} + x_{3} f_{3}$ 。

也就是说，如果我们强制选择 $f_{1}, \dots, f_{n}$ 这些向量（除了其中的零向量），接下来再选择一些边（以及他们对应的若干对向量 $(b_{u, v}, c_{u, v})$ ），那么选出的所有向量线性无关当且仅当这些边满足问题一中的条件。

怎么强制选择这些向量呢？可以考虑把所有 $b, c$ 向量对都投影到某个和他们正交的子空间上，这样就相当于已经选择这些向量了。

3. 转化的具体实现

首先，我们先（假装）求出这 $2 m$ 个向量 $b_{u, v} = e_{2 u - 1} - e_{2 v - 1}, c_{u, v} = e_{2 u} - e_{2 v - 1}$ .

接下来我们做投影。注意我们其实不需要正交，只需要把他们投影到某个子空间 $W$ 上，满足 $R^{2 n} = W \oplus Span (f_{1}, \dots, f_{n})$ （直和）即可。如果用 OIer 比较熟悉的说法，就是用这些 $f$ 对他们做高斯消元。

所以我们不妨设 $x_{c} = (1, - c)$ ，则 $f_{u} = e_{2 u - 1} - c_{u} e_{2 u}$ 。用 $f_{u}$ 对某个向量做高斯消元的话，就是把这个向量的第 $2 u$ 项系数加上 $c_{u}$ 倍的 $2 u - 1$ 项系数，然后把第 $2 u - 1$ 项系数置为 $0$ 。也就是 v[2*u] += c[u] * v[2*u - 1]; v[2*u-1] = 0;

（然后我们也可以把第 $2 u - 1$ 维直接删掉，反正所有向量的这一维都是 $0$ ）

4. 新问题的求解

接下来我们就要求解这个线性代数问题。

我们先定义一个记号：若 $S, T \subseteq {1, \dots, n}$ 是两个子集， $M$ 是 $n \times n$ 的矩阵，那么 $M_{S, T}$ 是只取行序号属于 $S$ 列序号属于 $T$ 的位置构成的 $| S | \times | T |$ 子矩阵。如果 $S =$ 全集或者 $T =$ 全集，我们也简记作 $M_{*, T}, M_{S, *}$ .

我们给出本文中四个不加证明的定理：

首先是解决线性拟阵配对问题的核心定理：

定理 [Lovász]：给定 $N$ 维空间里的 $M$ 对列向量，即 $2 M$ 个向量 $(b_{1}, c_{1}), (b_{2}, c_{2}), \dots, (b_{M}, c_{M})$ . 我们引入 $M$ 个未知数 $x_{1}, \dots, x_{M}$ ，并且在他们的分式域里考虑矩阵（即矩阵的每个分量可以是这些未知数的有理式）。定义矩阵 $Y = \sum_{i = 1}^{M} x_{i} (b_{i} \land c_{i}) .$ 其中 $b_{i} \land c_{i} = b_{i} c_{i}^{T} - c_{i} b_{i}^{T}$ 。则 $rank (Y) = 2 k$ ，其中 $k$ 是问题二的答案，即最多可以选出的线性无关向量对数。

$\land$ 读作 wedge。 其次是让我们摆脱掉矩阵分量里的未知数的引理：

定理 [Schwartz-Zippel 引理]：在模素数 $p$ 的意义下，若 $f (x_{1}, \dots, x_{n})$ 是不超过 $d$ 次的非零多元多项式，那么任取 $x_{1}, \dots, x_{n}$ 为 $0, \dots, p - 1$ 的随机值， $f (x_{1}, \dots, x_{n}) = 0$ 的概率不超过 $\frac{d}{p}$ .

然后是让我们之后算法的重点：

定理 [Harvey]：设 $M$ 是一个 $n \times n$ 可逆矩阵, $N = M^{- 1}$ 。若 $\tilde{M}$ 和 $M$ 长得很像——具体地说，有一个子集 $S$ ，使得除了 ${\tilde{M}}_{S, S} \neq M_{S, S}$ 之外 $\tilde{M}$ 的分量都和 $M$ 相等。记 $Δ = \tilde{M} - M$ 。那么

$\tilde{M}$ 可逆当且仅当 $I + Δ_{S, S} N_{S, S}$ 可逆。

若第一条成立，则 ${\tilde{M}}^{- 1} = N - N_{*, S} (I + Δ_{S, S} N_{S, S})^{- 1} Δ_{S, S} N_{S, *}$ ;

若第一条成立，则 $({\tilde{M}}^{- 1})_{T, T} = N_{T, T} - N_{T, S} (I + Δ_{S, S} N_{S, S})^{- 1} Δ_{S, S} N_{S, T}$ 。

最后是其实可能不必要的一个玩意。

定理：若 $A$ 是斜对称矩阵，即 $A^{T} = - A$ ，而 $rank (A) = r$ ，则存在 $| S | = r$ 使得 $A_{S, S}$ 可逆。具体地，只要 $A$ 中 $S$ 对应的行线性无关（即 $rank (A_{S, *}) = r$ ），就有 $A_{S, S}$ 可逆。

如果没有斜对称的条件，那么只能知道存在 $| S | = | T | = r$ 使得 $A_{S, T}$ 可逆。

对可能不熟悉的读者：秩 $rank$ 就是 Gauss 消元后剩下的非零行数。 $n \times n$ 的矩阵 $A$ 可逆等价于 $rank A = n$ 。

首先，我们可以取一个足够大的素数 $p$ ，比如 $998244353$ ，做模数。然后我们随机选 $m$ 个数 $x_{1}, \dots, x_{m}$ 当作“未知数”。根据 Schwartz-Zippel 引理，只要我们接下来处理的都是多项式（比如，算矩阵的秩相当于让某个子矩阵的行列式非 $0$ ），就基本上不会算错。

接下来，我们算出前面这个矩阵 $Y = \sum_{i = 1}^{M} x_{i} (b_{i} \land c_{i})$ ，并计算它的秩。但是光算出这个还不行，我们还需要找到方案。

为了方便，我们可以求出 $Y$ 的最大行向量线性无关组 $S_{0}$ ，那么由于 $Y$ 斜对称，就知道 $Y_{S_{0}, S_{0}}$ 可逆，因此只需要关心这些行列即可。也就是说接下来我们把所有 $b_{i}, c_{i}$ 向量里其他分量都扔掉，只留下 $S_{0}$ 对应的这些维分量。这样的话 $Y$ 就可逆，我们可以算出他的逆 $N = Y^{- 1}$ 。

我们找到方案的办法非常简单粗暴：对于每一对 $(b_{i}, c_{i})$ ，我们把 $x_{i} (b_{i} \land c_{i})$ 从 $Y$ 里减去，看看它是否还可逆。如果可逆，就删掉这一对。否则留下他（把 $x_{i} (b_{i} \land c_{i})$ 再加回去）。

但是这样的话每次判断是否可逆的复杂度，即 Gauss 消元，是 $O (n^{3})$ ，这肯定不行。怎么办呢？

如果是一般的情况，其实有 Sherman-Morrison-Woodbury 定理，可以做到每次 $O (n^{2})$ ，总复杂度 $O (m n^{2})$ 。不过本题中这个复杂度还是不行。

5. 最后的算法

我们发现其实在本题我们有一个很好的性质，那就是 $b_{i}$ 和 $c_{i}$ 最多只有最多两个分量非 $0$ 。所以 $x_{i} (b_{i} \land c_{i})$ 只影响一个很小的子矩阵 $Y_{S, S}$ 。

因此我们可以使用 Harvey 定理——等一下， $Y$ 本来也不一定可逆欸。不过没关系，回归正题，由于 $x_{i} (b_{i} \land c_{i})$ 只影响一个很小的子矩阵，可以快速用 Harvey 定理判断更新后是否仍然满秩（可逆），如果可逆还可以快速求出更新后的逆。不过这样，复杂度也还是 $O (m n^{2})$ 。怎么进一步优化呢？

我们可以采用分治算法：如果我们有一堆边，他们对应的 $S$ （也就是他们连接的顶点对应的那些维）并起来是 $T$ ，而且这个 $T$ 不是很大的话，我们就可以对他们统一更新：我们发现 Harvey 定理里判断去掉这些边之后是否合法只需要用到 $N_{T, T}$ 。所以我们可以只更新 $N_{T, T}$ ，不更新 $N_{*, *}$ 。待到这些边都更新完之后，我们再去重新更新全部的 $N$ 。

如果把这个思路用到极致，就可以写出一个分治算法（伪代码）：

SOLVE(R, C):
  记 S = R 并 C
  // 执行时，要求 N_(S, S) 是正确的
  // 执行后，已经尝试删除了所有一端点属于 R 而另一端点属于 C 的边

  IF |R| > 2 or |C| > 2:
   将 R 划分为 R1, R2
   将 C 划分为 C1, C2
   FOR (i, j) in { (1, 1), (1, 2), (2, 1), (2, 2) }:
     SOLVE(Ri, Cj)
     使用 Harvery 定理更新 N_(S, S)
  ELSE IF R = { 2u-1, 2u }, C = { 2v-1, 2v } 且存在左端点为 u 右端点为 v 的边:
      使用 Harvery 定理判断是否可以删掉这条边
      如果是，删掉
      // 不需要更新 N，因为递归上去之后会更新

注意伪代码里默认了结点 $u$ 对应的行列还是 $2 u - 1, 2 u$ 。事实上由于我们删了一些行列，可能未必如此。所以分治需要注意不要把一个点对应的两个列分开。笔者的实现里是直接对结点分治，而非对矩阵行列分治。

这样的话，复杂度为 $T (n) = 4 T (\frac{n}{2}) + O (n^{w})$ 。根据主定理可以计算出 $T (n) = O (n^{w})$ 。如果 $w = 2$ ，那么 $T (n) = O (n^{2} \log n)$ 。

代码

由于这个矩阵是 $2 n \times 2 n$ （如果大部分点都没有颜色并且答案几乎是满的），所以常数有点大。

笔者和出题人码风不通，没看懂他的实现。笔者本人的实现最后使用了 $O (n^{3})$ 算法，不过模数取用了 $10^{8} + 7$ ，这样可以在矩阵乘法时全部乘完加完最后取模，降低了很多常数。

#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>

typedef long long LL;
const int mod = 100000007;
const int N = 303;
const int NN = 606;
const int M = N * N / 2;

LL pow_mod(LL a, LL b) {
  LL ans = 1;
  for (; b; b >>= 1, a = a * a % mod)
    if (b & 1) ans = ans * a % mod;
  return ans;
}

inline void mul(int m, LL *A, LL x) {
  for (int i = 0; i < m; ++i) A[i] = A[i] * x % mod;
}

inline void add(int m, LL *A, const LL *B, LL x) {
  for (int i = 0; i < m; ++i) A[i] = (A[i] + x * B[i]) % mod;
}

inline void swap(int m, LL *A, LL *B) {
  for (int i = 0; i < m; ++i) std::swap(A[i], B[i]);
}

void rank(int n, int m, LL A[NN][NN], int &r, int *S) {
  // A[S][*] has rank r.
  static int id[NN];
  for (int i = 0; i < n; ++i) id[i] = i;
  for (int i = 0, j = 0; j < m; ++j) {
    int k = i;
    while (k < n && A[k][j] == 0) ++k;
    if (k == n) continue;
    if (i != k) {
      std::swap(id[i], id[k]);
      swap(m, A[i], A[k]);
    }
    S[r++] = id[i];
    LL t = pow_mod(A[i][j], mod - 2);
    for (int k = i + 1; k < n; ++k)
      add(m, A[k], A[i], -t * A[k][j] % mod);
    ++i;
  }
}

bool inv(int n, const LL A_[NN][NN], LL B[NN][NN]) {
  static LL A[NN][NN];
  // return true iff A is invertible.
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      A[i][j] = A_[i][j], B[i][j] = i == j;
  for (int i = 0; i < n; ++i) {
    int k = i;
    while (k < n && A[k][i] == 0) ++k;
    if (k == n) return false;
    if (i != k) {
      swap(n, A[i], A[k]);
      swap(n, B[i], B[k]);
    }
    LL t = pow_mod(A[i][i], mod - 2);
    mul(n, A[i], t); mul(n, B[i], t);
    for (int k = 0; k < n; ++k) if (k != i) {
      add(n, B[k], B[i], -A[k][i] % mod);
      add(n, A[k], A[i], -A[k][i] % mod);
    }
  }
  return true;
}

void mul(int n, int m, int r, const LL A[NN][NN], const LL B[NN][NN], LL C[NN][NN]) {
  static LL tmp[NN][NN];
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < r; ++j) tmp[i][j] = 0;
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < m; ++j)
      for (int k = 0; k < r; ++k)
        tmp[i][k] = tmp[i][k] + A[i][j] * B[j][k];
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < r; ++j) C[i][j] = tmp[i][j] % mod;
}

#define take(A, X, Y, FORMULA)                     \
  for (int i = 0; i < X; ++i)                   \
    for (int j = 0; j < Y; ++j)                 \
      A[i][j] = FORMULA

bool modify(const LL A[NN][NN], LL B[NN][NN], const LL A1[NN][NN],
            int s, const int *S, int t, const int *T) {
  // assuming B = A^(-1) and A[i][j] != A1 only when i, j in S
  // if A1 is singular, return false
  // else:
  // modify B such that B_(T, T) = A1^(-1)_(T, T)
  // only use A[S][S], A1[S][S], B[T][T].
  static LL D[NN][NN], C[NN][NN], P[NN][NN];
  take(D, s, s, A1[S[i]][S[j]] - A[S[i]][S[j]]);
  take(P, s, s, B[S[i]][S[j]]);
  mul(s, s, s, D, P, P);
  for (int i = 0; i < s; ++i) (++P[i][i]) %= mod;
  if (!inv(s, P, C)) return false;
  mul(s, s, s, C, D, C);
  take(P, t, s, B[T[i]][S[j]]);
  mul(t, s, s, P, C, C);
  take(P, s, t, B[S[i]][T[j]]);
  mul(t, s, t, C, P, C);
  for (int i = 0; i < t; ++i)
    for (int j = 0; j < t; ++j) {
      int ti = T[i], tj = T[j];
      B[ti][tj] = (B[ti][tj] - C[i][j]) % mod;
    }
  return true;
}

int C[N], L[N];
LL X[N][N];
bool link[N][N];
int n, nn;
int id[NN];
int t, T[NN];

void Add(LL A[NN][NN], int u, int v, int i, int j, LL p) {
  int x = L[u] + i, y = L[v] + j;
  if (C[u] && i) x = L[u], p = p * C[u] % mod;
  if (C[v] && j) y = L[v], p = p * C[v] % mod;
  if (id[x] >= 0 && id[y] >= 0) A[id[x]][id[y]] = (A[id[x]][id[y]] + p) % mod;
}

inline void Add(LL A[NN][NN], int u, int v, LL p) {
  // b = e(u0) - e(v0)
  // c = e(u1) - e(v1)
  // A += p (b wedge c)
  Add(A, u, u, 0, 1, p);
  Add(A, u, v, 0, 1, -p);
  Add(A, v, u, 0, 1, -p);
  Add(A, v, v, 0, 1, p);
  Add(A, u, u, 1, 0, -p);
  Add(A, u, v, 1, 0, p);
  Add(A, v, u, 1, 0, p);
  Add(A, v, v, 1, 0, -p);
}

LL Y[NN][NN], Z[NN][NN];

int ss[10], SS[10][NN];
LL AA[10][NN][NN], BB[10][NN][NN];

bool solve(int l1, int r1, int l2, int r2, int dep) {
#define copyS(U, V)                             \
  for (int i = 0; i < s; ++i)                   \
    for (int j = 0; j < s; ++j)                 \
      U[S[i]][S[j]] = V[S[i]][S[j]]
  if (r1 < l1 || r2 < l2) return false;
  int &s = ss[dep] = 0;
  int *S = SS[dep];
  LL (*A)[NN] = AA[dep], (*B)[NN] = BB[dep];
  for (int i = L[l1]; i < L[r1 + 1]; ++i) if (id[i] >= 0) S[s++] = id[i];
  for (int i = L[l2]; i < L[r2 + 1]; ++i) if (id[i] >= 0) S[s++] = id[i];
  std::sort(S, S + s); s = std::unique(S, S + s) - S;
  copyS(A, Y);
  copyS(B, Z);

  if (r1 == l1 && r2 == l2) {
    if (!link[l1][l2]) return false;
    Add(A, l1, l2, -X[l1][l2]);
    if (!modify(Y, B, A, s, S, 0, NULL)) return false;
    link[l1][l2] = false;
    Add(Y, l1, l2, -X[l1][l2]);
    return true;
  }

  int m1 = (l1 + r1) / 2, m2 = (l2 + r2) / 2;
  bool changed = false;
  for (int a = 0; a < 2; ++a)
    for (int b = 0; b < 2; ++b) {
      int u1 = l1, v1 = m1, u2 = l2, v2 = m2;
      if (a) u1 = m1 + 1, v1 = r1;
      if (b) u2 = m2 + 1, v2 = r2;
      if (solve(u1, v1, u2, v2, dep + 1)) {
        changed = true;
        // fprintf(stderr, "%d %d\n", ss[dep + 1], s);
        if (ss[dep + 1] != s) {
          if (false) {
            inv(s, Y, Z);
          } else {
            copyS(Z, B);
            modify(A, Z, Y, ss[dep + 1], SS[dep + 1], s, S);
          }
        }
        copyS(A, Y);
        copyS(B, Z);
      }
    }
  return changed;
}

bool vis[N];
int path[N][N], len[N], ans;

bool dfs(int x, bool rt) {
  vis[x] = true;
  if (C[x] && !rt) {
    ++ans;
    path[ans][len[ans]++] = x;
    return true;
  }
  for (int y = 1; y <= n; ++y)
    if ((link[x][y] || link[y][x]) && !vis[y]) {
      int t = dfs(y, false);
      if (t) {
        path[ans][len[ans]++] = x;
        return true;
      }
    }
  return false;
}

int main() {
  int m;
  scanf("%d%d", &n, &m);
  for (int i = 1; i <= n; ++i) {
    scanf("%d", &C[i]);
    L[i] = nn++;
    if (!C[i]) nn++;
  }
  L[n + 1] = nn;
  for (int i = 0; i < nn; ++i) id[i] = i;
  for (int i = 0; i < m; ++i) {
    int u, v;
    scanf("%d%d", &u, &v);
    if (u == v || link[u][v] || link[v][u]) continue;
    link[u][v] = true;
    X[u][v] = rand() % mod; // TODO: ERROR
    Add(Y, u, v, X[u][v]);
  }
  rank(nn, nn, Y, t, T);
  memset(id, -1, sizeof id);
  for (int i = 0; i < t; ++i) id[T[i]] = i;
  memset(Y, 0, sizeof Y);
  for (int u = 1; u <= n; ++u)
    for (int v = 1; v <= n; ++v)
      if (link[u][v]) Add(Y, u, v, X[u][v]);
  inv(nn, Y, Z);
  solve(1, n, 1, n, 0);
  for (int i = 1; i <= n; ++i)
    if (!vis[i] && C[i]) dfs(i, true);
  printf("%d\n", ans);
  for (int i = 1; i <= ans; ++i) {
    printf("%d", len[i]);
    for (int j = 0; j < len[i]; ++j)
      printf(" %d", path[i][j]);
    puts("");
  }
  return 0;
}

_rqy的博客

博客