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; } } } }
