Polynomial<T> — 多項式テンプレート

概要

Polynomial<T> は calx の多項式テンプレートクラスである。 係数型 T に double, Rational, Complex<double> など任意の数値型を使用できる。

内部表現は昇べきの順: coeffs_[i] が $x^i$ の係数。

• 任意の係数型 — double から Rational (厳密計算) や Complex<T> まで
• Horner 評価 — $O(n)$ で $p(x)$ を計算
• 多項式除算 — 商と剰余 (divmod)
• ユークリッド GCD — 多項式の最大公約因子
• 微分・積分 — derivative(), integral()
• 直交多項式 — Chebyshev, Legendre, Hermite, Laguerre, Bessel, Jacobi, Gegenbauer
• C++23 concepts — PolynomialCoeff で型制約

ビルド

必要なヘッダ

#include <math/core/Polynomial.hpp>

ヘッダオンリー。リンクするライブラリは不要。

コンストラクタ

Polynomial<double> p = {1.0, 2.0, 3.0};  // 1 + 2x + 3x^2
Polynomial<Rational> q = {Rational(1,2), Rational(3,4)};  // 1/2 + 3/4 x

アクセサ

算術演算子

スカラーとの演算: p + c, c * p, p / c もサポート。 複合代入: +=, -=, *=, /=。

除算と剰余

auto [q, r] = p.divmod(d);  // p = d * q + r, deg(r) < deg(d)

除算モード — operator/ の動作切り替え

Polynomial / Polynomial の戻り値型は、コンパイル時マクロで切り替えられる。

切り替え方法 1: ソースコード内で定義

// RationalFunction.hpp をインクルードする前に定義
#define CALX_POLY_DIV_RATIONAL
#include <math/core/RationalFunction.hpp>

Polynomial<double> p = {1, 0, 1};  // x^2 + 1
Polynomial<double> q = {1, 1};     // x + 1

auto rf = p / q;       // RationalFunction<double>: (x^2+1)/(x+1)
auto quot = p.quo(q);  // Polynomial<double>: x - 1 (商が必要な場合)

切り替え方法 2: CMake で設定

# CMakeLists.txt
option(CALX_POLY_DIV_RATIONAL
       "Polynomial / Polynomial returns RationalFunction" OFF)

if(CALX_POLY_DIV_RATIONAL)
    add_compile_definitions(CALX_POLY_DIV_RATIONAL)
endif()

この設定でプロジェクト全体の operator/ の動作が切り替わる。 どちらのモードでも quo() と divmod() は常に多項式除算として使える。

評価

p(x) の x は T と異なる型でもよい (例: Polynomial<double> に Complex<double> を代入)。

x に多項式を渡すと 多項式合成 $p(q(x))$ になる。 compose(p, q) と等価。

Polynomial<double> p({1, 0, 1});  // 1 + x^2
Polynomial<double> q({0, 1, 1});  // x + x^2
auto r = p(q);  // p(q(x)) = 1 + (x + x^2)^2 = 1 + x^2 + 2x^3 + x^4

合成後の次数は $\deg(p(q)) = \deg(p) \cdot \deg(q)$ となる。

微分・積分

代数的操作

因数分解・終結式

整数係数多項式の厳密な因数分解と、終結式・判別式を提供する。 ヘッダ math/core/Polynomial_factorization.hpp をインクルードする。

#include <math/core/Polynomial.hpp>
#include <math/core/Polynomial_factorization.hpp>
#include <math/core/FactoredInteger.hpp>

factorize() — 既約因子への分解

整数多項式 $f(x) \in \mathbb{Z}[x]$ を既約因子に分解する。

係数型 T は int または int64_t。

パイプライン:

1. content 抽出 — 全係数の GCD を定数因子として分離
2. 有理根定理 — $f(b/a) = 0$ を $b \mid a_0$, $a \mid a_n$ の範囲で探索し、$(ax - b)$ の 1 次因子を全て取り除く
3. Yun の無平方分解 — $\gcd(f, f')$ を用いて $f = \prod f_i^i$ ($f_i$ は無平方かつ互いに素)
4. Kronecker 法 — 各無平方因子について $\deg f / 2 + 1$ 点での値の約数の組合せから Lagrange 補間で因子候補を構成

// x^4 - 1 = (x - 1)(x + 1)(x^2 + 1)
Polynomial<int> f = {-1, 0, 0, 0, 1};
auto factors = factorize(f);

for (size_t i = 0; i < factors.size(); ++i) {
    std::cout << factors[i] << "^" << factors.exponent(i) << std::endl;
}
// (-1 + x)^1
// (1 + x)^1
// (1 + x^2)^1

上記は core 版 (ヘッダのみ、リンク不要) で、 Kronecker 法の組合せ上限 100,000 により $\deg f \lesssim 10$ が実用範囲。

高速版 (calx_factorize.lib)

calx_factorize をリンクすると、自動的に高速パイプラインに切り替わる:

1. $\mathrm{GF}(p)$ 因数分解 — Cantor-Zassenhaus (DDF + EDF)
2. Hensel lifting — 線形 + 多因子二分木 Hensel で $\mathbb{Z}[x]$ に引き上げ
3. 因子復元 — $r \leq 12$: Zassenhaus 組合せ、$r > 12$: van Hoeij 格子復元 (LLL)

Int 型は内部で int64_t に委譲 (全係数が 64-bit に収まる場合)。 Rational 型は分母の LCM で整数化してから Int パスに委譲する。

// CMake: target_link_libraries(your_target PRIVATE calx_factorize)

squareFreeFactorization() — Yun の無平方分解

Yun のアルゴリズム (1976) を $\gcd(f, f')$ の subresultant PRS で実装。 重複因子を冪指数付きで分離する。既約性は保証しない。

整数多項式 GCD と擬剰余

擬剰余と通常の剰余の違い

通常の剰余は係数体 (体上での割り算) を仮定するので、整数係数では使えない:

例: $f = x^2,\ g = 2x + 1$ を整数で割ろうとすると $q = \dfrac{x}{2} - \dfrac{1}{4}$ となり、商に分数が現れる ($\mathbb{Z}[x]$ では商も剰余も求まらない)。

擬剰余 (pseudo-remainder) は、先に $g$ の主係数の冪を掛けることで分数を回避する:

主係数 $\mathrm{lc}(g)$ を予め乗算しておくため、$g$ で割っても係数が整数のまま保たれる。

同じ例で: $\mathrm{lc}(g)^{2-1+1} = 2^2 = 4$ を掛けてから割ると

擬剰余 $r = 1$ (整数のみ)。

calx の pseudoRemainder() は intPolyGcd() (整数多項式 GCD) のために用意されている — 整数係数のまま GCD を計算する subresultant PRS の基本演算である。

終結式・判別式

// f(x) = x^2 - 2 の判別式 = 4·1·2 = 8
Polynomial<int> f = {-2, 0, 1};
std::cout << discriminant(f) << std::endl;  // 8

判別式が $0$ でなければ $f$ は無平方 (重根を持たない)。

直交多項式

漸化式により生成される古典的な直交多項式。

auto T5 = chebyshevT<double>(5);
std::cout << T5 << std::endl;
// 5x - 20x^3 + 16x^5

型特性

使用例

#include <math/core/Polynomial.hpp>
#include <iostream>
using namespace calx;

int main() {
    // 多項式の構築
    Polynomial<double> p = {1.0, -3.0, 2.0};  // 1 - 3x + 2x^2
    Polynomial<double> q = {1.0, 1.0};         // 1 + x

    // 四則演算
    std::cout << "p = " << p << std::endl;
    std::cout << "q = " << q << std::endl;
    std::cout << "p * q = " << p * q << std::endl;

    // 除算
    auto [quot, rem] = p.divmod(q);
    std::cout << "p / q = " << quot << ", remainder = " << rem << std::endl;

    // 評価 (Horner)
    std::cout << "p(2.0) = " << p(2.0) << std::endl;  // 3.0

    // 微分・積分
    std::cout << "p' = " << p.derivative() << std::endl;  // -3 + 4x
    std::cout << "∫p dx = " << p.integral() << std::endl;

    // GCD
    Polynomial<double> a = {-1.0, 0.0, 1.0};  // x^2 - 1 = (x-1)(x+1)
    Polynomial<double> b = {-1.0, 1.0};        // x - 1
    std::cout << "gcd(x^2-1, x-1) = " << gcd(a, b) << std::endl;

    // 直交多項式
    auto L3 = legendre<double>(3);
    std::cout << "P_3(x) = " << L3 << std::endl;
}

関連する数学的背景

以下の記事では、Polynomial クラスの基盤となる数学的概念を解説している。

• 多項式の表現と基本演算 — 密/疎表現・加減乗除
• 多項式の GCD 計算 — ユークリッド・部分終結式アルゴリズム
• 因数分解（有限体上） — Berlekamp・Cantor-Zassenhaus 法
• 因数分解（整数上） — Hensel 持ち上げ・LLL 格子基底縮小