noshi91のメモ

データ構造のある風景

ICPC Yokohama 2023 参加記

リハーサル

MLE を見ようとして大量にメモリを取ったら PC が落ちたりして面白かった。

本番

ある程度関わった問題を何となく時系列っぽく並べる

問題文と解説: Regional | ICPC 2023 Asia Yokohama Regional

B

親の顔より見た平均の式変形。 c の小ささを利用しようかと考えたが数秒で思いつかなかったので却下した。

D

ほとんどやるだけ。連続部分列の高速な一致判定が欲しくなったので  \Theta(n  ^ 2) で LCP のテーブルを作っておいた。

H

potato と一緒に考察する。 しばらく悩んでいたが、potato が「最小カットじゃないか?」と発言したのでそこから数秒で解けた。 最小カットで非自明な問題を久々に見た気がするが、すぐに解けるべきだった。

K

線分の端点が正方形の内側にあるようなクエリのことが脳から抜け落ちていて、 x 軸と  y 軸のそれぞれに対して座標を特定しようとするとクエリがギリギリ溢れるななどと考えていた。 その結果アルゴリズムも複雑になり、実装に手間取ったうえにバグらせてしまった。

J

最初の方で potato から「min plus の畳み込みは高速になるか?」と聞かれ、いかにも担当っぽい問題だったので貰う。 木上の最小費用流になることはすぐに分かり既知のアルゴリズムを適用しようとしたのだが、頂点の湧き出しに対するコストが  f m ^ 2 なので辺の本数が  2 乗になってしまう。 こういうときはアルゴリズムの動作を観察して適切に高速化すると (例えば等差数列加算とか) 上手く行くことが多いのだが、今回は等差数列の insert だったので解けそうで解けない。 得意分野のはずという気持ちでかなり時間を消費してしまったのだが、結局最後まで残ってしまう。 最小費用流なので多分何かしらの貪欲だろうということで potato Hemimor に任せているとそれっぽい貪欲が提案されたので実装してもらい、通った。

結局、等差数列の inset を愚直にやっても上手く計算量が抑えられていたらしい。 この分野に関連するアルゴリズムはかなり触っており自信があったので、ショックだった。

C

とりあえずポリアはする。

回転を同一視しない数え上げを考えると、カタラン数を少し弄ったようなものが出てくる。 ラグランジュ反転の思い出しに自身が無かったのと、最後の方まで畳み込みなしで押せる自信が無かったので FPS で書いておく。

周期が決まっている場合の数え上げだが、当初スタックに色を積んでいく方針で考えていたので  1 周期で全部消えるケースしかないんじゃないかと予想した。 この方針で実装するとサンプル  1 が合わないので見てみると、消えないケースがある。 サンプルすら合わない解法を実装し始めてしまうのは悪い癖なのだが直りそうもない。 しばらく考えると、隣り合う同じ色を消す方針を思いつき、奇数長回文が残る場合が必要十分と分かった。 これも FPS で立式できるので書いて出してみると TLE。

ここからの高速化の作業は悪夢だった。 まず FPS の立式を偶奇で分けることで長さを半分にする (TLE)。 FPS の inv を定数倍が良好な実装に変える (TLE)。 同じ列の DFT を繰り返し計算している部分を使いまわすようにする (TLE)。 その他あり得そうな部分をひたすら高速化して回る (TLE)。

結局この方針は十分に高速化された FFT でギリギリ通せるラインのものだったようだ。 SpeedStar にはどのみち負けていたのだが、食らいつけなくて残念。

総評

不調気味だった気がするが  3 位には落ちなかったし、 1 位は好調でも難しかっただろうと思うので、全完できなかったことだけが心残りになった。 海外に行くのが嫌いなので playoff 確定でげんなりしているが頑張りたい。

懇親会

例の大きい名札を付けてうろうろし、色々な人に話しかけてもらった。 初めて会う人もたくさんいたので嬉しい。

懇親会でアルゴリズムの話をする難しさをこの手のイベントの度に感じていて、「最近良い感じのアルゴリズムありましたか?」から切り出せる相手があまり多くない。 次の機会があれば、hotman さんに倣って懇親会で話せるアルゴリズムのお品書きを作っておくべきかもしれない。 fenwick tree の最悪計算量を半分にする方法や、競プロで使えそうにない最近のフローアルゴリズム群などについて話した。

