yukicoder contest 420 開催記

前回のコンテスト で重くし過ぎたことを反省してやや軽めにしようと思っていたんですが、結局後ろが大変になってしまいました。

言い訳をさせて頂くと、テスター作業を依頼する前に設定していた難易度は☆1.5 - 2.0 - 2.5 - 3.0 - 3.0 - 3.5 - 3.5 - 4.0でした。

A - Zero-Sum Submatrices

ギャグその①。

2022年に作っていた問題です。前回のF はこの問題からの派生です。

想定解の構築を直接思いつけなくても、解説の後半で説明している再帰的な構築はよく見る気がします。

B - Prime Sum

ギャグその②。

一番最初に考えた設定では構築ありで $X _ i \geq 1$ でした。 この場合は全ての $i$ に対して $X _ i = 1$ とすれば条件を満たすんですが、流石に度が過ぎるかと思い $X _ i \geq 2$ にしました。

この問題のように、与えられたペアを表す辺を張ってできるグラフで考えると上手くいく問題というのはよく見ますが、私は初見で結構感動した覚えがあります。cp-unspoiler とかは結構お気に入りの問題で、競プロの説明に使ったこともあります。

C - Minimize Inversions of Deque

前回の B に引き続いて転倒数に関する問題。

操作を後ろから考えると良い性質が見える問題というのもよくあって、例えば 前回の C もそうでした。

ABCのE, F辺りに置かれそうという印象です。

D - Decreasing Modulo Nim

結論がかなり面白くてお気に入りです。

一番最初に考えた設定は $m \to \infty$ としたものでした。 ただ、これだと通常の Nim の勝敗判定と全く同じになってしまったので、もう一捻り加えました。

明らかに $x$ に関する制約が厄介なので、初手で先手が $x$ を $0$ にしてしまうことを考えようというのが解説の 1 行目の気持ちです。そこから先の考察は比較的素直だと思います。

$x$ という global な (?) 値が存在するので、山毎の Grundy 数を求めて~、という解法は棄却する必要があります。

E - Constrained Permutation

区間スケジューリングを題材に作りました。この問題も性質が綺麗でお気に入りです。

以下のコードに示した貪欲法で $k$ が条件を満たすかを判定できます。

# ranges: 閉区間のリスト
def check(k, ranges):
    n = len(ranges)
    # 左端の昇順にソート
    ranges.sort_by_left_bound()
    # 右端の最小値を管理する優先度付きキュー
    pq = MinPriorityQueue()
    j = 0
    for v = k+1, k+2, ..., k+n:
        while j < n and ranges[j].left <= v:
            pq.push(ranges[j].right)
            j += 1
        # 条件1
        if pq.empty():
            return False
        # 条件2
        if pq.pop().right < v:
            return False
    return True

上記のコードにおける条件1は左端の(多重)集合にのみ依存しており、この部分で False が返らないような条件を考えているのが解説の初手に対応します。 条件2に関しては、直感的には $k$ が小さい方が有利そうであり、実際にその通りであることを解説で証明しています。

証明がやや難しいので、未証明で通した人が多そうです。

F - Trees on Graph Paper

出題方法にかなり悩んだ問題です。

初めは次数 1 の頂点がちょうど $K$ 個の良い全域木を数えたかったのですが、私が解けなかったので次数 1 の頂点の $x$ 座標の積の和という奇妙な設定になりました。DPの遷移が $x$ 座標に依存するようにしたのは Berlekamp–Massey 防止策です。

ナイーブな DP を考えると考えるべき情報が多く遷移を詰めるのが大変ですが、考察を進めると僅か 3 状態のシンプルなDPになります。

解説に余談として書いた良い全域木の個数とフィボナッチ数の関係が面白いです。OEISで解けるので出題は出来ませんが...

G - Generalized Hitting Set

解説に気合が入っているのでぜひ最後まで読んでほしい。

包除のような謎の係数 (この係数の正体は解説②で説明されます) が現れるのが面白いと思っています。 「上位集合系のゼータ変換 → 各点積 → 下位集合系のゼータ変換」というアルゴリズムも個人的には違和感があって面白いと思っています。解説の最後に変種として紹介している問題の解法は更に不思議な感じになっています。

$n = 24$ という過激な制約は、各 $i\in\lbrace 0,1,\ldots,n\rbrace,\ T\subseteq\lbrace 1,2,\ldots,n\rbrace$ に対して $\lvert S _ j \cap T\rvert = i$ を満たす $j$ の個数を数える DP 解 (空間 $\Theta(n 2 ^ n)$ / 時間 $\Theta(n ^ 2 2 ^ n)$) を TL,ML の両方で落とそうとしたためです。

ところで、問題準備中に PyPy に関して謎の現象に遭遇しました。 以下の2つのコードの違いは 40 行目あたりの del b の有無のみなんですが、実行時間が 800 ms 程度も変わっています。 原因が分かる人がいたら教えてください。

H - Sum of Products of Interval Length

唯一解法から作った問題です。

いわゆる(?)同じものの連続に関する包除ですが、その勉強がてら作りました。

積の和典型で上手くDPすれば分割統治FFTで $\Theta(n (\log n) ^ 2)$ 時間で解けるようです。

階乗 mod 素数

Many Factorials を準備したので宣伝を兼ねて。

$N! \bmod P$ ($P$ は固定された素数) をたくさんの非負整数 $N$ に対して評価する状況を考える。

$N$ の上限が小さい (例えば $N \leq 10 ^ 7$ が保証される) 場合

$N$ の上限を $M $ とする。

予め $N! \bmod P$ を全ての $N=0,1,\ldots,M $ に対して前計算しておくことで、各クエリを定数時間で処理できる。

$N$ の上限が大きい場合

本記事の本題。

方法① 埋め込み

$N \geq P$ ならば $N! \equiv 0 \bmod P$ なので、$N\lt P$ を仮定する。

