sangi の多項式は内部でどのように表現されていますか？

Polynomial は係数配列 std::vector による dense 表現を採用し、coeffs_[i] が x^i の係数となる昇べき順で格納されます。dense 表現は係数アクセス・加減・Horner 評価・FFT 乗算をすべて O(1) のインデックス参照で実行できるためです。極めて疎な多項式 (例: x^1000 + 1) ではメモリ効率が下がりますが、sangi の対象ワークロード (低〜中次数の解析・求根・因数分解) では dense が定数倍で勝ちます。

Horner 法とは何ですか？

多項式 p(x) = c_0 + c_1 x + ... + c_n x^n を p(x) = c_0 + x(c_1 + x(c_2 + ... + x c_n)) のように入れ子の積和に書き換え、n 回の乗算と n 回の加算で評価する方法です。sangi の Polynomial::operator()(x) はこれを採用し、計算量は O(n)。直接的な x^k の累乗計算は使いません。

sangi は多項式乗算で FFT を使いますか？

Polynomial のデフォルト乗算は O(nm) の素朴乗算です。係数が浮動小数点で次数が大きい場合は FFT 畳み込みで O(n log n) に下げる経路がありますが、Polynomial の対象規模 (典型的に次数数十〜数百) では定数倍を含めた素朴乗算が実装上勝つため、デフォルトは Toom/FFT を呼び出しません。整数多倍長乗算 (Int) は別経路で Toom-Cook と NTT を使います。

擬剰余 (pseudo-remainder) と通常の剰余の違いは？

通常の剰余は係数体上での除算を仮定するので Z[x] では分数が現れます。擬剰余は割られる側に lc(g)^(deg f - deg g + 1) を予め乗じてから割ることで結果を整数係数のまま保ちます。sangi の pseudoRemainder() は intPolyGcd() の subresultant PRS 用の基本演算として用意されています。

Polynomial 算術アルゴリズム

概要

sangi の Polynomial<T> は係数型 $T$ を任意の数値型 (double, Rational, Complex<T>, Int など) に取れる多項式テンプレートである。基本的な加減乗除と評価・微分・積分はすべてヘッダオンリーで提供される。

API リファレンスは Polynomial API を参照。ここでは内部表現と算術アルゴリズムの選択を扱う。

内部表現 (密表現, 昇べき順)

Polynomial<T> は係数を std::vector<T> による dense 配列で保持する。添字は昇べき順で、coeffs_[i] が $x^i$ の係数 $c_i$ に対応する。