終わりに

準備に携わったスタッフの皆様、ありがとうございました。

Aliens DP における二分探索の色々

とにかく使いたい人向け

スコアとペナルティがどちらも整数とする。 *1

ペナルティを  p にしたときの最適値を  g(p) と置く。

事実  1

ペナルティが  p の時の最適解がペナルティを受ける回数としてあり得るものは、区間  \lbrack l _ p , r _ p \rbrack になっている。

事実  2

最小化 (ペナルティの度にスコアに  + p ) なら、 r _ {p + 1} = l _ p = g(p + 1) - g(p) である。

最大化 (ペナルティの度にスコアに  - p ) なら、 r _ {p + 1} = l _ p = g(p) - g(p + 1) である。

使い方

 g(p) の計算の際にペナルティ回数の復元は不要で、 g(p + 1) - g(p) (最大化なら逆) を返せばいい。

詳しめの解説

概要

Aliens DP などで二部探索をするときに用いる複数の手法を紹介する。

手法 利点 欠点
ペナルティ回数の管理 分かりやすい。ペナルティ回数を求める際に計算時間がほとんど増えない問題では最速。 ペナルティ回数の管理が面倒。tie-breaking に気を使わないと壊れる。
傾き二分探索 回数の管理が不要な中では最も分かりやすく、三分探索より速い。 最速ではない。
三分探索 凸じゃないケースで誤魔化せる可能性が傾き二分探索より高いかもしれない。 傾き二分探索のほぼ下位互換。
黄金分割探索 回数の管理で計算時間が嵩む場合最速。 ほぼライブラリ化前提、その場で書くのは苦しい。

問題設定

下に凸な関数  f \colon \mathbb{Z} \to \mathbb{Z} \cup \lbrace + \infty \rbrace, ~ f \neq + \infty があり、 g \colon \mathbb{Z} \to \mathbb{Z} \cup \lbrace - \infty \rbrace, ~ p \mapsto \inf \limits _ x \left ( f(x) - px \right ) と定義する *2 g(p) は、 y = f (x) のグラフに傾き  p の接線を引いた際の  y 切片の値と解釈することもできる。 我々の目標は、任意の  p \in \mathbb{Z} に対して  g \left ( p \right ) が高速に計算できるとき、与えられた  \tilde x \in \mathbb{Z} に対する  f \left ( \tilde x \right ) を計算することである。

ペナルティ回数の管理

 \inf \limits _ x \left ( f(x) - px \right ) を達成する最小の  x x ^ {\ast} _ p とする。  g \left ( p \right ) の評価の際に  x ^ {\ast} _ p を同時に計算できる場合を考える。

 x ^ {\ast} _ p p に対して単調増加であるから、二分探索を用いて  x ^ {\ast} _ {\tilde p} \leq \tilde x \lt x ^ {\ast} _ {\tilde p + 1} を満たす  \tilde p を発見できる。 すると  \tilde x \in \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \tilde p x \right ) であるから、 g \left ( \tilde p \right ) = f \left ( \tilde x \right ) - \tilde p \tilde x から  f \left ( \tilde x \right ) の値が分かる。

 x ^ {\ast} _ p の定義を  \inf を達成する最大の  x とした場合も、 x ^ {\ast} _ {\tilde p - 1} \lt \tilde x \leq x ^ {\ast} _ {\tilde p} を満たす  \tilde p を発見すれば同様である。 一方で、 \min を達成する  x を何か  1 つ計算できるとした場合、この方法では  f \left ( \tilde x \right ) を特定できない *3

傾き二分探索

 \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - p x \right )  \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 1 \right ) x \right ) の共通部分は  1 元集合であり、 \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 0.5 \right ) x \right ) と一致する。 証明は容易なので省略する。

上記の事実から、 \left \lbrace x ^ {\ast} _ {p + 0.5} \right \rbrace = \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 0.5 \right ) x \right ) と定義すると以下の式が成り立つ。

\begin{align} g ( p ) &= f \left ( x ^ {\ast} _ {p + 0.5} \right ) - p x ^ {\ast} _ {p + 0.5} \cr g ( p + 1 ) &= f \left ( x ^ {\ast} _ {p + 0.5} \right ) - \left ( p + 1 \right ) x ^ {\ast} _ {p + 0.5} \end{align}

