Generating All the Minimal Separators of a Graph
問題
グラフ $G = (V, E)$ の頂点部分集合 $S \subseteq V$ が separator であるとは $G \setminus S$ が非連結となることをいう．包含関係の意味で極小の separator を minimal separator という．minimal separator を列挙せよ．
理論
Minimal separator は以下の生成的特徴づけをもつ．
Theorem 1. Minimal separator 全体 $\mathcal{S}$ は以下の手続きで生成できる．
- (Base Case) $x \in V$ をとる，$G \setminus (\{x\} \cup N(x))$ の連結成分 $C$ について $N(C) \in \mathcal{S}$．
- (Induction Case) $X \in \mathcal{S}$, $x \in X$ をとる．$G \setminus (X \cup N(x))$ の連結成分 $C$ について $N(C) \in \mathcal{S}$.
この特徴づけをそのままアルゴリズムとして実装すればよい．
実装
- 頂点部分集合に対する操作が多発するため BitSet を利用する．集合サイズが 64bit に収まることは仮定しない．
- すでに集合を出力したかどうかを判定するデータ構造が必要になる．ここでは素朴に HashSet を使う．頂点を文字と思って Trie を使うことで最悪計算量評価を得ることもできる．
頂点集合 $S$ で切ったときの連結成分を高速に計算できる必要がある．
for u in V:
W = V - N(u) - u;
for v in V:
C = bfs_on(v, W)
list.push(N(C))
for S in
代表元を変更できる Union Find
Union-Find データ構造では各頂点が自分の親を覚えて根を木の代表元に採用します．よくある実装では根の親ポインタは特に指すものがないので None にしたりグループのサイズを覚えたりしますが．ここに代表元を登録するようにすれば代表元を任意に選択できる Union-Find になります．
enum Link {
Parent(usize),
Representative(usize),
}
struct UnionFind {
link: Vec<Link>,
}
impl UnionFind {
fn set_representative(&mut self, item: usize, repr: usize) {
let root = self._find(item);
self.link[root] = Representative(repr);
}
fn _find(&mut self, mut item: usize) -> usize {
while let Link::Parent(parent) = self.link[item] {
item = if let Link::Parent(grand_parent) = self.link[parent] {
self.link[item] = Parent(grand_parent);
grand_parent
} else {
parent
}
}
item
}
fn find(&mut self, item: usize) -> usize {
let root = self._find(item);
if let Link::Representative(repr) = self.link[root] {
repr
} else {
unreachable!()
}
}
}
使い道
最大マッチングを求める Edmonds' Blossom Algorithm で使います．アルゴリズムでは奇サイクル (blossom) を見つけて縮約していきます．その際に最も根本の頂点を代表限に採用したくなるのでこのアルゴリズムを利用できます．
ネットワーク単体法の実装
ネットワーク単体法 (network simplex algorithm) は最小費用流問題に対する実用的に高速なアルゴリズムです．この記事ではそれなりに高速なネットワーク単体法をフルスクラッチで実装するための技術を紹介します．
ネットワーク単体法の概要
$G = (V, A)$ を有向グラフ，$c: A \to \mathbb{R}_{\ge 0}$ を枝容量，$w: A \to \mathbb{R}$ を枝費用とします．フロー $f: A \to \mathbb{R}$ について残余ネットワーク (residual network) を「フローをまだ流しうる枝からなるグラフ」として定義します．数式で書けば $G_f = (V, E_f)$ ただし $E_f = \{ a \in A : f(a) < c(a) \} \cup \{ \mathrm{rev}(a) \in A : f(e) > 0 \}$ です．ここで $\mathrm{rev}(a)$ は $a$ の逆枝をあらわします．
残余ネットワーク上の閉路 $C$ に沿ってフローを流すことで現在のフロー $f$ を更新できます．もし閉路 $C$ の枝費用和（ただし逆枝の費用は元枝の符号を変えたものとする）が負であればこの更新によってより費用の小さなフローが得られます（このような閉路を負閉路 (negative cycle) といいます）．逆に，残余ネットワークが負閉路をもたないとき $f$ が最適解であることが示せます．この性質に基づき，残余ネットワークの負閉路を消去し続けるアルゴリズムを負閉路消去アルゴリズム (negative cycle canel algorithm) といいます．ネットワーク単体法は負閉路消去アルゴリズムの一種です．
ネットワーク単体法では負閉路を効率的に発見するために全域木に注目します．フロー $ f $ が全域木解 (spanning tree soution) であるとは，ある（枝の向きを無視した）全域木 $T$ が存在して全ての $a \in A \setminus T$ の容量がタイト (tight)，すなわち $f(a) \in \{ 0, c(a) \}$ が成立することをいいます．全域木解について以下の結果が知られています．
- 最適解である全域木解が存在する．
- 最適でない全域木解についてある枝 $a \in A \setminus T$ が存在して $T$ 上で $e$ が誘導する閉路が負閉路をなす．
全域木解 $f$ と $a \in A \setminus T$ について，$a$ が誘導する閉路上にフローを流せるだけ流して得られる解もまた全域木解となります．これは閉路上にフローを流せるだけ流すということは閉路上のいずれかの枝 $a'$ の容量いっぱいまでフローを流すということなので，その枝を全域木から外すことができて，新しいフローに対応する全域木として $T + a - a'$ がとれることからわかります．
ネットワーク単体法は全域木解 $f$ および対応する全域木 $T$ を陽に管理し，負閉路を誘導する枝 $a$ を発見し，負閉路上にフローを流せるだけ流す，という手続きを繰り返すアルゴリズムです．擬似コードで書くと次のようになります．
while find_entering_arc():
pivot()
この擬似コードからも明らかですが，ネットワーク単体法の実装では以下の２点が重要となります．
- どうやって負閉路を誘導する枝 $e$ を見つけるか
- どうやって全域木を更新するか
負閉路を誘導する枝の見つけ方
以下では全域木を根付き木の形でもっているとします（詳細は後述）．
枝 $a = (u, v)$ を全域木に追加したときの閉路の費用を効率的に計算するために各頂点 $u$ に根から $u$ までの木に沿った費用和 $\pi(u)$ を保存しておきます．すると $(u, v)$ がなす閉路の費用は $w_{\pi}(a) = \pi(u) - \pi(v) + w(a)$ となります．$w_{\pi}(a)$ は簡約費用 (reduced cost) と呼ばれています．
簡約費用が負の枝を見つけて更新する手続きを繰り返せば有限回の手続きで最適解に到達するので，どの枝を選ぶかに設計の余地があります．実用的に高速と言われている基準がblock pivot ruleです．この基準では枝を適当なサイズ（典型的には $\sqrt{m}$）のブロックに分けてその中で最小のものを選びます．ただしブロック内に簡約費用が負のものがなかった場合は探索範囲を拡張します．
全域木の更新
全域木に必要な操作は以下の通りです．
- 子孫を走査できる（∵ $\pi$ を更新するため）．
- 親を見つけられる（∵ 新しく容量いっぱいになる枝 $e'$ を見つけるため）
- 枝を追加・削除できる．
これらの操作を見ると動的木を使うのが自然に見えますね．現行最速のアルゴリズムは素朴に探索します．
extended threaded index (XTI) と呼ばれる技術を用います．XTI は頂点に対して以下の情報を覚えます．
- thread_next: 深さ優先探索順で次の頂点．最後の頂点からは根に戻る．
- thread_prev: 上の逆．
- size: 自分以下の部分木の頂点数
- last: 自分以下の最後の子孫
子孫の走査は thread_next を last まで手繰ればできます．
頂点 x の部分木を削除するには
- thread_ next(thread_prev(x)) = thread_next(last(x))
- thread_next(last(x)) = x
でできます．
サイズの更新．
- v_out の ancestor に size(x) を足す
- v_in の ancestor から size(x) を引く
実装の解説
以下，説明のための擬似コードのかわりに Rust のコード片を貼ります．頂点と枝は usize で管理しており，必要なデータは配列 nodes[u], arcs[a] に保存されています．実際に動くコードは (TODO) にあります．
fn find_entering_arc(&mut self) -> Option<usize> { let mut min_cost = C::zero(); let mut min_arc = None; let arcs = (self.next_arc..self.arc_num).chain(0..self.next_arc); for (count, a) in arcs.enumerate() { if count >= self.block_size && min_arc.is_some() { self.next_arc = a; break; } let cost = self.arcs[a].state.sign() * self.residual_cost(a); if cost < min_cost { min_cost = cost; min_arc = Some(a); } } min_arc }
枝 $a_{in} = (u_{in}, v_{in})$ を追加したときにできる閉路から抜ける枝を探します．場合分けを減らすためにフローは $u_{in}$ から $v_{in}$ の向きに流れるものとしておきます．
let u_in = self.arcs[in_arc].tail; let v_in = self.arcs[in_arc].head; // Push a flow from `source` to `target`. let (source, target) = match self.arcs[in_arc].state { Lower => (u_in, v_in), Upper => (v_in, u_in), _ => unreachable!(), };
実際に閉路を探す手続きは $u_{in}$ と $v_{in}$ の最小共通先祖を探す手続きとほぼ同一です．各頂点の部分木の頂点数を比較して小さいほう（i.e., 深いほう） を親方向に進みます．探索中に出会った残余容量最小の枝を記憶しておきます．
let (join, delta, u_out, source_side) = { let mut delta = self.arcs[in_arc].upper; let mut u_out = usize::MAX; let mut source_side = true; let mut u = source; let mut v = target; while u != v { if self.nodes[u].size < self.nodes[v].size { let residual_capacity = match self.nodes[u].pred { Up(a) => self.arcs[a].flow, Down(a) => self.arcs[a].upper - self.arcs[a].flow, }; if residual_capacity < delta { delta = residual_capacity; u_out = u; source_side = true; } u = self.nodes[u].parent; } else { let residual_capacity = match self.nodes[v].pred { Up(a) => self.arcs[a].upper - self.arcs[a].flow, Down(a) => self.arcs[a].flow, }; if residual_capacity <= delta { delta = residual_capacity; u_out = v; source_side = false; } v = self.nodes[v].parent; } } (u, delta, u_out, source_side) };
続いてサイクルに沿ってフローを流せるだけ流します．
if delta.is_positive() { match self.arcs[in_arc].state { Lower => self.arcs[in_arc].flow += delta, Upper => self.arcs[in_arc].flow -= delta, _ => unreachable!(), }; let mut u = source; while u != join { self.arcs[self.nodes[u].pred.unwrap()].flow -= self.nodes[u].pred.sign() * delta; u = self.nodes[u].parent; } let mut v = target; while v != join { self.arcs[self.nodes[v].pred.unwrap()].flow += self.nodes[v].pred.sign() * delta; v = self.nodes[v].parent; } }
ここから全域木の更新に入ります．
一番簡単な場合として，追加した枝 $a_{in}$ 自身がタイトになることがあります．これは $u_{in}$ から $v_{in}$ にフローをもっと流せたという場合です．このときは単純に $a_{in}$ の状態を更新して終了します．
if u_out == usize::MAX { self.arcs[in_arc].state = self.arcs[in_arc].state.flip(); return; }
枝 $a_{out} = (u_{out}, v_{out})$ を抜ける枝とします．
簡単なケースとして $u_{in} = u_{out}$ の場合があります．
$u_{in}$ と $u_{out}$ が同じ側にあるようにしておきます．
頂点 u の子孫を走査するのは簡単です．
end = last(u)
do {
# do something
u = next(u)
} while u != end
頂点 x 以下の部分木の切除
# x を含む部分木の更新 next(last(x)) = x # x を含まない部分木の更新 next(prev(x)) = next(last(x))
その他
ネットワーク単体法では縮退 (degeneracy) とよばれる現象が多発します．
枝をランダムな順序に並べ替えておくと高速化します．
Least Common Ancestor O(d) time and O(1) space
木の頂点 $u$, $v$ について最小共通先祖 (LCA) を計算するシンプルなアルゴリズムを紹介します．複数の頂点の組に対して LCA を求めるのであれば，木を前処理をしてクエリあたり $O(\log n)$ なり $O(1)$ で求めるのが定石ですが，ここではそうでなく，線形時間かけて求めることを考えます．木は親ポインタ parent: &[Option<usize>] によって与えられるとします．
O(n) 時間 O(n) 空間アルゴリズム
すぐに思いつくアルゴリズムは頂点に訪問済みかどうかをあらわすフラグをもたせて各頂点から親方向に走査するものでしょう．このアルゴリズムは領域の確保に $O(n)$ 時間かかる上に $O(n)$ 空間を消費するので非効率的です．
fn lca(parent: &[Option<usize>], u: usize, v: usize) -> usize { let mut x = u; let mut y = v; let mut seen = vec![false; parent.len()]; loop { seen[x] = true; if let Some(p) = parent[x] { x = p; } else { break; } } loop { if seen[y] { return y; } if let Some(p) = parent[y] { y = p; } else { unreachable!(); } } }
O(d) 時間 O(1) アルゴリズム ver 1
別のアルゴリズムとして $O(d)$ 時間かけて $u$ と $v$ の深さを求め，深い方だけを親方向に走査して深さを揃えてから同時に親に進んでいくものがあります．このアルゴリズムは全体で $O(d)$ 時間 $O(1)$ 空間を達成するので最適です．
fn depth(parent: &[Option<usize>], u: usize) -> usize { let mut x = u; for d in 0.. { if let Some(p) = parent[x] { x = p; } else { return d; } } unreachable!() } fn lca(parent: &[Option<usize>], u: usize, v: usize) -> usize { let mut x = u; let mut y = v; let d_x = depth(parent, x); let d_y = depth(parent, y); for _ in d_y..d_x { x = parent[x].unwrap(); } for _ in d_x..d_y { y = parent[y].unwrap(); } while x != y { x = parent[x].unwrap(); y = parent[y].unwrap(); } x }
O(d) 時間 O(1) 空間アルゴリズム ver 2
上述のアルゴリズムは最適ですが，こんな簡単な手続きのためにこんなコード量を使いたくないところです．そこで以下では別の最適なアルゴリズムを紹介します．
fn lca(parent: &[Option<usize>], u: usize, v: usize) -> usize { let mut x = u; let mut y = v; while x != y { x = parent[x].unwrap_or(v); y = parent[y].unwrap_or(u); } x }
このアルゴリズムでは，変数 $x$ は頂点 $u$ から親方向に進み，根に到達したら頂点 $v$ に戻って同じ手続きを繰り返します．最初に根に至るまでを１巡目，次に根に至るまでを２巡目と呼ぶことにします．変数 $y$ は頂点 $v$ から同様の手続きを行います．以下このアルゴリズムが LCA を求めることを証明します．
$a$ を $u$ と $v$ の LCA とします．もし $u$ と $v$ の深さが等しければ１巡目の走査によって $a$ が発見されます．$u$ と $v$ の深さが異なるとしましょう．このとき変数 $x$ が頂点 $a$ に到達するステップ数は となります．同様に変数 $y$ が頂点 $a$ に到達するステップ数は
となります．これらは等しいのでそれぞれ２巡目の走査で $a$ が発見されます．
Modulo Square Root
正整数 $ m $ と整数 $ a \in [0, m) $ について $ x^2 = a \mod m $ の解を求めます．
方針
$ m = p_1^{e_1} \cdots p_k^{e_k} $ と素因数分解し，各素数べきについて $ x^2 = a \mod p_i^{e_i} $ を解いて，中国剰余定理を用いて解を構成します．
素因数分解と中国剰余定理部分は解説しません．
前処理
$ x^2 = 0 \mod p^e $ の場合，$ x^2 = \alpha p^e$ と書いて両辺の $ p $ 因子の数を比較すると $x = b p^{\lceil e / 2 \rceil}$ と書けることがわかります．この形の数はすべて方程式を満たすのでこれらが解となります．解の個数は $ \sqrt{p^e} $ です．
$ x^2 = a \mod p^e $ ($ a \neq 0 $) の場合，$ a = b p^k $ と $ p $ 因子を括り出して $ x^2 = b p^k + \alpha p^e$ と書き，両辺の $ p $ 因子の数を比較すると $ k $ が奇数のとき解が存在しないこと・$ k $ が偶数のとき $ x = y p^{k / 2} $ とおけることがわかります．したがって以下 $y^2 = b \mod p^e$ ただし $ \mathrm{gcd}(b, p) = 1 $ について考えれば十分です．
Hensel Lift
一般に $ y^2 = b \mod p^e $ を満たす $y $ は任意の $ 0 \le k \le e $ について $ y^2 = b \mod p^k $ を満たします．このことを手がかりに，まずは低い次数の素数べきでの解を求め，それを高い次数の素数べきに「持ち上げる」ことを考えます．
$l \le k$ とします．$ y_k^2 = b \mod p^k $ を満たす $ y_k $ をとり，$ y_{k+l}^2 = b \mod p^{k+l} $ を満たす $ y_{k+l} $ のうち $ y_{k+l} = y_k + t p^k $ と書けるものを探します．直接計算により
$$ y_{k+l}^2 = y_k^2 + 2 y_k t p^k + t^2 p^{2 k} = b + \left( \frac{y_k^2 - b}{p^k} + 2 y_k t \right) p^k \mod p^{k+l} $$
となるので $ (y_k^2 - b)/p^k + 2 y_k t = 0 \mod p^{l} $ なる $ t $ を求めればよいことがわかります；式中に現れる分数は割り切れることに注意してください．この構成（を一般化したもの）は Hensel Lift と呼ばれています．$l \approx k$ にとることで $O(\log e)$ 回の Hensel Lift で解を構成できます．
以下 $ p $ が奇素数の場合と $ 2 $ の場合にわけて具体的な手続きを与えます．
奇素数べき
$ y_1^2 = b \mod p $ を満たす $y_1$ は Tonelli-Shanks や Cipolla のアルゴリズムによって見つけられます．この解を初期値として Hensel Lift を行います．
$ y_k^2 = b \mod p^k $ なる $y_k$ をもっているとします．もし $y_k$ が $ p $ と互いに素であれば $2 y_k$ は法 $p^l$ のもとでモジュロ逆元 $(2 y_k)^{-1}$ をもち，これを用いて表せる $ t = ((y_k^2 - b)/p^k) \times (2 y_k)^{-1} $ が Hensel Lift の条件を満たす $ t $ となります．したがって $ y_{k+l} = y_k + t p^k = y_k + (y_k^2 - b) (2 y_k)^{-1} \mod p^{k+l} $ に解を持ち上げることができました．この式から $ y_{k+l} $ もまた $ p $ と互いに素であるのでこの手続きを繰り返すことができ，任意次数まで解を持ち上げることができます．
2 べき
$ k \le 2 $ のときは全探索で解が求まります．具体的には以下がわかります：
- $ y^2 = b \mod 2 $ は前処理より $b = 1$，解は $x = 1$．
- $ y^2 = b \mod 4 $ は前処理より $b = 1$ または $b = 3$．全通り試すことで $b = 1$ の場合にのみ解をもち，解は $y = 1$ と $y = 3$．
$ k \ge 3 $ を考えます．$ s = (y_k^2 - b)/2^k $ とおいて Hensel Lift の式を眺めると $ s + 2 y_k t = 0 \mod 2^l $ であればよいとわかります．この式から $ s $ は 2 の倍数でなければならないので $ s' = s / 2 $ とおけて $ 2 (s' + y_k t) = 0 \mod 2^l $ となり，$l = 1$ では $t = 0$, $ = 1$ がともに条件を満たし，$l \ge 2$ では $ s' + y_k t = 0 \mod 2^l $ および $ s' + y_k t = 2^{k - 1} \mod 2^l $ に帰着します．いずれも $ y_k $ は $ 2 $ と互いに素なのでモジュロ逆元がとれて $ t $ が計算できます．
なお，この構成から $k \ge 3$ のとき解が 0 個または 4 個であることがわかります．
実装例
ソースコードの全体はテンプレート周りが長いので一部のみを貼ります．
pub trait SqrtMod: Int { fn sqrt_mod_all(self, modulo: Self) -> Vec<Self> { let solve = |p: Self, e: u32| -> Vec<Self> { let p_e = p.clone().pow(e); let a = self.clone() % p_e.clone(); return if a.is_zero() { let mut solutions = Vec::new(); let mut curr = Self::zero(); let step = p.pow((e + 1) / 2); while curr < p_e { solutions.push(curr.clone()); curr = curr + step.clone(); } solutions } else { let (b, l) = a.trim_factor(p.clone()); if l.is_odd() { Vec::new() } else { let mut solutions = Vec::new(); for x in b.sqrt_mod_prime_power(p.clone(), e - l) { let begin = x * p.clone().pow(l / 2); let step = p.clone().pow(e - l / 2); let mut curr = begin.clone(); loop { solutions.push(curr.clone()); curr = curr.add_mod(step.clone(), p_e.clone()); if curr == begin { break; } } } solutions.sort_unstable(); solutions.dedup(); solutions } }; }; let mut modulos = Vec::new(); let mut solutions = Vec::new(); for (p, e) in modulo.prime_factors() { solutions.push(solve(p.clone(), e)); modulos.push(p.pow(e)); } solutions .into_iter() .map(|vec| vec.into_iter()) .multi_cartesian_product() .map(move |data| { let (a, _) = chinese_remainder_theorem(data.into_iter().zip(modulos.clone())).unwrap(); a }) .collect() } fn sqrt_mod_two_power(self, e: u32) -> Vec<Self> { let two = Self::one() + Self::one(); let four = two.clone() + two.clone(); let eight = four.clone() + four.clone(); if e == 1 { return vec![Self::one()]; } if e == 2 { if !(self % four).is_one() { return Vec::new(); } else { return vec![Self::one(), Self::one() + Self::two()]; } } if !(self.clone() % eight).is_one() { return Vec::new(); } let mut ys = vec![Self::one(), Self::one() + Self::two()]; let mut p_k = four; let mut k = 2; while k < e { let l = k.min(e - k); let p_l = if k == l { p_k.clone() } else { two.clone().pow(l) }; let p_kl = p_k.clone() * p_l.clone(); let mut next = Vec::new(); for y in ys { let s = (y.clone() * y.clone() - self.clone()) / p_k.clone(); if s.is_odd() { continue; } if l == 1 { next.push(y.clone()); next.push(y.clone() + p_k.clone()); } else { let s = (s / two.clone()) % p_l.clone(); let inv_y = (y.clone() % p_l.clone()).inv_mod(p_l.clone()).unwrap(); let t = (-s.clone() * inv_y.clone()) % p_l.clone(); let new_y = (y.clone() + t * p_k.clone()) % p_kl.clone(); next.push(if new_y.is_negative() { p_kl.clone() + new_y } else { new_y }); let t = ((p_l.clone() / two.clone() - s) * inv_y.clone()) % p_l.clone(); let new_y = (y.clone() + t * p_k.clone()) % p_kl.clone(); next.push(if new_y.is_negative() { p_kl.clone() + new_y } else { new_y }); } } ys = next; k += l; p_k = p_kl; } ys } fn sqrt_mod_odd_prime_power(self, p: Self, e: u32) -> Vec<Self> { if legendre_symbol(self.clone(), p.clone()) != 1 { return Vec::new(); } let mut y = (self.clone() % p.clone()) .sqrt_mod_tonelli_shanks(p.clone()) .unwrap(); let mut p_k = p.clone(); let mut k = 1; while k < e { let l = k.min(e - k); let p_l = if k == l { p_k.clone() } else { p.clone().pow(l) }; let p_kl = p_k.clone() * p_l.clone(); let two_y = (y.clone() + y.clone()) % p_l.clone(); let inv_df = two_y.inv_mod(p_l.clone()).unwrap(); let f = (y.clone() * y.clone() - self.clone()) % p_kl.clone(); y = (y - f * inv_df) % p_kl.clone(); if y.is_negative() { y = y + p_kl.clone(); } k += l; p_k = p_kl; } vec![y.clone(), p_k - y] } fn sqrt_mod_prime_power(self, p: Self, e: u32) -> Vec<Self> { if p.is_even() { self.sqrt_mod_two_power(e) } else { self.sqrt_mod_odd_prime_power(p, e) } } fn sqrt_mod(self, modulo: Self) -> Option<Self> { self.sqrt_mod_all(modulo).pop() } } impl<T: Int + Signed> SqrtMod for T {}
Sorensen's Incremental Prime Sieve
Sorensen's Incremental Sieve は素数を順番に列挙するアルゴリズムです．このアルゴリズムは incremental, compact という 2 つの性質をもっています．incremental というのは $n$ までの素数を列挙したとき $n+1$ に対する処理が高速（$ t(n) $ 回の算術演算）であること，compact というのは $n$ まで素数を列挙するのに $O^*(\sqrt{n})$ のメモリしか使わないことを指します．Sorensen's Incremental Sieve は $O(\log n / \log \log n)$-incremental time，および $O(\sqrt{n} \log n)$-bit spaceのアルゴリズムです．
方針
素数をバケットのリストで管理します．次に処理する数 $n$ に対して $\sqrt{n}$ 未満の素数 $p$ を列挙し，$(-n) \mod p$ 番目のバケットに挿入します．このとき 0 番目のバケットが空であることと $n$ が素数または素数の平方数であることが同値となるので素数判定が $O(1)$ 回の算術演算でできます．また，数を $n$ から $n+1$ に変化させたときのバケットたちの変化は簡単に管理できます．
バケットたちの変化を具体例で見ましょう．$n = 26$ とします．このとき考える素数は $p = 2, 3, 5$ であり，それらによる $n$ の剰余は $0$, $2$, $1$ です．したがってそれらを 0 番目, 1 番目, 4 番目のバケットに挿入することで (2) (3) () () (5) が得られます．ここから順にカウントアップしていくとバケットたちの様子は次のようになります．
- 26: (2) (3) () () (5)
- 27: (3) (2) () (5) ()
- 28: (2) () (3,5) () ()
- 29: () (2,3,5) () () ()
- 30: (2,3,5) () () () ()
- 31: () (2) (3) () (5)
- 32: (2) (3) () (5) ()
- 33: (3) (2) (5) () ()
これを見ればわかる通り，「各ステップで全ての数はバケットを右から左に 1 つずつ移動する．ただし 0 番目のバケットの数は $p - 1$ に移動する」という法則で変化しています．あとはこれを実装するだけです． すべてのバケットを同時に移動させるのは添字にオフセットを加えればよいので，普通のバケットの配列を用いて実装できます．
実装例
#[derive(Debug)] pub struct IncrementalSieveSorensen { n: usize, s: usize, r: usize, stacks: Vec<Vec<usize>>, pos: usize, } impl IncrementalSieveSorensen { pub fn new(start: u64) -> Self { let start = start.max(2) as usize; let r = (start as f64).sqrt() as usize + 1; let mut is_prime = vec![true; r]; let mut stacks = vec![vec![]; r + 2]; for p in 2..r { if is_prime[p] { let j = (p - start % p) % p; stacks[j].push(p); for q in (p * p..r).step_by(p) { is_prime[q] = false; } } } Self { n: start, s: r * r, r, stacks, pos: 0, } } fn is_next_prime(&mut self) -> bool { let m = self.stacks.len(); let mut is_prime = true; while let Some(p) = self.stacks[self.pos].pop() { self.stacks[(p + self.pos) % m].push(p); is_prime = false; } if self.n == self.s { if is_prime { self.stacks[(self.r + self.pos) % m].push(self.r); is_prime = false; } self.r += 1; self.s = self.r * self.r; } self.n += 1; self.pos += 1; if self.pos == m { self.pos = 0; self.stacks.push(vec![]); self.stacks.push(vec![]); } is_prime } } impl Iterator for IncrementalSieveSorensen { type Item = u64; fn next(&mut self) -> Option<u64> { while !self.is_next_prime() {} Some(self.n as u64 - 1) } }
Lenstra-Lenstra-Lovasz's Lattice Basis Reduction
概要
整数ベクトル $B := (b_1, \ldots, b_m) \in \mathbb{Z}^{m \times n} $ が与えられたとき，これらの整数線型結合で作られる集合 $L(B) := \{ x_1 b_1 + \ldots + x_m b_m : x_1, \ldots, x_m \in \mathbb{Z} \} \subseteq \mathbb{Z}^n $ を lattice といいます．Lattice $ L $ を生成する $ B $ を $ L $ の基底といいます．一般に $ L $ に対して複数の基底が存在します．例えば $ ( (1, 0), (0,1) ) $ と $ ( (1,2), (3,5) ) $ はどちらも $\mathbb{Z}^2$ を生成します．
Lattice における基本的な問題はできるだけシンプルな基底を求めることです．特に以下の問題たちが基本的です．
- 最短ベクトル問題．$L(B)$ のうち，ユークリッドノルムが最も小さいベクトルを求めよ．
- 最短基底問題．$L(B)$ の基底 $B'$ のうち，ユークリッドノルムが（適当な意味で）最も小さなものを求めよ．
Lenstra-Lenstra-Lovasz (1982): https://doi.org/10.1007/BF01457454 の Lattice-Basis Reduction Algorithm (LLL) はこれらの問題の近似解を求める多項式時間アルゴリズムです．近似率は $2^{\Omega(m)}$ とあまりよろしくないのですが，対象となるベクトルの本数が少なければ十分有用な解を得ることができます．
LLL の応用
無理数の有理数近似
円周率 $\pi = 3.14159265...$ を近似する有理数として $22/7 = 3.1428...$ や $355/113 = 3.14159292...$ が知られています．これらの数値を求める手法として連分数展開が知られていますが，LLL を用いて求めることもできます．
十分大きな整数 $ M $ を固定してベクトル $b_1 = (1, 0, \lfloor M \pi \rfloor)$, $b_2 = (0, 1, -M)$ を考えます．これらが生成する lattice から $(a, b, s)$ が取れたとしましょう．これは $a \lfloor M \pi \rfloor - b M = s$ を意味するため，$ a \pi - b \approx 0$ となり，$\pi$ の近似値として $b / a$ が得られました．
無理数を解にもつ代数方程式
黄金比 $ \gamma = 1.61803398... $ は代数方程式 $1 + x - x^2 = 0$ を満たすことを知っていますが，この方程式を LLL を用いて求めることができます．
十分大きな整数 $ M $ を固定してベクトル $b_1 = (1, 0, 0, M)$, $b_2 = (0, 1, 0, \lfloor M \gamma \rfloor)$, $b_3 = (0, 0, 1, \lfloor M \gamma^2 \rfloor)$ を考えます．これらが生成する lattice から $(a, b, c, s)$ が取れたとしましょう．これは $ a M + b \lfloor M \gamma \rfloor + c \lfloor M \gamma^2 \rfloor = s $ を意味するため，$ \gamma $ が $ a + b x + c x^2 \approx 0 $ を満たすことがわかりました．
小さな Bezout 係数
整数 $ a_1, \ldots, a_m $ が与えられたとき，$ g := \mathrm{gcd}(a_1, \ldots, a_m) = x_1 a_1 + \cdots + x_m a_m $ を満たす $ (x_1, \ldots, x_m) $ を Bezout 係数といいます．Bezout 係数は拡張ユークリッド互除法で求まることが知られています．Bezout 係数は一意ではありません．例えば $ \mathrm{gcd}(3141, 5926, 5358) = 1 $ について $(-1881, 997, 0)$ はひとつの Bezout 係数ですが，他にも $(81, -14, -32)$ という係数の小さな Bezout 係数が存在します．一般に，ノルム最小の Bezout 係数を求める問題は NP 困難ですが，近似でよければ LLL を用いて求められます．
記号的な無限大 $ M $ を導入して $ b_1 = (1, 0, \ldots, 0, a_1 M), \ldots, b_m = (0, \ldots, 0, 1, a_m M) $ を考えます．$ M $ は無限大なのでノルムの小さな基底において最終成分の値はできる限り小さくなるはずです．具体的にはどれか１つのベクトルだけ $g M $，残りは $0$ となります．$ g M $ があらわれるベクトルの成分は構成から Bezout 係数となり，LLL の性質から小さな値をもつことが期待できます．
整数計画問題
変数の数が定数個の整数計画問題も LLL を使うと多項式時間で解くことができます． 一般に，「格子点を含まない多面体は細長い（太さが次元にのみ依存する定数になる）」ということが証明できます (Khinchin's Flatness Theorem)．多項式時間で多面体から格子点を見つけるもしくは細長い軸を求めることができるので (*)，後者の場合にその軸に沿って多面体をスライスして分割統治を行うことで整数計画問題の実行可能性問題を解くことができます．このとき分枝の個数は次元にのみ依存する定数になり，分枝の深さは次元に等しくなるので，全体として探索木のサイズが定数となります．(*) の手続きのなかで LLL が用いられます．
理論的に価値の大きな結果ですが，実用的には通常の分枝限定法のほうが遥かに早いと言われています．
LLL の実装
LLL は基底に対して適当な基準で行基本変形を繰り返すタイプのアルゴリズムです．基準をうまく選ぶことによって手続き中でノルムが減少していくことが示せます．正当性の証明は複雑ですが，実装自体は Gauss の消去法やら拡張互除法やらとそれほど大きく違わないシンプルなものとなります．
Wikipedia.en: https://en.wikipedia.org/wiki/LLL_algorithm に擬似コードがあります．これは Lenstra-Lenstra-Lovasz によるオリジナル実装に近いものですが，有理数を管理しなければいけないため定数倍がよくありません．有理数を必要としない実装が de Weger (1987): https://doi.org/10.1016/0022-314X(87)90088-6 によって提案されています．以下はそれを実装したものです．
注意として，LLL を有効に使うには原則として多倍長演算が必要です．これは行基本変形中に係数が増大するためですが，それでも高々多項式程度に収まるのが LLL の偉いところです．
use num::traits::{NumRef, RefNum}; use num::{Integer, Signed}; pub trait FracRound { fn frac_round(&self, other: &Self) -> Self; } impl<T> FracRound for T where T: Clone + Integer + Signed + std::fmt::Debug, { fn frac_round(&self, other: &Self) -> Self { let sign = self.signum() * other.signum(); let num = self.abs(); let den = other.abs(); let (div, rem) = num.div_rem(&den); let value = match (rem.clone() + rem.clone()).cmp(&den) { std::cmp::Ordering::Less => div, std::cmp::Ordering::Equal => { if sign.is_positive() { div } else { div + T::one() } } std::cmp::Ordering::Greater => div + T::one(), }; sign * value } } pub struct LatticeBasisReductionDeWeger<T> { b: Vec<Vec<T>>, u: Vec<Vec<T>>, lambda: Vec<Vec<T>>, d: Vec<T>, } impl<T> LatticeBasisReductionDeWeger<T> where T: Clone + Integer + Signed + NumRef + From<i8> + std::fmt::Debug, for<'a> &'a T: RefNum<T>, { fn dot(a: &[T], b: &[T]) -> T { a.iter() .zip(b.iter()) .fold(T::zero(), |sum, (a, b)| sum + a * b) } pub fn new<S: Clone + Into<T>>(matrix: Vec<Vec<S>>) -> Self { let m = matrix.len(); let mut b = Vec::new(); for row in matrix { b.push(row.into_iter().map(|x| x.into()).collect::<Vec<_>>()); } let mut c = b.clone(); let mut u = vec![vec![T::zero(); m]; m]; for i in 0..m { u[i][i] = T::one(); } let mut d = vec![T::one(); m]; let mut lambda = vec![vec![T::zero(); m]; m]; for i in 0..m { for j in 0..i { lambda[i][j] = Self::dot(&b[i], &c[j]); for p in 0..c[i].len() { c[i][p] = &d[j] * &c[i][p] - &lambda[i][j] * &c[j][p]; if j >= 1 { c[i][p] = &c[i][p] / &d[j - 1]; } } } d[i] = Self::dot(&c[i], &c[i]); if i >= 1 { d[i] = &d[i] / &d[i - 1]; } } Self { b, u, lambda, d } } fn reduce(&mut self, k: usize, l: usize) { if T::from(2) * self.lambda[k][l].abs() > self.d[l] { let r = self.lambda[k][l].frac_round(&self.d[l]); for p in 0..self.b[k].len() { self.b[k][p] = &self.b[k][p] - &r * &self.b[l][p]; } for p in 0..self.u[k].len() { self.u[k][p] = &self.u[k][p] - &r * &self.u[l][p]; } for j in 0..l { self.lambda[k][j] = &self.lambda[k][j] - &r * &self.lambda[l][j]; } self.lambda[k][l] = &self.lambda[k][l] - &r * &self.d[l]; } } fn swap(&mut self, k: usize) { self.b.swap(k, k - 1); self.u.swap(k, k - 1); for j in 0..k - 1 { // swap lambda[k-1][j] and lambda[k][j] let (a, b) = self.lambda.split_at_mut(k); std::mem::swap(&mut a[k - 1][j], &mut b[0][j]); } for i in k + 1..self.lambda.len() { let t = &self.lambda[i][k - 1] * &self.d[k] - &self.lambda[i][k] * &self.lambda[k][k - 1]; self.lambda[i][k - 1] = (&self.lambda[i][k - 1] * &self.lambda[k][k - 1] + if k >= 2 { &self.lambda[i][k] * &self.d[k - 2] } else { self.lambda[i][k].clone() }) / &self.d[k - 1]; self.lambda[i][k] = t / &self.d[k - 1]; } self.d[k - 1] = (if k >= 2 { &self.d[k - 2] * &self.d[k] } else { self.d[k].clone() } + &self.lambda[k][k - 1] * &self.lambda[k][k - 1]) / &self.d[k - 1]; } pub fn run(&mut self) { let m = self.b.len(); let mut k = 1; while k < m { self.reduce(k, k - 1); let lhs = if k >= 2 { &self.d[k - 2] * &self.d[k] - &self.lambda[k][k - 1] * &self.lambda[k][k - 1] } else { &self.d[k] - &self.lambda[k][k - 1] * &self.lambda[k][k - 1] }; let rhs = &self.d[k - 1] * &self.d[k - 1]; if T::from(4) * lhs < T::from(3) * rhs { self.swap(k); if k >= 2 { k -= 1; } } else { for i in (0..k - 1).rev() { self.reduce(k, i); } k += 1; } } } }
