AtCoder Regular Contest 163 C,D問題メモ

AtCoder Regular Contest 163

C - Harmonic Mean

• 以下の条件を全て満たす数列 $A=(A_1,A_2,...,A_N)$ が存在するか判定し、存在する場合は一例を示せ
  • 項数は $N$
  • 逆数の総和が1（$\displaystyle \sum_{i=1}^{N}\frac{1}{A_i}=1$）
  • $A_i$ は全て相異なる
  • $A_i \le 10^9$
• 1つの入力につき $T$ ケース与えられる
• $1 \le T \le 500$
• $1 \le N \le 500$

$N=5$ のサンプルの答え 3 4 5 6 20 をOEISで検索してみると、 A216975 が見つかり、「エジプト式分数」というキーワードを得る。

• エジプト式分数
• エジプト分数（単位分数の和）に関する4つの話題 | 高校数学の美しい物語

作り方に決まったものはないが、単位分数1つを2つに置き換える簡潔な方法として、以下がある。

• $\dfrac{1}{a}=\dfrac{1}{a+1}+\dfrac{1}{a(a+1)}$

$N=1$ の時は $A=(1,)$、$N=2$ の時は解無しで、 $N=3$ 以降は（かぶりと、上限に気をつければ）再帰的に置き換えていくことで必ず構築できそう。

$N=3$ の時の解 $A=(2,3,6)$ からはじめて、被らないように変換していきたい。

小さい方から変換していくと、かぶりの管理が難しそう。
（ある時点では $a+1$ が $A$ に存在するので $a$ の変換ができなかったが、 $a+1$ が変換されて無くなったので $a$ も変換できるようになった、という差し戻しが連鎖的に発生する）

そこで、大きい方から変換していくことにする。

• 2つのリスト $L=[2,3,6],Fixed=[]$ を用意する
• 以下を、$|L|+|Fixed|=N$ になるまで繰り返す
  • $L$ から最大の数 $x$ を取り出す
  • $x+1$ や $x(x+1)$ が、既にFixedに存在したり、$10^9$ を超える場合、$x$ をFixedに追加する
  • 存在しない場合、$L$ に $x+1$ と $x(x+1)$ を追加する

この方法だと一度変換に失敗した数はずっと失敗する（差し戻しが発生しない）ので、 次に何に対して変換を行えばよいかの管理が楽。

$N=500$ で試してみて問題なかったので、制約の範囲内ではとりあえず大丈夫。
（$N=31623$ まではできたので、結構余裕がある）

Python3

t = int(input())
LIMIT = 1_000_000_000
for _ in range(t):
    n = int(input())
    if n == 1:
        print('Yes')
        print(1)
        continue
    if n == 2:
        print('No')
        continue
    ans = [2, 3, 6]
    fixed = set()
    l = 3
    while l < n:
        p = ans.pop()
        q = p + 1
        r = p * q
        if q in fixed or r in fixed:
            fixed.add(p)
            continue
        if r > LIMIT:
            fixed.add(p)
            continue
        ans.append(q)
        ans.append(r)
        l += 1

    ans += sorted(fixed)
    print('Yes')
    print(*ans)

D - Sum of SCC

• $N$ 頂点のトーナメントグラフのうち、含まれる昇辺がちょうど $M$ 本であるものを「良いグラフ」とする
  • トーナメント：無向完全グラフの $\dfrac{N(N-1)}{2}$ 本の辺に、それぞれ向きを付けて有向グラフとしたもの
  • 昇辺（勝手に命名）：頂点 $u→v$ に辺があるとき、頂点番号が $u \lt v$ であるもの
• 全ての良いグラフについて、その強連結成分の個数の総和を、$\mod{998244353}$ で求めよ
• $1 \le N \le 30$

公式Editorialの解法。

トーナメントでは、全ての頂点間に有向辺があるので、強連結成分を縮約した上でのトポロジカルソートを考えると、その順番は一意に決まる。
自分より下流にある頂点への辺は、全て 上流→下流 の向きになっている。

  ②  →  ①  →  ③
 ④⑤ --------→  ⑥

逆に、「2つの頂点集合 $A,B$ があって、$A$ の頂点と $B$ の頂点を結ぶ全ての辺は $A→B$ の向きに張られている」場合、 $B$ から $A$ に行くことは不可能であり、強連結が途切れているといえる。

このような「強連結成分の境界」に着目する。
境界の個数を数えると、それに良いグラフの個数 $\binom{N(N-1)/2}{M}$ を加えたものが答えとなる。

主客転倒の考え方を用いると、以下のようになる。

• 頂点を空でない2グループ $A,B$ に分け、全ての辺が $A→B$ の向きに張られるような「分け方×辺の張り方」のパターン数を数える
• そのうち昇辺の個数が $M$ であるものの総和が、全ての良いグラフに対する強連結成分の境界の個数である