ここから  x ^ {\ast} _ {p + 0.5} = g ( p ) - g ( p + 1 ) が分かるため、 g の評価ができるだけで  x ^ {\ast} _ {p + 0.5} を得ることができる。 後は二分探索を用いて  x ^ {\ast} _ {\tilde p - 0.5} \leq \tilde x \leq x ^ {\ast} _ {\tilde p + 0.5} を満たす  \tilde p を発見すれば、ペナルティ回数の管理の節と同様である。

本節の内容は  \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - p x \right ) = \left \lbrack x ^ {\ast} _ {p - 0.5} , x ^ {\ast} _ {p + 0.5} \right \rbrack を意識すると理解しやすいかも知れない。

また、ペナルティ回数の管理を用いて  x ^ {\ast} _ {p + 0.5} を直接求めることで本節の内容を実行しても良い。 その場合、 \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 0.5 \right ) x \right ) = \mathop{ \mathrm{arg\,min} } \limits _ x \left ( 2 f(x) - \left ( 2 p + 1 \right ) x \right ) を用いると整数の範囲で計算できる。

三分探索

 g の定義より、任意の  p について  g(p) \leq f \left ( \tilde x \right ) - p \tilde x である。 これを移行して  f \left ( \tilde x \right ) \geq g(p) + \tilde x p としたのち両辺を  p についての  \sup をとると、 f \left ( \tilde x \right ) \geq \sup \limits _ p \left ( g(p) + \tilde x p \right ) を得る。  f \left ( \tilde x \right ) \lt +\infty として、 f の凸性から  g(p) = f \left ( \tilde x \right ) - p \tilde x を満たす  p が存在するため前述の式は等号で成立し、

\begin{align} f \left ( \tilde x \right ) = \sup \limits _ p \left ( g(p) + \tilde x p \right ) \end{align}

である。  g は線形関数の  \inf で定義されるから上に凸であり、従って  p \mapsto g (p) + \tilde x p も上に凸であるから、三分探索を用いて最大化すればよい。

傾き二分探索の節の内容は、 p \mapsto g (p) + \tilde x p の最大化を傾きの二分探索で行っているのと本質的に同じである。 よって、傾き二分探索の節の内容を押し進めて三分探索に行き着くことも可能である。

黄金分割探索

三分探索の節の内容をそのまま黄金分割探索にする。

問題設定の変種

 f \colon \mathbb{R} \to \mathbb{R} \cup \lbrace + \infty \rbrace が下に凸かつ狭義単調増加で、 g (p) = \inf \limits _ x \left ( f(x) - px \right ) が高速に計算できるとする。 与えられた  k に対し、 f \left ( \tilde x \right ) = k を満たす  \tilde x を求めよ。

 f が離散ではなく連続になっていることに注意。 これは  f \left ( \tilde x \right ) = k となる  \tilde x が存在した方が議論が簡潔になるためであり、離散の設定なら線形に連続に拡張してしまえば良い。

ペナルティ回数の管理, 傾き二分探索

 x ^ {\ast} _ {p} が分かれば  g (p) = f \left ( x ^ {\ast} _ p \right ) - p x ^ {\ast} _ p から  f \left ( x ^ {\ast} _ p \right ) が分かるため、二分探索を行うことができる。

三分探索, 黄金分割探索

まず、 f を平行移動することで  k = 0 とする。 すると、 f(x) のグラフと  x 軸との交点の  x 座標を求める問題となる。

 f x 軸との交点における傾き  \tilde p を求めることを考える。  h \colon \mathbb{R} _ {\gt 0} \to \mathbb{R} \cup \lbrace + \infty \rbrace h(p) = - \dfrac{g (p) }{p} で定める。  h(p) は、 f (x) に引いた傾き  p の接線の  x 切片である。  f が単調かつ凸であるから、 h は単峰で  \tilde p において最小値を取り、三分探索を適用することで求まる。 このとき  \tilde x = h \left ( \tilde p \right ) である。

正確には、 f x 軸との交点において微分可能とは限らないため、 \tilde p は任意の劣勾配ということになる。