適当なブロックサイズ $B$ を決めて $(iB)! \bmod P\ (i=0,1,\ldots,\lfloor (P-1) / B\rfloor$ を時間を掛けて計算し、その結果をソースコードに直接埋め込む。各クエリは前計算の結果として得られている $(\lfloor N/B\rfloor B)!$ からの差分だけを計算することで $O(B)$ 時間で処理できる。

なお、定数倍改善として、$\lfloor N/B\rfloor B$ と $(\lfloor N/B\rfloor + 1) B$ のうちより $N$ に近い方からの差分を計算するという方法が考えられる。$(\lfloor N/B\rfloor + 1) B \geq P$ となるケースに注意する必要があるが、この場合は $(P - 1)! \equiv -1 \pmod{P}$ (ウィルソンの定理 - Wikipedia) を用いて $(P - 1)!$ からの差分を計算するとよい。

この方法は手軽だが、次の点に注意する必要がある。

  • $P$ が実行時に与えられる場合は使えない
  • ソースコード長の制限により $B$ を小さく取れないことがある

方法② $(iB)! \bmod P$ の実行時計算

本節の内容は(現在は削除されている)min_25さんの記事で紹介されていた手法を元にしている。

$(iB)! \bmod P$ を実行時に計算することで方法①の弱点を克服する。

以下では Shift of Sampling Points of Polynomial (以下評価点シフトと呼ぶ) を道具として認める。詳細については Shift of Sampling Points of Polynomial | cp-library-cpp などを参照。

簡単のため $B = 2 ^ K$ とする。また $\displaystyle f _ i(x) \coloneqq \prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i x + j) = (2 ^ i x + 1) (2 ^ i x + 2) \cdots (2 ^ i x + 2 ^ i - 1)$ と定める。

$f _ K(0), f _ K(1), \ldots, f _ K(\lfloor P/2 ^ K\rfloor)$ を計算することが目標である。

さて、$f _ i$ と $f _ {i + 1}$ の間には次の関係が成り立つ。

$$\begin{aligned} f _ {i + 1} (x) &= \prod _ {j = 1} ^ {2 ^ {i + 1} - 1} (2 ^ {i + 1} x + j) \\ &= \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ {i + 1} x + j)\right)\cdot (2 ^ {i + 1} x + 2 ^ i) \cdot \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ {i + 1} x + 2 ^ i + j)\right) \\ &= \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i (2x) + j)\right)\cdot 2 ^ i (2x + 1) \cdot \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i (2x + 1) + j)\right) \\ &= f _ i(2x) \cdot f _ i(2x + 1) \cdot 2 ^ i (2x + 1). \end{aligned}$$

つまり $f _ i(0), \ldots, f _ i(2 ^ i - 1)$ が既知であると仮定すれば、次の手続きにより $f _ {i + 1}(0), \ldots, f _ {i + 1}(2 ^ {i + 1} - 1)$ を計算することができる。

  1. 評価点シフトを用いて $f _ i(2 ^ i), f _ i(2 ^ i + 1), \ldots, f _ i(2\cdot 2 ^ {i + 1} - 1)$ を計算する。
    • 評価点シフトを用いるための前提条件として、$f _ i(0), f _ i(1), \ldots, f _ i(\deg f _ i)$ が既知でなければならないことに注意する。今回は $\deg f _ i = 2 ^ i - 1$ より条件を満たしていることを確認できる。
  2. 各 $j = 0, 1, \ldots, 2 ^ {i + 1} - 1$ に対する $f _ {i + 1}(j)$ を、$f _ {i + 1}(j) = f _ i(2j) \cdot f _ i(2j + 1) \cdot 2 ^ i (2j + 1)$ により得る。

1 は $O(i 2 ^ i)$ 時間で、2 は $O(2 ^ i)$ 時間で可能である。結局 $f _ 0(0) = 1$ から初めて上記の手続きを繰り返すことで $O(\sum _ {i = 0} ^ {K - 1} i 2 ^ i) = O(K 2 ^ K)$ 時間で $f _ K(0), \ldots, f _ K(2 ^ K - 1)$ を得ることができる。

$2 ^ K \leq \lfloor P/2 ^ K\rfloor$ の場合、最後に $f _ K(0), \ldots, f _ K(2 ^ K - 1)$ に対してもう一度評価点シフトを用いることで $f _ K(2 ^ K), f _ K(2 ^ K + 1), \ldots, f _ K(\lfloor P/2 ^ K\rfloor)$ を得ることができる。これは $O((P / 2 ^ K) \log P)$ 時間で可能である。

結局、$2 ^ K$ を $\sqrt{P}$ 付近に取れば、前計算は $O(\sqrt{P} \log P)$ 時間で行うことができる。

実用においては、$K$ をより小さく取って前計算に時間を掛け、$1$ 回のクエリ処理に掛かる時間を短くするとよい (実験結果 を参照)。

方法③ $(iB)! \bmod P$ の実行時計算 + オフライン処理による高速化 (1)

クエリをオフラインで処理してよいと仮定する。

$(iB)! \bmod P$ を前計算しておけば、各クエリは $\displaystyle (\lfloor N/B\rfloor B)! \cdot N(N-1)\cdots(N-(N\bmod B)+1)$ として計算できる。$(\lfloor N/B\rfloor B)!$ は前計算の結果から既知なので、各 $N$ に対する $N(N-1)\cdots(N-(N\bmod B)+1)$ を高速に計算したい。

つまり、$\displaystyle f(a,b)\coloneqq \prod _ {i = 0} ^ {b - 1} (a - i)$ と定義すれば、$f(N, N \bmod B)$ が計算したいものである。