主客転倒すると、全ての辺が $A→B$ の向きである時点で $A,B$ 間に確実に境界が1つ存在し、 他に $A$ や $B$ が複数の強連結成分に分けられようが関係ないところが計算を楽にする。

頂点番号の小さい順に $A,B$ どちらに入れるかを確定させていくDPで求められる。
具体的な遷移や実装は公式Editorial参照。

ABC-F～G などに出てくるような解法を繰り返す・組み合わせることでも、地道にやっていけば解けなくはない。

以下のDPをおこないたい。

• $DP[i,j,k] = $ トポロジカルソートの上流側から、頂点を $i$ 個確定し、昇辺が $j$ 本の場合の（k=0:パターン数, k=1:のべ連結成分数）

$j$ で数える昇辺の範囲は、確定した $i$ 個に少なくとも一方の端点を持つ辺とする。 （未処理の $N-i$ 個間に張られる辺は対象外ということ）

                     未処理頂点
○→○ → ○←○ → （○○○）
↑×↓ →  ↘↗  →
○←○ →   ○   →

|--- i個 j本-------|

このような状態から、次に来る強連結成分を未処理頂点から切り出すことで遷移していく。
別途、以下の情報があればよい。

• $DP2[p,q,r] = $ 残り頂点 $p$ 個中から $q$ 個選び、新規の昇辺が $r$ 本できるパターン数
  • 必ず $q$ 個の頂点全てが強連結になるように辺の向きを決める
  • $r$ 本を数える範囲は、$q$ 個同士を結ぶ辺 ＋ $q$ 個から未処理頂点への辺

残り頂点の番号が具体的に分からなくても、残り頂点内での辺が昇辺かどうかは相対的な順序で決まるので、$p,q$ さえ同じなら同じになるところがポイント。

もらうDPで考えて、$DP[i-q,j-r]$ から新規で $q$ 頂点、$r$ 昇辺の強連結成分を1つ決めて $DP[i,j]$ にする方法は、

• $DP[i,j,0]+=DP[i-q,j-r,0] \times DP2[N-i,q,r]$
• $DP[i,j,1]+=(DP[i-q,j-r,0]+DP[i-q,j-r,1]) \times DP2[N-i,q,r]$

このようになる。

普通に $i,j,q,r$ でループを回すと、$i,q$ は頂点数 $O(N)$、$j,r$ は辺数 $O(N^2)$ なので $O(N^6)$ となってしまうが、 もらうDPの形にすると $j,r$ を固定する部分は畳み込みができるので、$O(N^4 \log{N})$ になる。

さて、では解けたかというとそうでもなく、$DP2$ を求めるのが難しい。

$DP2$ を、以下の2つに分けて求める。

• $DP3[q,r]=q$ 頂点のトーナメントのうち、全体が強連結で、昇辺が $r$ 本のものの個数
• $DP4[p,q,s]=$ 残り頂点 $p$ 個中から $q$ 個選び、「$q$ 個と未処理頂点間で」新規の昇辺が $s$ 本できるパターン数

この2つがあれば、$DP2$ は $p,q$ ごとに $r,s$ を添字とした畳み込みにより $O(N^4 \log{N})$ で求まる。

DP4は、“スマートな解法”と似た感じで、小さい頂点から「$q$ 個に含まれる側と、未処理側のどちらに入れるか」を確定させるDPでできる。$O(N^4)$

DP3は、なんかOEISをあさると A339590 が見つかる。
さらにそこからこんな 論文 が見つかる。

それによると、DP3は以下の要領で計算できる。

$u$ を形式的な変数とした形式的冪級数 $t_n(u)$ を、以下で定義する。$\binom{n}{k}_q$ はq-二項係数とする。

• $\displaystyle t_n(u)=(1+u)^{\binom{n}{2}}-\sum_{k=1}^{n-1} \binom{n}{k}_u(1+u)^{\binom{n-k}{2}}t_k(u)$

• メモ: q-二項係数 - noshi91のメモ

$t_1(u)=1,t_2(u)=0$ として、$n$ の小さい方から計算できる。すると、$[u^r]t_q(u)$（※ $t_q(u)$ における $u^r$ の係数）が、$DP3[q,r]$ に相当する。

これで全ての情報が計算できて、解けました。めでたしめでたし。

• $\displaystyle t_n(u)=(1+u)^{\binom{n}{2}}-\sum_{k=1}^{n-1} \binom{n}{k}_u(1+u)^{\binom{n-k}{2}}t_k(u)$

$(1+u)^{\binom{n}{2}}$ は普通の二項係数を作る式なので、$u^0,u^1,u^2,...$ の係数は「$n$ 頂点のトーナメントで昇辺が $0,1,2,...$ 個のものの個数」を表す。