注意

  • 二分探索や三分探索の初期値は、問題毎に考えて行うしかない。
  •  g(p) = - \infty の場合  x ^ {\ast} _ p などが定まらないが、これも問題毎に処理する。ただし  f の定義域が有界区間であることが多く、その場合  g は常に有限になる。また、 g(p) = - \infty となるのは  p が大きすぎるか小さすぎるかのときなので、 x ^ {\ast} _ p を適切に  \pm \infty にすることでも回避できる。

余談

本稿で定義した  g の符号を反転した関数は  f の凸共役と呼ばれ、 f ^ {\bullet} で表される。 凸共役に関する各種の定理を認めると、本稿の内容は大幅に簡略化される。

関連する記事

Aliens DP で辺の本数が区間で指定される場合 - noshi91のメモ

*1: スコアが実数ならペナルティも実数にしないと駄目で、tie-breaking は自由にやってよい

*2: p がペナルティというよりはボーナスになっているが、argmin が単調減少ではなく単調増加の方が見た目が分かりやすいためである。

*3: 正確には、次節の内容と組み合わせると特定できる

省メモリな区間 add 区間 min

概要

長さ  n の数列  a に対する区間 add と区間 min を時間計算量  \Theta ( \log ( n ) ) で処理する、空間計算量  cn のデータ構造を説明する。 ただし、一連の操作における  \lvert a _ i \rvert の最大値を  X としたとき、 \lbrack -2X , 2X \rbrack の範囲の数を管理するために必要な空間計算量を  c としている。

例: 長さ  10 ^ 5 の数列に  \pm 10 ^ 9 の加算クエリが  10 ^ 6 回ほど行われる状況では、 64 bit 整数  10 ^ 5 個のメモリで処理できる。

構造

簡略化のため、 n 2 の冪の場合をまず考える。 Segment Tree と同様に、配列を用いて二分木の構造を作る。 まず、各ノードに対して対応する区間内の最小値を考える (この値は実際には保持しない)。 続いて、各ノードに対して左子  - 右子を計算し、これを保持する。 葉には子が無いので、データは保持しない。 すると下図のような構造となる。

この右側のデータに加えて全体の最小値も持つことで、Segment Tree に対するものと似たようにクエリの処理が可能となる。 例えば区間 add では、データを適切に更新しつつ、着目しているノードの区間の最小値の増分を返すような関数を設計すると良い。 区間  \lbrack l , r )  x を加算するクエリは以下のようになる

  1. 着目しているノードの区間 \lbrack l , r ) に完全に含まれているなら、 x を返す。
  2. 着目しているノードの区間 \lbrack l , r ) と交差しないなら、 0 を返す。
  3. 左子と右子に対して再帰的に加算クエリを行う。それぞれの返り値から着目しているノードの値を更新し、返り値を計算する。

 n 2 の冪でない場合も、空間計算量が  2n の Segment Tree と同様に完全二分木でない二分木構造を作ればよい。 ノードの個数は  n - 1 になるため、全体の最小値の管理と併せて  n 個の整数を管理することになる。

上限付き単調増加列の数え上げ

概要

  • 長さ  N の広義単調増加な整数列  0 \leq A _ 0 \leq A _ 1 \leq \dots \leq A _ {N - 1} \leq M が与えられたとき、 0 \leq B _ i \leq A _ i を満たす広義単調増加な整数列  B の個数  \bmod P \Theta ( (N + M) ( \log (N + M) ) ^ 2 ) の時間計算量で計算できる。
  • 上記の設定に加えて  B _ i の下限も与えられる場合も、同じ計算量で計算できる。
  • (, ), ? からなる長さ  N の文字列が与えられたとき、?( または ) に置き換えて得られる文字列であって valid な括弧列になるものの個数  \bmod P \Theta ( N ( \log (N) ) ^ 2 ) の時間計算量で計算できる。
  • 見た目の異なる  2 つのアルゴリズムがある。どちらが高速かは未検証。

単調増加列の数え上げ

列の数え上げを経路の数え上げで言い換える。 例えば  N = 5, M = 6, A = (0, 2, 2, 3, 5) の場合、以下の図において赤い点から右か上に進んで青い点に至る経路の個数と等しい。

経路の数え上げなので、 \Theta(NM) 掛けた単純な DP でひとまず答えは出ることが分かる。 さらに分割統治法の都合上、左下から右上への経路の数え上げというのを強めて以下の問題を考える。

上図のような階段状の盤面が与えられ、赤い点にはそれぞれ何らかの値が書き込まれている。 ここから経路数え上げの DP を行ったとき、最終的にそれぞれの青い点に書き込まれる値を求めよ。

