immut/dense の設計

設計目標

DensePolynomial[A] は値としての一変数多項式です。自由に共有でき、== で比較でき、環の演算子で組み合わせられる正規形の係数ベクトルです。係数の大半が非ゼロである多項式を対象としており、その場合は次数までのすべての係数を格納するのが最も単純で、かつ最も高速なレイアウトです。

数学的背景

$R$ 上の一変数多項式は有限台の列 $(c_0, c_1, \dots)$ であり、$f = \sum_i c_i x^i$ と書きます。次数は $\deg f = \max\{ i \mid c_i \neq 0 \}$ で、$\deg 0 = -\infty$ とします（None として返されます）。加法は係数ごとに行い、乗法はコーシー積です。

$$(fg)_k = \sum_{i + j = k} f_i \, g_j .$$

$R$ が可換環なら $R[x]$ も可換環であり、$R$ が（UInt のように）半環にすぎないなら $R[x]$ は半環です。次数は次を満たします。

$$\deg(f + g) \le \max(\deg f, \deg g), \qquad \deg(fg) \le \deg f + \deg g ,$$

2 つ目の式で等号が成り立つのは、先頭係数どうしの積が非ゼロのときちょうどであり、$R$ が零因子を持たなければ常に成り立ちます。

設計上の判断

正規形: 末尾のゼロを取り除く

問題。 列 $(1, 2)$ と $(1, 2, 0, 0)$ は同じ多項式を表します。両方を格納できるとすると、==、degree、length、ハッシュのすべてが最初に正規化を行わなければなりません。

選択。 すべてのコンストラクタとすべての演算は、返す前に末尾のゼロを取り除き、零多項式は空のベクトルとします。この切り詰めにより、多項式と $c_{n-1} \neq 0$ である語 $(c_0, \dots, c_{n-1})$（および空語）との間に全単射が定まるので、次が成り立ちます。

$$\texttt{p == q} \iff p = q \text{ in } R[x], \qquad \texttt{p.length()} = \deg p + 1 .$$

構築時だけでなく すべての 演算の後に切り詰めることが重要なのは、上の次数の不等式が狭義になりうるからです。$\mathbb{Z}/2^{32}$ である Int では $(2^{16} x + 1)^2 = 2^{32} x^2 + 2^{17} x + 1 = 2^{17} x + 1$ となり、先頭係数が消えます。切り詰めた結果は正しく次数 $1$ を持ちます。

ベクトルは永続的な @immut/vector.Vector であり、from_coefficients は入力をコピーするので、呼び出し側が後で自分の配列を変更しても多項式は変わりません。

演算ごとに最小限の係数制約

問題。 型全体に A : Ring のような単一の制約を課すと、有用な係数型が除外されてしまいます。UInt には符号反転がなく、単位元のない乗法を持つ型もあります。

選択。 各関数は、自分が使う luna-generic の能力の最小集合を明示します。加法には Eq + AddMonoid（切り詰めのための Eq と Zero）が必要で、乗法はそれに Mul を加えますが One は不要です。符号反転には Neg が必要ですが Mul は不要で、One を必要とするのは variable、one、pow、karatsuba だけです。型自体には制約がないため、zero()、length()、degree() は任意の A で動作します。

既定は筆算式の乗法

* は格納された係数に対する二重ループでコーシー積を計算し、長さ $m$ と $n$ に対して $mn$ 回の係数の乗算と加算を行ってから切り詰めます。典型的な用途の大半を占める短い多項式では、これはどの再帰的手法よりも高速で、Neg を必要とせず、厳密な係数型に対しては厳密です。

明示的なメソッドとしての Karatsuba

問題。 オペランドが長いと、$O(mn)$ がボトルネックになります。

導出。 各オペランドを $k$ で分割します。$y = x^k$、$\deg a_0, \deg b_0 < k$ として $a = a_0 + a_1 y$、$b = b_0 + b_1 y$ とおくと、

$$\begin{aligned} ab &= a_0 b_0 + (a_0 b_1 + a_1 b_0)\, y + a_1 b_1 \, y^2, \\ a_0 b_1 + a_1 b_0 &= (a_0 + a_1)(b_0 + b_1) - a_0 b_0 - a_1 b_1 , \end{aligned}$$