各クエリで $N \bmod B$ が全て等しい場合は比較的容易であり、その値を $R$ として、$\displaystyle \prod _ {i = 0} ^ {R - 1} (x - i)$ を $R$ 点ずつに分けて多点評価すると $O( (R + Q) (\log R) ^ 2)$ 時間となる ($Q $ はクエリの数)。

一般にはこのような性質は成り立たない。しかし、$0\leq i\leq b$ なる整数 $i$ に対して $f(a, b) = f(a, i) \cdot f(a - i, b - i)$ が成り立つことを用いて $N \bmod B$ を $2$ 冪の和に分解することで、$b$ の種類数を高々 $\lfloor \log _ 2 (B-1) \rfloor$ に減らすことができる。

具体的には、例えば $N = 49,\ N \bmod B = 13 = 2 ^ 0 + 2 ^ 2 + 2 ^ 3$ のとき $f(49,13) = f(49, 2 ^ 0) \cdot f(49 - 2 ^ 0, 2 ^ 2) \cdot f(49 - 2 ^ 0 - 2 ^ 2, 2 ^ 3)$ のように分解する。

つまり各 $t = 0, 1, \ldots, \lfloor \log _ 2 (B-1) \rfloor$ に対して $\displaystyle \prod _ {i = 0} ^ {2 ^ t - 1} (x - i)$ を高々 $Q$ 個の点で評価すればよく、計算量は $\displaystyle O\left( \sum _ {i = 0} ^ {\lfloor \log _ 2 (B-1)\rfloor - 1} (2 ^ i + Q) i ^ 2 \right) = O(Q (\log B) ^ 3 + B (\log B) ^ 2)$ となる。

$B$ を $\sqrt{P}$ 付近に取ると、前計算と併せて全体で $O(Q (\log P) ^ 3 + \sqrt{P} (\log P) ^ 2)$ 時間となる。

この方法でも、実用上は $B$ をより小さく取って多点評価の回数を減らした方が速い (実験結果 を参照)。

方法④ $(iB)! \bmod P$ の実行時計算 + オフライン処理による高速化 (2)

方法③では多点評価を $O(\log B)$ 回に分けて行ったが、まとめて $O((Q + \sqrt{P}) (\log P) ^ 2)$ 時間で行うことができる (参考: yukicoder Paint and Fill 解法詳細解説 | Mathenachia)。

この方法では $B$ の大きさによる実行時間への影響は小さい (実験結果 を参照)。

実験

$P = 998244353,\ Q = 5\times 10 ^ 5$ のランダムケースで $K = \log _ 2 B$ を動かして実行時間 (秒) を計測した。なお、③' は除算を用いない高速な多点評価を用いて③を高速化したものである。コードの詳細については factorial mod prime · GitHub を参照。

方法\$K$ $7$ $8$ $9$ $10$ $11$ $12$ $13$ $14$ $15$
0.80 0.53 0.46 0.56 0.90 1.6 3.0 6.0 12
1.0 1.1 1.3 1.7 2.2 2.7 3.4 4.1 4.9
③' 0.90 0.78 0.83 1.0 1.3 1.6 2.0 2.5 3.4
1.6 1.4 1.3 1.3 1.3 1.4 1.6 1.7 2.1

$P$ が小さいため方法②が一番速いという結果になった。また、③や④はクエリをオンラインで処理できない一方で②はオンラインで処理ができるため、競技の文脈では (大抵の場合 $P \lt 2 ^ {31}$ なので) 方法②を使用するのが無難だと思われる。

実装

最も高速であった方法②で $K=9$ としたものの実装を以下に示す。

judge.yosupo.jp

AtCoder Regular Contest 136 F - Flip Cells

未証明要素を多分に含む怪しい解法で通してしまった

問題

atcoder.jp

解法

ちょうど $n$ 回の操作で終了する確率を $f _ n$ とする。$f _ n$ を直接計算するのは難しいので、終了条件を無視して $n$ 回の操作を行った後に終了条件を満たしている確率 $g _ n$ を考える。また、終了条件を満たしている状態から終了条件を無視して $n$ 回の操作を行った後に再び終了条件を満たしている確率 $h _ n$ を考える。

このとき、$f,g,h$ の母関数 $F,G,H$ について $FH = G$ が成り立つ。

$H = 1$ (この $H$ は $h$ の母関数ではなく盤面の行数である) の場合を考える。$c _ 1$ 個の 1 がある状態から (終了条件を無視して) $n$ 回の操作を行うことで $c _ 2$ 個の 1 がある状態になる確率は、以下で定める行列 $T = (t _ {i, j}) _ {0\leq i\leq W,\ 0\leq j\leq W}$ に対して ${T ^ n} _ {c _ 2, c _ 1}$ と表せる。

$$t _ {i, j} = \begin{cases} j / W & \text{if } i = j - 1 \\ 1 - j / W & \text{if } i = j + 1 \\ 0 & \text{otherwise} \end{cases}$$

$H$ が一般の場合に話を戻すと、上記考察より $i$ 行目に操作した回数 $t _ i$ を固定して数えることで、$g _ n$ は次で表せる。なお、初期盤面で $i$ 行目に存在する 1 の数を $B _ i$ とした。

$$ g _ n = \dfrac{n!}{H ^ n} \sum _ {t _ 1 + \cdots + t _ H = n} \prod _ {i = 1} ^ H \dfrac{1}{t _ i !} {T ^ {t _ i}} _ {A _ i, B _ i}. $$

さて、Cayley–Hamilton の定理より、$T$ の特性多項式を $p(x) \coloneqq \det(xI - T) = \sum _ {k = 0} ^ {W + 1} p _ k x ^ k$ とすれば、任意の整数 $n \gt W$ および $0\leq i, j\leq W + 1$ に対して次が成り立つ。

$${T ^ n} _ {i, j} = - \sum _ {k = 0} ^ W {T ^ {n - (W + 1 - k)}} _ {i, j} \cdot p _ k.$$