これを分割統治法で解く。 まず、下辺の中心から上に直進し、盤面に突き当たったら右に直進すると、下図のように長方形領域が切り取られる。(盤面を簡略化して三角形風に描いている)

そして、左下の領域において DP を計算する (赤矢印)。 これは再帰的な呼び出しで実現できる。 次に、中央の長方形部分において DP を計算する (青矢印)。 これは寄与が二項係数で表されるので、丁寧に式変形をすると畳み込みになることが分かる。 最後に、右上の領域において再帰的に DP を計算する (緑矢印)。

時間計算量は  \Theta( (N + M) ( \log(N + M) ) ^ 2) である。

下限もついている場合

上図のように盤面を分割すると、分割されたそれぞれの領域を前述したアルゴリズムで順々に計算していくことができる。 時間計算量は再び  \Theta( (N + M) ( \log(N + M) ) ^ 2) となる。

括弧列の数え上げ

(, ), ? からなる長さ  N の文字列が与えられる。?( または ) に置き換えて得られる文字列であって valid な括弧列になるものの個数を数えよ。

文字列に含まれる ? の個数を  n とする。 長さ  n の数列  a を考え、 a _ i を、 i 個目の ?( に置き換えるなら  1, ) に置き換えるなら  -1 と定義する。 すると、 1, -1 から成る数列  a が valid な括弧列に対応する条件は「 \sum _ {i = 0} ^ {k} a _ i \geq t _ k という条件が各  0 \lt k \lt n についてあり、 k = n については等式で成り立っている」といった形のものになる。 これを図で表してみると、以下のような図における経路の数え上げとなる。

×印で禁止されている場所以外にも、赤点から青点への経路として通ることが不可能な場所がある。 そのような無駄な部分を削除すると、下図のようになる。

これは回転すると単調増加列の数え上げと同じ形である。

亜種のアルゴリズム

本節は、https://codeforces.com/contest/1770/problem/G この問題の公式解説を私が再解釈して記述したものである。

単調列の数え上げも、括弧列の数え上げも、上図のような経路の数え上げに帰着できるのであった。 括弧列の節で言及したとおり、上の方に到達不能な部分があるが、無限に広がっていると思う方が都合がよいのでそのままにしている。 この図では右下か右上かに進むことを繰り返す経路が得られるが、これが右と右上になるように盤面を線形変換すると下図のようになる。

ここで、分割統治のために一般化した問題を考える。

上図の赤辺部分に値が書き込まれている。右か右上に進めるとして経路数え上げの DP をしたとき、青辺部分に書き込まれる値を求めよ。

赤辺を青辺の下端の高さで上下に分け、それぞれについて青辺への寄与を計算する。 赤辺上部の青辺への寄与は単に二項係数で表されるので、畳み込みで計算できる (緑矢印)。 赤辺下部の青辺への寄与は分割統治法により計算する (水色矢印)。

このアルゴリズムの優秀な部分は、 n 2 冪になるように調整すれば、畳み込みの部分で使いまわしが効くようになることである。 記事中前半で紹介したアルゴリズムでは、 N, M 2 冪にしても畳み込みの部分では中途半端な値が出てきてしまうので使いまわしができない場所がある。 しかし上手くやればできるかもしれないし、全体を通してどちらが高速なのかは未検証。 また、このアルゴリズムを上限と下限のある単調列の数え上げに適用できるかも不明である。

メモ: q-二項係数

q-二項係数や Gaussian binomial coefficient などと呼ばれる概念がある。

q-Binomial Coefficient -- from Wolfram MathWorld

Gaussian binomial coefficient - Wikipedia

Q-analogs, q-Lucas theorem and q-Catalan numbers

詳細は上記のページなどを参照すると良い

基本事項

非負整数  m, n に対し以下の記号を定義する。