$$p(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_n x^n$$

密表現 (dense) と疎表現 (sparse)

多項式の表現には大きく 2 つの選択肢がある:

密表現 (dense): 全係数 (零も含めて) を昇べき順に配列で持つ。係数アクセス・加減・Horner 評価が $O(1)$ のインデックス操作。
疎表現 (sparse): 非零項のみを (指数, 係数) のペアで持つ。$x^{10^6} + 1$ のような疎多項式に有利。

sangi は密表現を採用する。対象ワークロード (求根・因数分解・微分積分・直交多項式) では次数 $n$ が数十〜数百のオーダーであり、全係数を連続メモリに置く密表現のほうがキャッシュ局所性・ベクトル化・FFT 経路への移行のいずれでも有利になる。極端に疎な多項式は RationalFunction や MultivariatePolynomial の疎形式で扱える。

昇べき順を選ぶ理由

$x^i$ の係数を添字 $i$ にそのまま対応させることで:

$p(x) \cdot x^k$ がベクトルの先頭挿入 (左シフト) ではなく末尾 push_back になり償却 $O(1)$ になる
FFT/NTT 多項式乗算と同じインデックス規約になり、変換層が不要
微分 $p'$ が「添字 $i$ の係数 $\cdot i$ を添字 $i-1$ に置く」という単純な操作で書ける

降べき順 (Horner の教科書記法) は出力時にのみ用い、内部では一切使わない。

加減算

2 つの多項式 $p = \sum c_i x^i$ と $q = \sum d_i x^i$ の和は、同じ次数の係数同士を足すだけである:

$$(p + q)(x) = \sum_{i=0}^{\max(\deg p, \deg q)} (c_i + d_i) x^i$$

計算量: $O(\max(n, m))$

結果の coeffs_ は両者の最大長に揃え、加算後に末尾の零係数を normalize() で取り除いて次数を確定する。引き算も同様で、$+$ を $-$ に置き換えるだけで実装は対称。

乗算 — 素朴・Karatsuba・FFT

多項式の乗算は係数列の畳み込みとして定義される:

$$(p \cdot q)(x) = \sum_{k=0}^{n+m} \left( \sum_{i+j=k} c_i d_j \right) x^k$$

素朴乗算 (デフォルト)

$O(nm)$ の二重ループによる直接畳み込み。$n, m$ が数十〜数百のオーダーでは定数倍が小さく実装上最速。 Polynomial::operator* はこの経路をデフォルトで使う。

計算量: $O(nm)$ の係数乗算と加算

Karatsuba 多項式乗算

多項式を上位・下位に分割して 4 回ではなく 3 回の半サイズ乗算に帰着する。係数長 (係数自体の bit 数) と次数 (項数) の両方が大きい場合に効く。 Karatsuba は項数に対する $O(n^{\log_2 3}) = O(n^{1.585})$ であり、係数 $T$ が Int や Rational のとき項数が大きくなると素朴乗算より速くなる。

sangi の Int 乗算 (係数長方向の Karatsuba) と、項数方向の多項式 Karatsuba は独立に効く。係数長が大きく次数も大きい場合、両者の積として全体計算量が抑えられる。

FFT 多項式乗算

多項式の係数列を「信号」と見なし、$\hat p \cdot \hat q$ の点ごと積を逆 FFT で戻す方法。浮動小数点係数の場合、長さ $N \geq n + m + 1$ の FFT で:

$$p \cdot q = \mathcal{F}^{-1}\!\left( \mathcal{F}(p) \odot \mathcal{F}(q) \right)$$

計算量: $O(n \log n)$

$T = $ double や Complex<double> で次数が十分大きい場合に有効だが、Polynomial の典型的な利用 (次数数百以下) では FFT の定数倍と前処理コストが素朴乗算に対して優位を取りにくい。整数係数を厳密に保ちたい場合は NTT (数論変換) を使うことで丸め誤差を回避できる。

整数の多倍長乗算で sangi が NTT を使っている理由・実装の詳細は Int 乗算 — NTT を参照。多項式と整数では「長さ」の意味が異なる (前者は項数、後者はワード数) ため、閾値や経路選択は別に管理される。

関連記事: 多項式の基本演算 / 多項式の表現

Horner 法による評価

多項式 $p(x) = c_0 + c_1 x + \cdots + c_n x^n$ の値を素朴に $x^k$ を別々に計算して足すと $O(n^2)$ あるいは累乗の事前計算で $O(n)$ の乗算と $O(n)$ の加算が必要になる。 Horner 法は入れ子の形に書き換えることで、ちょうど $n$ 回の乗算と $n$ 回の加算で評価する:

$$p(x) = c_0 + x \bigl( c_1 + x \bigl( c_2 + x ( \cdots + x \, c_n) \bigr) \bigr)$$

sangi の Polynomial::operator()(x) はこの形を採用する。ループは降べき順に展開する:

T result = coeffs_[n];
for (size_t i = n; i-- > 0; )
    result = result * x + coeffs_[i];
return result;

計算量: $n$ 回の乗算 + $n$ 回の加算 (= $O(n)$ の係数演算)

x が多項式や行列の場合

operator()(x) の x は係数型 $T$ と異なる型でもよい。 $x$ に Polynomial<T> を渡せば多項式合成 $p(q(x))$ になり、行列を渡せば行列多項式 $p(A)$ になる。 Horner 形は乗算の結合性のみを要求するので、これらすべてで同じコードが使える。

Estrin 法 (並列性が高い変種)

evaluateEstrin(x) は Horner と異なり、$x^2, x^4, x^8, \ldots$ をまず計算し、ペアごとに $c_{2i} + c_{2i+1} x$ をまとめてから累乗で組み合わせる。乗算回数は Horner より少し増えるが、依存チェーンが浅く、SIMD やパイプライン並列で高速化できる。 Horner との切り替えは利用者が明示的に選ぶ。

関連記事: 多項式の基本演算 — Horner 法

合成除算 (synthetic division)

除数が 1 次式 $x - r$ の場合、多項式除算は係数列に対する単純な前向き漸化式に縮約される。これを合成除算 (Ruffini-Horner) と呼ぶ:

$$p(x) = (x - r) \cdot q(x) + p(r)$$

$q(x) = \sum_{i=0}^{n-1} b_i x^i$ の係数は次の漸化式で得られる:

$$b_{n-1} = c_n, \quad b_{i-1} = c_i + r \cdot b_i \quad (i = n-1, \ldots, 1)$$

そして剰余は $p(r) = c_0 + r \cdot b_0$。この漸化式はHorner 評価そのものであり、評価と除算が同じ手続きで実現される。

計算量: $O(n)$

有理根定理で求めた根の候補 $r$ から $(x - r)$ を取り除くとき、sangi の因数分解パイプラインはこの合成除算を使う。汎用の divmod() は $O(n \cdot m)$ の長除法であり、除数が 1 次の場合は内部的に合成除算経路に分岐する。

polyder / polyint と自動微分

多項式の微分・積分は係数の単純な変換で計算できる。

微分 (polyder)

$$p(x) = \sum_{i=0}^n c_i x^i \;\Rightarrow\; p'(x) = \sum_{i=1}^n i \cdot c_i \, x^{i-1}$$

sangi の derivative() は係数配列を 1 つ縮め、各要素に添字を掛ける $O(n)$ 操作。 $n$ 階微分 derivative(n) は係数 $c_i$ に $i!/(i-n)!$ を掛けて先頭 $n$ 要素を捨てるだけで一括で計算できる。

積分 (polyint)

$$\int p(x)\,dx = C + \sum_{i=0}^n \frac{c_i}{i+1} x^{i+1}$$

sangi の integral() は積分定数 $C = 0$ で固定し、係数配列の先頭に $0$ を挿入してから各要素を $1/(i+1)$ 倍する。 $T = $ Rational なら厳密、$T = $ double なら浮動小数。

自動微分との関係

$p(x_0)$ と $p'(x_0)$ を同時に求めたい場面では、Horner 法を二重に走らせるテクニックが使える:

T fx = coeffs_[n], dfx = T{0};
for (size_t i = n; i-- > 0; ) {
    dfx = dfx * x + fx;       // p' を進める
    fx  = fx  * x + coeffs_[i]; // p を進める
}
// fx = p(x), dfx = p'(x)

これは 1 次の前進モード自動微分 (forward-mode AD) を Horner で展開したものと同等で、 Newton 法・Halley 法の内部評価で使われる ($p(x)/p'(x)$ を 1 パスで計算)。

API としては derivative() で多項式を取り出してから operator() を 2 回呼ぶこともできるが、反復求根の内側ループでは Horner 二重展開のほうが配列確保を避けられて高速。

擬除算と擬剰余

整数係数多項式の GCD や因数分解では、商や剰余に分数が現れることを避けたい。通常の多項式除算は係数体 ($\mathbb{Q}$ や $\mathbb{R}$) を仮定するため、$\mathbb{Z}[x]$ の中だけでは閉じない。

擬除算 (pseudo-division) は除数の主係数の冪を先に乗じることで剰余を整数係数に保つ操作である:

$$\mathrm{lc}(g)^{\deg f - \deg g + 1} \cdot f \;=\; q \cdot g + r, \quad \deg r < \deg g$$

例: $f = x^2,\ g = 2x + 1$ を $\mathbb{Z}[x]$ で擬除算する。 $\mathrm{lc}(g)^{2-1+1} = 4$ を掛けてから割ると:

$$4 x^2 = (2x - 1)(2x + 1) + 1$$

擬商 $q = 2x - 1$、擬剰余 $r = 1$ がいずれも整数。

サブレザルタント PRS

素朴な擬剰余列 (PRS) は反復のたびに主係数の冪が積み重なり、係数が指数的に爆発する。 サブレザルタント PRS (subresultant pseudo-remainder sequence) は、毎ステップで適切な「共通因子」を除くことで係数の伸びを多項式に抑える。 sangi の intPolyGcd() はこのアルゴリズムで実装されており、 pseudoRemainder() はその基本演算として公開されている。

関連記事: 多項式の GCD 計算 / 終結式

設計判断

密表現 + 昇べき順を統一する理由

多項式は計算機代数 (CAS)・数値解析・信号処理のいずれでも基本オブジェクトであり、表現規約が分かれると相互運用のコストが大きい。 sangi は全モジュールで密表現 + 昇べき順に統一しており、 Polynomial・直交多項式・Chebyshev 補間・FFT・求根のいずれも同じ coeffs_ 規約を共有する。

乗算経路を素朴に固定する理由

多項式乗算で Karatsuba/FFT が定数倍を含めて素朴乗算に勝つのは、 $T$ が浮動小数で次数が数百以上、または $T$ が Int/Rational で次数が数十以上の領域。対象ワークロード (求根・因数分解の内部演算、Chebyshev/Legendre の漸化生成) はその閾値以下に収まることが多いため、 Polynomial の operator* はデフォルトで素朴 $O(nm)$ を使う。必要なら FFT 畳み込みは math/fft 経由で呼べる。

参考文献

Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997.
von zur Gathen, J. and Gerhard, J. Modern Computer Algebra. 3rd ed. Cambridge University Press, 2013.
Brown, W. S. (1971). "On Euclid's algorithm and the computation of polynomial greatest common divisors". Journal of the ACM, 18(4), 478–504.