これは、$i, j$ を固定して数列 $a ^ {(i, j)}$ を ${a ^ {(i, j)}} _ n\coloneqq {T ^ n} _ {i, j}$ と定めれば、$a ^ {(i, j)}$ が $W + 2$ 項間の線形漸化式を持つことを表している。

また、$a ^ {(i, j)}$ の母関数は、ある $W$ 次以下の $q ^ {(i, j)}(x)$ が存在して $\dfrac{q ^ {(i, j)}}{\mathrm{rev}(p)}$ と表せる。$\mathrm{rev}(p)$ は $p$ の係数列を反転したものであり、$\mathrm{rev}(p)(x) = x ^ {W+1} p(1/x)$ である。$q ^ {(i, j)}$ の計算については、$q ^ {(i, j)}(x) = \left(\mathrm{rev}(p)(x) \cdot \sum _ {n = 0} ^ W {T ^ n} _ {i, j} x ^ n\right) \bmod x ^ {W + 1}$ に従って行えばよい。

以上より、$g _ n$ は次のように表せる。

$$g _ n = \dfrac{n!}{H ^ n} \sum _ {t _ 1 + \cdots + t _ H = n} \prod _ {i = 1} ^ H \dfrac{1}{t _ i !} \lbrack x ^ {t _ i}\rbrack \dfrac{q ^ {(A _ i, B _ i)}}{\mathrm{rev}(p)}. \tag{1}$$

$h _ n$ についても全く同様に考えることで次のように表せることが分かる。

$$h _ n = \dfrac{n!}{H ^ n} \sum _ {t _ 1 + \cdots + t _ H = n} \prod _ {i = 1} ^ H \dfrac{1}{t _ i !} \lbrack x ^ {t _ i}\rbrack \dfrac{q ^ {(A _ i, A _ i)}}{\mathrm{rev}(p)}. \tag{2}$$

適当な整数 $D$ を決めて、式 $(1),(2)$ に基づいて $G, H$ を $\mathrm{mod}\ x ^ D$ で計算すると、$F \bmod x ^ D$ を計算することができる。

さて、$f _ n$ が $K + 1\ (K\leq \lfloor D/2\rfloor)$ 項間線形漸化式を持つと予想する。この予想が正しければ、Berlekamp–Massey のアルゴリズムによって $F(x) = \dfrac{P(x)}{Q(x)}$ なる $P, Q$ を計算できるので、求める期待値を $F'(1) = \dfrac{P'(1) Q(1) - P(1) Q'(1)}{Q(1) ^ 2}$ として計算できる (本当?収束などの条件を確認していないのでかなり怪しい)。

いくらか実験すると $K = HW$ が成り立ちそうだったので、$D = 2HW + \varepsilon$ 程度で定めると AC が取れる。

提出

Submission #46252219 - AtCoder Regular Contest 136

感想

やってることが何もかも怪し過ぎて駄目だろという感じ。

これは余談なんだけど、持っている Berlekamp–Massey が Library Checker で通っていたのにバグっていて、大惨事になってしまった。Hack Case がどういうものかは分かったので取り敢えず issue を投げておいた。

AtCoder Regular Contest 118 F - Growth Rate

問題

atcoder.jp

解法

数列 $X$ を数える代わりに整数列 $D = (D _ 1, \ldots, D _ {N + 1})$ を $D _ 1 \coloneqq X _ 1,\ D _ i \coloneqq X _ i - A _ {i - 1} X _ {i - 1}\ (i \geq 2)$ と定めて $D$ を数える。

$D$ が次の条件をすべて満たすことが、$X$ が条件を満たすための必要十分条件である。

  1. $\displaystyle \sum _ {i = 1} ^ {N + 1} D _ i \prod _ {j = i} ^ N A _ j \leq M $
  2. $D _ 1\geq 1$
  3. $D _ i \geq 0\ (i=2,3,\ldots,N+1)$

$D _ 1$ だけ範囲が異なるのが嫌なので改めて $D _ 1\gets D _ 1 - 1$ および $\displaystyle M\gets M - \prod _ {i = 1} ^ N A _ i$ とすれば、$D$ が次の条件をすべて満たすことが、$X$ が条件を満たすための必要十分条件である。

  1. $\displaystyle \sum _ {i = 1} ^ {N + 1} D _ i \prod _ {j = i} ^ N A _ j \leq M $
  2. $D _ 1\geq 0\ (i=1,2,\ldots,N+1)$

$\displaystyle P _ i\coloneqq\prod _ {j = i} ^ N A _ j$ と定めれば、これは次のように書き直すことができる。

  1. $\displaystyle \sum _ {i = 1} ^ {N + 1} D _ i P _ i \leq M $
  2. $D _ 1\geq 0\ (i=1,2,\ldots,N+1)$

$D$ を桁 DP の要領で数える。より具体的には $\displaystyle D _ i = \sum _ {j = 0} ^ {i - 1} d _ {i, j} P _ {j + 1}$ および任意の $j=1,2,\ldots,i-1$ に対して $0\leq d _ {i, j} \lt A _ j$ を満たすように $D _ i$ を桁 $(d _ {i, 0} \ldots, d _ {i, i - 1})$ へと分解して考える。そして、$j$ 桁目 $d _ {j + 1,j},d _ {j + 2,j},\ldots, d _ {N + 1, j}$ を同時に決めることを $j$ の降順に行っていく。以降、場合分けが煩雑になることを避けるため、$i\leq j$ なる $i, j$ に対して $d _ {i, j} = 0$ とする。

さて、いま、$j \gt t\ (t \geq 0)$ なる $j$ に対する $d _ {\ast, j}$ を既に決めたとする。残りの $j\leq t$ に対する $d _ {\ast, j}$ が満たすべき条件は次のように書ける。