\begin{align} \lbrack n \rbrack _ q &:= 1 + q + q ^ 2 + \cdots + q ^ {n - 1} = \frac{1 - q ^ n}{1 - q} \cr\cr \lbrack n \rbrack _ q ! &:= \lbrack 1 \rbrack _ q \lbrack 2 \rbrack _ q \cdots \lbrack n \rbrack _ q \cr\cr \binom{m}{n} _ q &:= \frac{\lbrack m \rbrack _ q !}{\lbrack n \rbrack _ q ! \lbrack m - n \rbrack _ q !} = \frac{(1 - q ^ m)(1 - q ^ {m - 1}) \cdots (1 - q ^ {m - n + 1})}{(1 - q)(1 - q ^ 2) \cdots (1 - q ^ n)} \end{align}

  •  \binom{m}{n} _ q q多項式となる (分数が割り切れる)
  • 有限体  \mathbb{F} _ q 上の  n 次元ベクトル空間の  k 次元部分空間の個数は  \binom{n}{k} _ q
  •  0 a 個、 1 b 個並べた数列であって転倒数が  k のものの個数は  \lbrack q ^ k \rbrack \binom{a + b}{a} _ q
  • 各要素が  0 以上  m 以下で、長さが  n、合計が  k の広義単調増加列の個数は  \lbrack q ^ k \rbrack \binom{n + m}{n} _ q
  • 各要素が  0 以上  m 未満で、長さが  n、合計が  k の狭義単調増加列の個数は  \lbrack q ^ {k - n (n - 1) / 2} \rbrack \binom{m}{n} _ q

\begin{align} \prod _ {i = 0} ^ {n - 1} (1 + q ^ i t) = \sum _ {k = 0} ^ {n} q ^ {k (k - 1) / 2} \binom{n}{k} _ q t ^ k \end{align}

 n _ 0, k _ 0 \lt d とする。

\begin{align} \binom{n _ 1 d + n _ 0}{k _ 1 d + k _ 0} _ q \equiv \binom{n _1}{k _ 1} \binom{n _ 0}{k _ 0} _ q \pmod {\Phi _ d(q)} \end{align}

 \binom{n _ 1}{k _ 1} の部分は通常の二項係数であることに注意。  \Phi _ d d 番目の円分多項式

 q-二項係数は、通常の二項係数で数えられる対象を更に何らかのパラメータで分類して、それを母関数として表したものと考えることも出来る。 言い換えると、 q = 1 を代入すると通常の二項係数と一致する。 例えば転倒数の例では、全ての転倒数について数列の個数を足し合わせれば当然  \binom{a + b}{a} 個になるから  \lbrack q ^ k \rbrack \binom{a + b}{a} _ q が正解だろうと推測できて、 \lbrack q ^ k \rbrack \binom{a}{b} _ q のようなものが答えになることはあり得ない。

計算について

  •  q が定数なら、通常の二項係数のように階乗とその逆元のテーブルを作れば高速に二項係数が計算できる。ただし、 n q \mathbb{F} _ p における位数より大きい場合、q-Lucas の定理を使う必要がある。
  •  q不定元として、 \binom{n}{r} _ q k 次まで計算することは FPS の exp と log を用いれば  \mathrm{O}(k \log(k) ) で行える。

証明

 0 a 個、 1 b 個並べてできる全ての数列を考える。 これらの転倒数の分布を母関数として表す。すなわち

\begin{align} I (q) := \sum _ {k = 0} ^ {\infty} \left \lvert \left \lbrace 0 \text{ を } a \text{ 個、} 1 \text{ を } b \text{ 個並べてできる、転倒数が } k \text{ の数列} \right \rbrace \right \rvert q ^ k \end{align}

と定義する。

 \lbrace 1, 2, \dots, a + b \rbrace の順列の全てについて同様に、その転倒数の分布を母関数として表す。

\begin{align} P (q) := \sum _ {k = 0} ^ {\infty} \left \lvert \left \lbrace \lbrace 1, 2, \dots, a + b \rbrace \text{ の順列のうち、転倒数が } k \text{ であるもの} \right \rbrace \right \rvert q ^ k \end{align}

 P(q) 2 通りの方法で計算する。

計算方法  1

挿入 DP のように考える。 \lbrace 1, 2, \dots, t - 1 \rbrace までの並びを決めたとき、次に  t をどこに差し込むか考えると、 t 通りの選択肢がありそれぞれで転倒数が  0, 1, \dots, t - 1 増加する。 これは母関数に  1 + q + \dots + q ^ {t - 1} = \lbrack t \rbrack _ q を掛けることと同義だから、 P(q) = \lbrack a + b \rbrack _ q ! である。

計算方法  2