そこから、内部を2つ以上の強連結に分けられるものを引いている。

重複しないよう、一番上流の強連結成分を $k$ 頂点と固定して、残りを自由に決めている感じ。

• $t_k(u)$：一番上流の $k$ 頂点を強連結成分にするパターン数
• $(1+u)^{\binom{n-k}{2}}$：残り $n-k$ 頂点同士の辺の決め方

なので、Σの中でのあと1項 $\binom{n}{k}_u$ が、「$n$ 頂点から $k$ 個選ぶパターン数、ただし選んだ頂点→選ばれなかった頂点間にできる昇辺の個数毎」を示していることになる。

これって、$DP4$ と一緒じゃん！ 知らないうちに同じことをしてしまっていたらしい。

noshi氏のサイトによると、q-二項係数は以下の性質があるらしい。

• 0を $a$ 個、1を $b$ 個並べた数列であって転倒数が $k$ のものの個数は $[q^k]\binom{a+b}{a}_q$

「昇辺の個数ごとに $n$ 頂点から $k$ 個選ぶパターン数」も、この問題に言い換えられることが確認できる。

Python3

import numpy as np

n, m = map(int, input().split())
MOD = 998244353


# https://arxiv.org/pdf/1909.01550.pdf

def convolve(f, g):
    fft_len = 1
    true_len = len(f) + len(g) - 1
    while fft_len < true_len:
        fft_len <<= 1

    Ff = np.fft.rfft(f, fft_len)
    Fg = np.fft.rfft(g, fft_len)
    Fh = Ff * Fg

    h = np.fft.irfft(Fh, fft_len)
    h = np.rint(h).astype(np.int64)

    return h[:true_len]


def convolve_mod(f, g, p):
    f1, f2 = np.divmod(f, 1 << 15)
    g1, g2 = np.divmod(g, 1 << 15)

    a = convolve(f1, g1) % p
    c = convolve(f2, g2) % p
    b = (convolve(f1 + f2, g1 + g2) - (a + c)) % p

    h = (a << 30) + (b << 15) + c
    return h % p


def precompute_factorials(n, MOD):
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f = f * m % MOD
        factorials.append(f)
    f = pow(f, MOD - 2, MOD)
    finvs = [1] * (n + 1)
    finvs[n] = f
    for m in range(n, 1, -1):
        f = f * m % MOD
        finvs[m - 1] = f
    return factorials, finvs


facts, finvs = precompute_factorials(500, MOD)

qb = [[[1]]]
for l in range(1, n + 1):
    row = [[1]]
    for k in range(1, l):
        a = qb[l - 1][k]
        b = qb[l - 1][k - 1]
        p = len(a)
        q = len(b)
        r = max(p + k, q)
        res = b + [0] * (r - q)
        for i in range(p):
            res[i + k] += a[i]
        row.append(res)
    row.append([1])
    qb.append(row)

ttt = [
    [0],
    [1],
    [0],
    [0, 1, 1, 0],
]

for l in range(4, n + 1):
    ll = l * (l - 1) // 2
    f = facts[ll]
    result = np.array([f * finvs[i] * finvs[ll - i] % MOD for i in range(ll + 1)], dtype=np.int64)
    for k in range(1, l):
        if k == 2:
            continue
        a = qb[l][k]
        ll2 = (l - k) * (l - k - 1) // 2
        f2 = facts[ll2]
        b = [f2 * finvs[i] * finvs[ll2 - i] % MOD for i in range(ll2 + 1)]
        c = ttt[k]
        res = convolve_mod(a, b, MOD)
        res = convolve_mod(res, c, MOD)
        result -= res
    result %= MOD
    ttt.append(result)

dp = [[[0] * (m + 1) for _ in range(n + 1)] for _ in range(n + 1)]
dp[0][0][0] = 1

for i in range(n):
    for j in range(i + 1):
        for k in range(m + 1):
            dp[i + 1][j + 1][k] += dp[i][j][k]
            if k + j <= m:
                dp[i + 1][j][k + j] += dp[i][j][k]
dp = np.array(dp, dtype=np.int64)

dp2 = np.zeros((n + 1, m + 1), np.int64)
dp3 = np.zeros((n + 1, m + 1), np.int64)
dp2[0, 0] = 1
dp3[0, 0] = 0
for i in range(1, n + 1):
    for j in range(1, i + 1):
        tmp = convolve_mod(dp[n - (i - j)][j], ttt[j], MOD)[:m + 1]
        tmp2 = convolve_mod(tmp, dp2[i - j], MOD)[:m + 1]
        tmp3 = convolve_mod(tmp, dp3[i - j], MOD)[:m + 1]
        dp2[i] += tmp2
        dp3[i] += tmp3 + tmp2
    dp2[i] %= MOD
    dp3[i] %= MOD

print(dp3[n, m])