$$\sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ t d _ {i, j} \dfrac{P _ {j + 1}}{P _ {t + 1}} \leq \left\lfloor \dfrac{M - \sum _ {i = 1} ^ {N + 1} \sum _ {j = t + 1} ^ N d _ {i, j} P _ {j + 1}}{P _ {t + 1}} \right\rfloor$$

$0\leq d _ {i, j} \lt A _ j$ の条件より $\sum _ {i = 1} ^ {N + 1} \sum _ {j = t + 1} ^ N d _ {i, j} P _ {j + 1}\lt (N+1) P _ {t + 1}$ であるから、ある整数 $x\ (0\leq x\lt N + 1)$ が存在して上式は次のように表せる。

$$\sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ t d _ {i, j} \dfrac{P _ {j + 1}}{P _ {t + 1}} \leq \left\lfloor \dfrac{M}{P _ {t + 1}} \right\rfloor - x$$

以降 $\left\lfloor \dfrac{M}{P _ j} \right\rfloor$ という形が良く現れるので、これを $M _ j$ と定める。以上より、次のような動的計画法を考えることができる。

$$\begin{aligned}\mathsf{dp}(t, x) \coloneqq{}& \text{全ての $i\gt j \leq t$ なる $i, j$ に対する $d _ {i, j}$ の決め方であって、} \\ & \text{$\displaystyle \sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ t d _ {i, j} \dfrac{P _ {j + 1}}{P _ {t + 1}} \leq M _ {t + 1} - x$ であるようなものの個数}\end{aligned}$$

このとき、答えは $\mathsf{dp}(N, 0)$ である。

まず、$\mathsf{dp}(0, x)$ の計算を考える。これは $\displaystyle \sum _ {i = 1} ^ {N + 1} d _ {i, 0} \leq M _ 1 - x$ なる $d _ {1,0},\ldots,d _ {N+1,0}\geq 0$ の決め方の個数なので、$\displaystyle \mathsf{dp}(0, x) = \binom{N + 1 + M _ 1 - x}{N + 1}$ である。

続いて $t\gt 0,\ 0\leq x\leq N$ に対する $\mathsf{dp}(t, x)$ の計算を考える。$x$ として $0\leq x\leq N$ を満たすものだけを考えれば十分なのは先述の通りである。

$M _ t = \dfrac{M _ {t + 1} - (M _ {t + 1} \bmod A _ t)}{A _ t}$ より $R _ t \coloneqq M _ {t + 1} \bmod A _ t$ と定めると $M _ {t + 1} - x = A _ t M _ t + R _ t - x$ である。

また $\displaystyle \sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ t d _ {i, j} \dfrac{P _ {j + 1}}{P _ {t + 1}} = \sum _ {i = t + 1} ^ {N + 1} d _ {i, t} + A _ t \sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ {t - 1} d _ {i, j} \dfrac{P _ {j + 1}}{P _ t}$ であるから、条件 $\displaystyle \sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ t d _ {i, j} \dfrac{P _ {j + 1}}{P _ {t + 1}} \leq M _ {t + 1} - x$ は次のように整理できる。

$$\sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ {t - 1} d _ {i, j} \dfrac{P _ {j + 1}}{P _ t} \leq M _ t - \left\lceil\dfrac{\left(\sum _ {i = 1} ^ {N + 1} d _ {i, t}\right) + x - R _ t}{A _ t}\right\rceil.$$

$\displaystyle s\coloneqq \sum _ {i = t + 1} ^ {N + 1} d _ {i, t}$ と定めるとさらに次のように表せる。

$$\sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ {t - 1} d _ {i, j} \dfrac{P _ {j + 1}}{P _ t} \leq M _ t - \left\lceil\dfrac{s + x - R _ t}{A _ t}\right\rceil.$$

従って、$\displaystyle \sum _ {i = t + 1} ^ {N + 1} d _ {i, t} = s,\ 0\leq d _ {i, t}\lt A _ t$ なる整数列 $(d _ {t + 1, t}, d _ {t + 2, t},\ldots,d _ {N + 1, t})$ の個数を $f(s)$ とすると、次が成り立つ。

$$\mathsf{dp}(t, x) = \sum _ {s = 0} ^ {(N-t+1)(A _ t - 1)} f(s) \cdot \mathsf{dp}(t - 1, \lceil (s + x - R _ t) / A _ t \rceil).$$

$0\leq x\leq N,\ 0\leq s\leq (N - t + 1)(A _ t - 1),\ 0\leq R _ t\leq A _ t - 1$ の下で、確かに再び $0\leq \lceil (s + x - R _ t) / A _ t \rceil\leq N$ が成り立っていることに注意する。

さて、$0\leq y\leq N$ なる整数 $y$ に対して $y = \lceil (s + x - R _ t) / A _ t \rceil$ を満たす $s$ の範囲は $(y A _ t + R _ t - x) - A _ t\lt s\leq (y A _ t + R _ t - x)$ である。

従って、$f(l,r)\coloneqq \sum _ {s = l + 1} ^ r f(s)$ と定めると次が成り立つ。