となります。2 行目は分配法則だけから従うので、任意の環で成り立ちます。4 回の積が半分のサイズの 3 回の積に置き換わります。長さ $n$ に対するコストを $T(n)$、加算とシフトのコストを $c\,n$ とすると、

$$T(n) = 3\,T(n/2) + c\,n \;\Rightarrow\; T(n) = c\,n \sum_{j=0}^{\log_2 n} \left(\tfrac{3}{2}\right)^j = O\!\left(n \cdot \left(\tfrac32\right)^{\log_2 n}\right) = O\!\left(n^{\log_2 3}\right) \approx O(n^{1.585}).$$

選択。 karatsuba は長いほうの長さの半分で分割し、3 つの積について再帰し、短いほうのオペランドの係数が $32$ 個以下になると * にフォールバックします。その範囲では再帰のオーバーヘッドが節約される乗算を上回るからです。減算には Neg が、$y$ によるシフト scale(k, One::one()) には One が必要です。これを * の実装にせず独立したメソッドとしているのは、* がより弱い制約と予測可能なコストを保てるようにするためです。テストで、しきい値を超える範囲で両者が一致することを確認しています。

Horner 法による評価

$p(a)$ は次のように計算します。

$$p(a) = c_0 + a\bigl(c_1 + a\bigl(c_2 + \cdots + a\,(c_{n-1})\bigr)\bigr),$$

ゼロから始めて、格納された係数ごとに 1 回の乗算と 1 回の加算を行います。Horner 法は、係数の前処理なしに一般の多項式を評価する方法の中で乗算回数が最少です。¹1 次数 4 以下での最適性は Ostrowski が（1954 年）、すべての次数での最適性は Pan が（1966 年）証明しました。次数 $n$ の一般の多項式を評価するアルゴリズムは、少なくとも $n$ 回の乗算と $n$ 回の加算を必要とします。 必要なのは AddMonoid + Mul だけなので任意の係数型で評価でき、累乗 $a^i$ を作らないため、固定幅整数では中間値が小さく保たれ、浮動小数点では条件が良くなります。

合成は多項式における評価

substitute(q) は同じ Horner ループを多項式の算術で実行し、$p(q) = \sum_i c_i q^i$ を計算します。$R$ が可換なら、これは評価準同型 $\mathrm{ev}_q : R[x] \to R[x]$, $x \mapsto q$ であり、$R$ を固定し $x$ を $q$ に送る唯一の環準同型です。単項式上では

$$\mathrm{ev}_q(x^i \cdot x^j) = q^{i+j} = q^i q^j = \mathrm{ev}_q(x^i)\,\mathrm{ev}_q(x^j),$$

であり、両辺は双線形なので $\mathrm{ev}_q(fg) = \mathrm{ev}_q(f)\,\mathrm{ev}_q(g)$ および $\mathrm{ev}_q(f + g) = \mathrm{ev}_q(f) + \mathrm{ev}_q(g)$ が成り立ちます。合成は結合的で $p(q(r)) = (p(q))(r)$ となります。両辺は $x$ 上で一致する準同型だからです。

$\deg p = n$、$\deg q = d$ のとき、Horner 法の第 $j$ ステップは次数 $jd$ の多項式に $q$ を掛けるので、筆算式のコストは $\sum_{j=1}^{n} (jd + 1)(d + 1) = O(n^2 d^2)$ であり、結果の次数は高々 $nd$ です。

ℕ からの標準写像による形式微分

derivative は $D\bigl(\sum_i c_i x^i\bigr) = \sum_{i \ge 1} (i \cdot 1_R)\, c_i\, x^{i-1}$ を計算し、整数 $i$ を NatHomomorphism::from_nat で $R$ に写します。$D$ は加法的でライプニッツ則を満たします。単項式上では

$$\begin{aligned} D(x^i x^j) &= (i + j)\, x^{i+j-1} \\ &= i\,x^{i-1} x^j + x^i\, j\,x^{j-1} \\ &= D(x^i)\,x^j + x^i\,D(x^j), \end{aligned}$$

であり、$D(fg) = D(f)\,g + f\,D(g)$ の両辺は $(f, g)$ について双線形なので、この規則はすべての多項式に拡張されます。$i$ は標数を法として取られるので、標数 $p$ では代数学のとおり $D(x^p) = p\,x^{p-1} = 0$ となります。luna-generic は Float、Double、BigInt に対する NatHomomorphism を提供しています。固定幅整数の係数は、そのようなインスタンスが与えられるまで微分できません（このトレイトは上流で FromNat に置き換えられつつあります）。