まずは  1, 2, \dots, a が順列のどこにあるかを決める。 1, 2, \dots, a 内の並び順はまだ固定しないでおく。 選択肢は  \binom{a + b}{a} 通りあるが、それぞれについて既に確定した転倒の数を数えるとその母関数は定義より  I(q) と一致する。 次に  1, 2, \dots, a の並びと、 a + 1, a + 2, \dots, a + b の並びをそれぞれ決定する。 その結果新たに確定する転倒の数の母関数は計算方法  1 と同様に考えればそれぞれ  \lbrack a \rbrack _ q !, \lbrack b \rbrack _ q ! である。 よって  P(q) = I(q) \lbrack a \rbrack _ q ! \lbrack b \rbrack _ q !

 以上の結果を合わせると  \lbrack a + b \rbrack _ q ! = I(q) \lbrack a \rbrack _ q ! \lbrack b \rbrack _ q ! となり、 I(q) = \binom{a + b}{a} _ q が従う。

例題

cp-unspoiler

線形漸化式を持つ数列の連続する項を高速に計算する

概要

 d 項間の線形漸化式を持つ数列の第  k 項から  k + d 項までを  \operatorname{O} \left( d \log \left( d \right) \log \left( k \right) \right) の時間計算量で計算するアルゴリズムを説明する。 繰り返し二乗法と多項式の余り付き除算を用いる方法と比べると、定数倍が良いと思われる。

問題設定

 P \left( x \right), Q \left( x \right) \in \mathbb{K} \left\lbrack x \right\rbrack \left\lbrack x ^ 0 \right\rbrack Q \left( x \right) \neq 0, \operatorname{deg} \left( Q \right) = d, \operatorname{deg} \left( P \right) \lt d を満たすとする。

この記事内で多項式の除算をしたら、形式的ローラン級数 \mathbb{K} \left( \left( x \right) \right) 上で考えているものとする。 これは形式的冪級数 \mathbb{K} \left\lbrack \left\lbrack x \right\rbrack \right\rbrack と概ね同じだが、 c x ^ {- 2} のように負の指数の項も有限個持ってよいとするものである。

アルゴリズムの大枠

maspy による記事 [多項式・形式的べき級数](3)線形漸化式と形式的べき級数 | maspyのHP の議論から、以下の事実が分かる。

  •  P \left( x \right) / Q \left( x \right) x ^ k 以降の係数を並べた数列の母関数は、 d 次未満の多項式  S を用いて  S \left( x \right) / Q \left( x \right) と表される。
  •  S \left( x \right) = x ^ {- k} P \left( x \right) \pmod {Q \left( x \right)} である。

これを用いて、以下のように計算する。

  1.    1 / Q \left( x \right) を母関数とする数列の第  \left\lbrack k, k + d \right) 項を計算する (後述)。
  2.   これに  Q \left( x \right) を掛けることで、 P \left( x \right) = 1 の場合の  S \left( x \right) が求まる。これは  x ^ {- k} \pmod {Q \left( x \right)} である。
  3.    S \left( x \right) = x ^ {- k} P \left( x \right) \pmod {Q \left( x \right)} より  S \left( x \right) が分かる。
  4.   FPS の逆元を計算するアルゴリズムを用いて  S \left( x \right) / Q \left( x \right) を計算すれば、元の数列の第  k 項以降を好きなだけ計算できる。

 1 / Q \left( x \right) x ^ k から  x ^ {k + d - 1} までの係数の計算

 n = k + d - 1 とすると、 \left( x ^ n Q \left( x \right) \right) ^ {- 1} x ^ {- d + 1} から  x ^ 0 までの係数を求めればよい。

 n = 0 の場合

 \left\lbrack x ^ 0 \right\rbrack Q \left( x \right) \neq 0 より  \left\lbrack x ^ 0 \right\rbrack \left( x ^ n Q \left( x \right) \right) ^ {- 1} = \left( \left\lbrack x ^ 0 \right\rbrack Q \left( x \right) \right) ^ {- 1} であり、それ以外は  0 となる。

 n \gt 0, n = 2m の場合

 V \left( x ^ 2 \right) = Q \left( x \right) Q \left( - x \right) と置き、再帰的に  \left( x ^ m V \left( x \right) \right) ^ {- 1} x ^ {- d + 1} から  x ^ 0 までの係数を求める。 すると  \left( x ^ {2m} V \left( x ^ 2 \right) \right) ^ {- 1} x ^ {- 2d + 2} から  x ^ {0} までの係数が求まる。 奇数次の係数は必ず  0 になるから、 x ^ {- 2d + 1} から  x ^ {0} までの係数が求まったことになる。 ここに  Q \left( - x \right) を掛けることで、 Q \left( - x \right) \left( x ^ {2m} Q \left( x \right) Q \left( - x \right) \right) ^ {- 1} = \left( x ^ n Q \left( x \right) \right) ^ {- 1} x ^ {- d + 1} から  x ^ 0 までの係数が求まる。

 n \gt 0, n = 2m + 1 の場合

 V \left( x ^ 2 \right) = Q \left( x \right) Q \left( - x \right) と置き、再帰的に  \left( x ^ m V \left( x \right) \right) ^ {- 1} x ^ {- d + 1} から  x ^ 0 までの係数を求める。 すると  \left( x ^ {2m + 1} V \left( x ^ 2 \right) \right) ^ {- 1} x ^ {- 2d + 1} から  x ^ {- 1} までの係数が求まる。 偶数次の係数は必ず  0 になるから、 x ^ {- 2d + 1} から  x ^ {0} までの係数が求まったことになる。 ここに  Q \left( - x \right) を掛けることで、 Q \left( - x \right) \left( x ^ {2m + 1} Q \left( x \right) Q \left( - x \right) \right) ^ {- 1} = \left( x ^ n Q \left( x \right) \right) ^ {- 1} x ^ {- d + 1} から  x ^ 0 までの係数が求まる。