$$\mathsf{dp}(t, x) = \sum _ {y = 0} ^ {N} f((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x) \cdot \mathsf{dp}(t - 1, y).$$

区間 $((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x\rbrack$ の長さはちょうど $A _ t$ なので、$f((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x)$ は、$0$ 以上 $A _ t$ 未満の整数からなる長さ $N - t + 2$ の数列であって、総和が $y A _ t + R _ t - x$ であるようなものの個数に一致する。

この数え上げは有名問題で、包除原理を用いることで次を得る。

$$\begin{aligned}& f((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x) \\ &{}=\sum _ {i = 0} ^ {N - t + 2} (-1) ^ i \binom{N - t + 2}{i} \binom{(y A _ t + R _ t - x) - i A _ t + (N - t + 1)}{N - t + 1}.\end{aligned}$$

式の簡単のため $u\coloneqq N - t + 1$ と定めると、$\mathsf{dp}(t, x)$ は次のように表せる。

$$\begin{aligned} \mathsf{dp}(t, x) &{}= \sum _ {y = 0} ^ {N} \sum _ {i = 0} ^ {u + 1} (-1) ^ i \binom{u + 1}{i} \binom{(y A _ t + R _ t - x) - i A _ t + u}{u} \cdot \mathsf{dp}(t - 1, y)\\ &{}=\sum _ {y = 0} ^ {N} \mathsf{dp}(t - 1, y) \sum _ {i = 0} ^ {u + 1} (-1) ^ i \binom{u + 1}{i} \binom{(R _ t - x + (y - i) A _ t) + u}{u} .\end{aligned}$$

さて、$i\gt y$ のとき $R _ t + (y - i) A _ t \lt 0$ より $\displaystyle \binom{(R _ t - x + (y - i) A _ t) + u}{u} = 0$ である。また、$i\gt u + 1$ のとき $\displaystyle\binom{u + 1}{i} = 0$ である。従って、上式における内側の和で $i$ が動く範囲は $0$ 以上 $y$ 以下の範囲としてよい。

つまり、多項式 $a(z) = \displaystyle \sum _ {i = 0} ^ N (-1) ^ i \binom{u + 1}{i} z ^ i$ および $b(z) = \displaystyle \sum _ {i = 0} ^ N \binom{(R _ t - x + i A _ t) + u}{u} z ^ i$ について、$\lbrack z ^ y\rbrack (ab) = \displaystyle \sum _ {i = 0} ^ {u + 1} (-1) ^ i \binom{u + 1}{i} \binom{(R _ t - x + (y - i) A _ t) + u}{u}$ が成り立つ。

計算量を考える。$a(z), b(z)$ は各 $t,x$ に対してそれぞれ $O(N)$ 時間で計算できる。$b(z)$ の計算についてもう少し説明すると、これは $\displaystyle \binom{(R _ t - x + i A _ t) + u}{u} = \binom{(R _ t - (x - 1) + i A _ t) + u}{u}\cdot \dfrac{(R _ t - x + i A _ t) + u}{(R _ t - (x - 1) + i A _ t) + 1}$ を用いて SWAG で積をスライドさせながら求めることでならし $O(N)$ 時間が達成される。

$a(z)$ と $b(z)$ の積の計算は FFT を用いれば $O(N\log N)$ 時間である。これを全ての $t, x$ について行うので、全体 $O(N ^ 3 \log N)$ 時間となる。

本問題では $N$ は最大 $1000$ 程度にもなるので、このままでは実行時間制限に間に合わせることは難しい。

そこで、$A _ i \neq 1$ なる $i$ が高々 $\lfloor \log _ 2 M\rfloor$ 個しかないことに注目する。

$A _ t = 1$ の場合を考えると、任意の $x=0,1,\ldots,N$ について $\mathsf{dp}(t, x) = \mathsf{dp}(t - 1, x)$ が成り立つことが上記の議論から分かる。従って、$A _ t = 1$ の場合の遷移は $O(N)$ 時間で可能であるから、全体の計算時間は $O(N ^ 2 \log N \log M)$ となる。

これでも実行時間制限に間に合わせることは非常に厳しいが、$a(z)$ の FFT 結果をメモ化するなどの定数倍高速化を行うことで正答を得ることができる。

実装

atcoder.jp

AtCoder Regular Contest 124 F - Chance Meeting

問題

atcoder.jp

解法

雑ですが画像だけで伝わってほしい。

問題文の $H,W$ からそれぞれ $1$ ずつ引いて縦の移動回数がちょうど $H$ 回、横の移動回数がちょうど $W$ 回になるよう調整しています。

画像中の「非交差」という表現は不適切かもしれなくて、正確には「始終点以外で同じマスにいることがないようなような移動」です。

補足

$\displaystyle \sum _ {n = 0} ^ \infty \binom{2n}{n} x ^ n = \dfrac{1}{\sqrt{1 - 4x}}$ は、カタラン数の母関数 $\displaystyle \sum _ {n = 0} ^ \infty \dfrac{1}{n + 1} \binom{2n}{n} x ^ n = \dfrac{1 - \sqrt{1 - 4x}}{2x}$ を $x$ 倍して微分することで得られる。

関係無いが、ここから有名な等式 $\displaystyle \sum _ {k = 0} ^ n \binom{2k}{k} \binom{2(n-k)}{n-k} = 4 ^ n$ が得られる (誰かがこの等式の組合せ的な解釈についてしゃべってた気がするけど忘れた)。

なお、A002420 - OEIS によると $\displaystyle\lbrack x ^ n\rbrack \sqrt{1 - 4x} = \dfrac{1}{1 - 2n}\binom{2n}{n}$ らしい。

誤差無し Barrett Reduction

$\mathrm{mod}\ M (\lt 2 ^ {31})$ 上で様々な計算を行うことを考える。このとき、よくある状況として、たくさんの $0\leq a,b\lt M $ なる整数組 $(a,b)$ に対して $a \star b \bmod M\ (\star\in\lbrace +,-,\times\rbrace)$ を計算することになる。

$\star\in\lbrace +,-\rbrace$ の場合は次のようにして容易に剰余演算を回避できる。

// 1 <= M < 2^31
uint32_t M = ...;
// 0 <= a,b < M の仮定下で a+b mod M を計算
uint32_t add(uint32_t x, uint32_t y) {
    uint32_t c = a + b;
    return c >= M ? c - M : c;
}
// 0 <= a,b < M の仮定下で a-b mod M を計算
uint32_t sub(uint32_t x, uint32_t y) {
    uint32_t c = a + (M - b);
    return c >= M ? c - M : c;
}

$a\times b \bmod M $ を剰余演算無しで計算するのは容易ではないが、Barrett Reduction はこれを実現する方法の $1$ つとして知られている。1 2

状況を少し一般化して、ある正整数 $M $ が固定されており、たくさんの非負整数 $x$ に対して $x \bmod M $ を計算することを考える。

また $x$ の上界 $X\ (\geq M)$ が定まっているものとする。つまり、以下では $0 \leq x\leq X$ を仮定する。$a\times b\bmod M $ の例では $X$ として $(M - 1) ^ 2$ 以上の値を与えることになる。

Barrett Reduction では、初めに $M,X$ のみに依存する (剰余演算を含む) 簡単な前計算を行うことで、各クエリを剰余演算を用いずに処理する。

本記事では、Barrett Reduction における商の近似値の誤差を $0$ にすることで条件分岐を解消する手法について説明する。なお、この手法が高速化に寄与するかは不明である (誰か実験してほしい)。

記号

定数と変数の区別を明確にするため、以下では $x$ に依存しない値を大文字で定義し、$x$ に依存する値を小文字で定義する。

通常の (?) Barrett Reduction

目標は $r\coloneqq x\bmod M $ を計算することである。$r = x - M\cdot \lfloor x/M \rfloor$ より、$q\coloneqq \lfloor x/M \rfloor$ を求められれば十分である。

Barrett Reduction では次のような方針で $q$ を計算する。

  1. $q$ の近似値 $q'$ を計算する。
  2. $q'$ を微調整して $q$ を得る。

Step 1. 近似値 $q'$ の計算

適切に正整数 $K$ を固定して $L\coloneqq \lceil 2 ^ K / M\rceil$ とする 3。$L$ は二進法における $1/M $ の小数第 $K + 1$ 位以下を切り上げて $2 ^ K$ 倍した値である。$K$ は (大文字で定義していることからも分かるように) $x$ に依存しない値であり、その決め方は後述する。$L$ も $x$ に依存しない値であることに注意する。

$q'\coloneqq \lfloor xL / 2 ^ K \rfloor$ を $q$ の近似値として、その誤差を評価する。任意の実数 $a$ に対して $a\leq \lceil a\rceil \lt a + 1$ および $a - 1\lt \lfloor a\rfloor \leq a$ が成り立つので、次が従う。

$$\begin{alignedat}{3} 2 ^ K / M & {}\leq{} & L & {}\lt{} & 2 ^ K / M + 1, \\ \lfloor xL/2 ^ K \rfloor & {}={} & q' & {}\leq{} & xL/2 ^ K. \end{alignedat}$$

これらを併せて次の式 $(1)$ を得る ($L,x,M $ は全て $0$ 以上であることに注意)。

$$q = \lfloor x/M \rfloor \leq q' \lt x/M + x/2 ^ K.\tag{1}$$

従って、$x\leq 2 ^ K$ ならば $q'\in \lbrace q,q+1\rbrace$ が成り立つ。つまり、$K \coloneqq \lceil \log _ 2 X \rceil$ と定めれば必ず $q'\in \lbrace q,q+1\rbrace$ が成り立つ。

実装に当たっては以下の点に注意せよ:

  • $L$ を計算するためには除算を行う必要があるが、これは $x$ に依存しないので前計算して使いまわすこと。クエリ毎に毎回 $L$ を除算で計算すると Barrett Reduction を用いる意味が無くなるので注意。
  • $q'$ の計算に当たって、$2 ^ K$ での切り捨て除算はシフト演算を用いること。具体的には q_prime = (x * L) >> K と計算すること。これも除算を使うと Barrett Reduction を用いる意味が無くなってしまう。(最適化でシフト演算に置き換えてくれるかもしれないけど)

Step 2. $q$ の計算

$q'\in \lbrace q, q+1\rbrace$ であるから、$r' \coloneqq x - q' M $ とすれば次のように $q$ を計算できる。同時に $r$ も計算できる。

  • $r'\lt 0$ の場合: $q = q' - 1,\ r = r' + M $ である。
  • $r' \geq 0$ の場合: $q = q',\ r = r'$ である。

$r'$ の値によって条件分岐が生じることに注意する。

計算途中に現れる値の最大値

オーバーフロー回避のために、計算途中に現れうる値の最大値 $U$ を見積もっておく。これは $U = XL = X\lceil 2 ^ K / M\rceil$ である。

例えば $X\lt 2 ^ {63}$ の場合は $K = 64$ とすれば $K \geq \lceil \log _ 2 X\rceil$ を満たす。このとき $U\lt 2 ^ {127}$ であるから、計算途中の値を $128$ bit 整数型で持てばオーバーフローは起こらない。

$r'$ は負の値を取り得るので符号無し整数型を用いて実装する場合は十分注意すること。なお、$x$ と $q'M $ の大小比較で条件分岐すればこの問題を避けることができる。

(本当は負方向のオーバーフロー回避のために最小値も確認すべきだが、省略する)

誤差無し Barrett Reduction

実は、上記アルゴリズムにおいて $K$ をさらに大きく取ることで $q'$ は $q$ と必ず一致するようにできる。これによって Step 2 を省略でき、即ちクエリ処理時に条件分岐が一切現れないようにできる。

$q'$ と $q$ が必ず一致するような $K$ の十分条件を考える。式 $(1)$ より、全ての $0\leq x\leq X$ について $x / M + x/2 ^ K \leq \lfloor x/ M \rfloor + 1$ を満たすことが十分である。

$x / M $ の小数部分 $t\coloneqq x/M - \lfloor x/M\rfloor$ について $0\leq t\leq 1 - 1/M $ が成り立つ。従って、$X / 2 ^ K \leq 1/M $ 即ち $K \geq \lceil \log _ 2 (MX)\rceil$ は十分条件である。

modint

冒頭の話に戻り、Barrett Reduction を modint 実装に用いることを考える。

固定された $M\ (\lt 2 ^ {31})$ とたくさんの $0\leq a,b\lt M $ なる $a,b$ に対して $a\times b\bmod M $ を計算したいのであった。$0\leq a\times b \leq (M - 1) ^ 2$ であるから $X = (M - 1) ^ 2\ (\lt 2 ^ {62})$ とすれば全ての入力をカバーすることができる。

また、$K = 96$ と定めると、$M \lt 2 ^ {31},\ X\lt 2 ^ {62}$ の下で誤差無し Barrett Reduction の条件 $K \geq \lceil \log _ 2 (MX)\rceil$ を満たす。

計算途中に現れる値の最大値は $U = (M - 1) ^ 2\lceil 2 ^ {96} / M\rceil$ であり、次が成り立つ。

$$\begin{aligned} U & \lt{} (M - 1) ^ 2 (2 ^ {96} / M + 1) &\quad &(\because \lceil 2 ^ {96} / M\rceil \lt 2 ^ {96} / M + 1)\\ & ={} 2 ^ {96} M - 2\times 2 ^ {96} + 2 ^ {96} / M + (M - 1) ^ 2 \\ & \lt{} 2 ^ {127}. &\quad &(\because 1\leq M\lt 2 ^ {31}) \end{aligned}$$

つまり、計算途中の値は $128$ bit 整数型で持つことができる。

modint 実装例

以上を踏まえて modint を実装すると次のようになる (※動作確認はしてません。雰囲気を伝えるだけのものと解釈してください。)

class modint {
    struct barrett {
        static constexpr int K = 96;
        uint32_t M;
        __uint128_t L;
        barrett(uint32_t M) : M(M), L(((__uint128_t(1) << K) - 1) / M + 1) {}
        // a*b mod M
        uint32_t mul(uint32_t a, uint32_t b) const {
            uint64_t c = uint64_t(a) * b;
            uint64_t q = (c * L) >> K;
            return c - q * M;
        }
    };
    static inline barrett _bt{1};
    uint32_t _val;
public:
    // ... (コンストラクタやメソッドなど)

    static void set_mod(int M) {
        assert(1 <= M);
        _bt = barrett(M);
    }

    modint& operator*=(const modint & b) {
        _val = _bt.mul(_val, b._val);
        return *this;
    }
};


  1. $M $ がコンパイル時定数である場合は a * b % M などと書いても最適化によって剰余演算を回避するバイナリが吐かれる (と思われる)。
  2. 他の方法としては、例えば Montgomery 乗算など
  3. バレットの還元アルゴリズム - Wikipedia のように $L\coloneqq \lfloor 2 ^ K / M\rfloor$ とした場合に誤差無し Barrett Reduction を作ろうとするとどうなるかを考えると面白いかもしれません。

rep マクロ

#define rep(i, n) for (int i = 0; i < n; ++i)

例えば rep(i, k & 1) とすると for (int i = 0; (i < k) & 1; ++i) と等価であるため、意図しない動作を引き起こします: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

これを避けるには #define rep(i, n) for (int i = 0; i < (n); ++i) とすればよいです: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

#define rep(i, l, r) for (int i = (l); i < (r); ++i)

$l$ や $r$ は 32 bit 整数で表せないほど大きな整数かもしれません。

#define rep(i, l, r) for (long long i = (l); i < (r); ++i)

$l$ や $r$ は 64 bit 整数で表せないほど大きな整数かもしれません。

#define rep(i, l, r) for (__int128_t i = (l); i < (r); ++i)

速度的にあまり嬉しくなさそう。

#define rep(i, l, r) for (decltype(r) i = (l); i < (r); ++i)

r は参照型かもしれません: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

#define rep(i, l, r) for (std::remove_reference_t<decltype(r)> i = (l); i < (r); ++i)

rconst 修飾されているかもしれません: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

#define rep(i, l, r) for (std::decay_t<decltype(r)> i = (l); i < (r); ++i)

std::decay_t<decltype(r)>l を表現できない場合にバグります: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

こうすると signedness さえ一致していれば正しく動く: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

signedness が一致していない場合は unsigned の方が選ばれるので、rep(i, -1, 3U) などとするとバグる: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

signedness が一致しない場合にコンパイルエラーを吐かせることもできる: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

追記

上記実装で壊れるケースを頂きました:

r の計算が比較のたびに走って計算量が悪化することがあります。

#define rep(i, l, r) for (rep_int_t<std::decay_t<decltype(l)>, std::decay_t<decltype(r)>> i = (l), end_i = (r); i < end_i; ++i)

r を先に評価することで回避できる: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

ただし、end_i という変数が既に外側で使われていた場合に変数の shadowing が起こる: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

対処として考えられそうなのは次の 2 種類。

  1. 衝突しなさそうな変数名を使う
  2. rep 用の構造体を作る

例えば 2 であれば次のように実装できる: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

この場合、ループ中に i を変更したときの挙動が今までとは異なることに注意する: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

この挙動の違いは非自明なので、iconst で宣言してあげるのがよさそう: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

構造体にする利点

次のように同じ rep マクロで引数の個数によって挙動を使い分けたいとする。

  • rep(i, n) : $i = 0, 1, \ldots, n - 1$ の順に走査
  • rep(i, l, r) : $i = l, l +1, \ldots, r - 1$ の順に走査

c - Overloading Macro on Number of Arguments - Stack Overflow のような方法もあるが、構造体であればコンストラクタを増やして次のように簡単に書ける: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

一番初めに挙げた優先順位の問題は起こらないはず (多分...)

コードがかなり長くなってしまうのが嫌なら l, r の型の一致を強制することで短くできて、次のようになる: [C++] gcc HEAD 14.0.0 20230712 (experimental) - Wandbox

まとめ

この記事の結論としては次の 2 つをおすすめする。