二分累乗法による累乗

pow(e) はパッケージ内部の pow_nat に委譲します。pow_nat は exp を半分にしながら不変条件 $\text{state} \cdot \text{factor}^{\text{exp}} = p^e$ を保ちます。$\lfloor \log_2 e \rfloor$ 回の二乗と、$e$ の立っているビットごとに高々 1 回の state への乗算を行うので、多項式の乗算は高々 $2\lfloor \log_2 e \rfloor + 1$ 回です。オペランドが大きくなっていくため最後の二乗が支配的であり、筆算式の乗法では総コストは $O\bigl((e\,\deg p)^2\bigr)$ です。

pow(0) はゼロを含むすべての $p$ に対して one() を返すので、$0^0 = 1$ です。これは @arithmetic.PowNatChecked の規約であり、DensePolynomial は Ok(self.pow(e)) を返すことでこれを実装しています。

コンテナのための構造的な順序

導出された Compare は、まず格納された長さで、次に定数項から係数ごとに順序付けます。切り詰めにより長さは次数に 1 を足したものなので、次数の低い多項式が先に並びます。この順序は多項式をソート済みマップやセットのキーにするためのものであり、環の順序ではなく、+ や * とも両立しません。

正しさ / 不変条件

• 正規形。 すべての公開演算の後、格納された最後の係数は非ゼロであり、ゼロは空です。したがって等価性は多項式としての等しさです。
• 環の法則。 係数型が可換環の法則を満たすとき、DensePolynomial[A] もそれを満たします。演算は正規形に対する教科書どおりの式だからです。プロパティテストで加法の単位元と正規化の冪等性を確認しています。
• 一致性。 karatsuba(a, b) == a * b であり、substitute は $\mathrm{ev}_q$ であり、derivative はライプニッツ則を満たし、pow(e) は $e$ 回の積に等しくなります。
• 値のセマンティクス。 レシーバや引数を変更するメソッドはなく、コンストラクタは入力配列をコピーします。
• 計算量（長さ $m, n$）: +、-、neg、scale は $O(m + n)$、* は $O(mn)$、karatsuba はしきい値を超える均衡したオペランドに対して $O(n^{\log_2 3})$、eval は $O(n)$、substitute は $O(n^2 d^2)$、pow(e) は $O(\log e)$ 回の乗算です。

採用しなかった代替案

• * の内部での Karatsuba。 すべての乗算に Neg + One を課し、短い積のコスト特性も変えてしまいます。明示的なメソッドにしておけば、選択は呼び出し側に委ねられます。
• FFT や数論変換による乗算。 これらは係数型に 1 の冪根や適切な法を必要としますが、ジェネリックな A はそれを提供しません。
• 次数を別に格納し、末尾のゼロを許すこと。 一部の演算で切り詰めを省けますが、等価性とハッシュが正規化に依存するようになります。演算ごとに 1 回切り詰めるのは安価で、この問題自体をなくせます。
• 疎な一変数型。 $x^{1000} + 1$ のような多項式ではここでは領域が無駄になります。そうした多項式には、要素数 1 の指数ベクトルを持つ SparsePolynomial を使ってください。

境界

• 一変数のみです。多変数多項式は TermPolynomial と SparsePolynomial です。
• 除算、GCD、因数分解、根の探索、補間はありません。
• 厳密な切り詰め: 係数が取り除かれるのは、それがゼロと == のときだけです。Float や Double の係数では、打ち消し合うべきなのに丸め誤差を含む先頭係数は残るため、格納された次数が数学的な次数を上回ることがあります。
• 固定幅整数の係数はラップアラウンドするので、多項式は $\mathbb{Z}/2^k$ 上のものとなり、上で示したように零因子と次数の低下が生じます。
• pow は UInt の指数を取ります。負の累乗はありません。

Footnotes

1. 次数 4 以下での最適性は Ostrowski が（1954 年）、すべての次数での最適性は Pan が（1966 年）証明しました。次数 $n$ の一般の多項式を評価するアルゴリズムは、少なくとも $n$ 回の乗算と $n$ 回の加算を必要とします。 ↩