これで、時間計算量が  \operatorname{\Theta} \left( d \log \left( d \right) \log \left( k + d \right) \right)アルゴリズムが得られた。  n \leq d になった時点で再帰を打ち切り FPS の逆元を計算するアルゴリズムを適用するようにすれば、 \operatorname{\Theta} \left( d \log \left( d \right) \log \left( k / d + 2 \right) \right) \subseteq \operatorname{O} \left( d \log \left( d \right) \log \left( k \right) \right) とすることもできる。

また、詳細は省略するが、アルゴリズム中の様々な畳み込みは DFT の性質を利用して定数倍高速化できる。 例えば、 Q \left( x \right) の DFT が分かっているとき  Q \left( - x \right) の DFT を  \operatorname{\Theta}(d) の時間計算量で計算することができる。

本記事のアルゴリズムと関連するアルゴリズム群: メモ: Bostan-Mori 法で計算できるものまとめ - noshi91のメモ

中心からずれた二項係数の母関数

 \displaystyle F _ k := \sum _ {i = 0} ^ {\infty} \binom{k + 2i}{i} x ^ i = \sum _ {i = 0} ^ {\infty} \binom{k + 2i}{k + i} x ^ i = \sum _ {i = 0} ^ {\infty} \frac{(k + 2i)!}{i!(k + i)!} x ^ i

と定義する。

  •  k = 0

 k = 0 は中心二項係数と呼ばれていて、 \displaystyle F _ 0 = \frac{1}{\sqrt{1 - 4x} } である。これは天下りに与えられたものをテイラー展開で確認する以外の導出方法を知らない。

  •  k = 1

パスカルの三角形を思い浮かべると、 F _ 1 2 つ足し合わせて  F _ 0 が作れることに気付く。すなわち  F _ 0 = 1 + 2x F _ 1 であり、ここから  \displaystyle F _ 1 = \frac{1 - \sqrt{1 - 4x} }{2x\sqrt{1 - 4x}} である。

  •  k が一般の場合

 k = 1 の場合と似たようにして、 F _ k F _ {k + 2} を足し合わせることで  F _ {k + 1} が作れる。すなわち  F _ {k + 1} = F _ k + xF _ {k + 2} が成り立つ。これは  3 項間漸化式なので特性方程式  f = 1 + xf ^ 2 を解くと  \displaystyle f = \frac{1 \pm \sqrt{1 - 4x} }{2x} となる。あとは  F_ 0, F _ 1 から係数を決定すれば一般項が分かるのだが、実は以下のようになる。

\begin{align} F _ k = \frac{1}{\sqrt{1 - 4x} } \left( \frac{1 - \sqrt{1 - 4x} }{2x} \right) ^ k \end{align}