Root Finding — 求根アルゴリズム

概要

sangi の求根モジュールは 3 つのカテゴリに分かれる。

• 1D 求根 — $f(x) = 0$ を満たす実数 $x$ を求める。二分法から Brent 法まで 10+ 種類
• 多項式求根 — $p(x) = 0$ の全根 (複素根含む) を求める。低次は公式解、高次は反復法
• N-D 求根 — $\mathbf{F}(\mathbf{x}) = \mathbf{0}$ を満たすベクトル $\mathbf{x}$ を求める

ビルド

// 全モジュールまとめてインクルード
#include <math/roots/root_finding.hpp>

// 個別インクルードも可能
#include <math/roots/root_finding_1d.hpp>   // 1D 求根
#include <math/roots/root_finding_nd.hpp>   // N-D 求根
#include <math/roots/polynomial_roots.hpp>  // 多項式求根

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

結果型と収束判定

RootFindingResult<R>

1D 求根関数、 多項式根反復関数、 N-D 求根関数の 統合戻り値型。 math/roots/root_finding_base.hpp 定義。 テンプレートパラメータ R がスカラー型 T なら 1D 用、 ベクトル型 V なら N-D 用に自動分岐する。

template<typename R>
struct RootFindingResult {
    using error_type = /* R がスカラーなら R、 ベクトルなら R::value_type */;

    std::optional<R>  root;           // 成功時: 根の近似値 (スカラーまたはベクトル)、 失敗時: std::nullopt
    bool              converged;      // 収束したかどうか
    size_t            iterations;     // 実際に行った反復回数
    error_type        error_estimate; // 推定誤差 (絶対値、 収束時は最終ステップ幅 等)
};

1D 求根は RootFindingResult<double> を、 N-D 求根は RootFindingResult<Vector<double>> を返す。

使い方:

auto r = brent_method(f, a, b);
if (r.converged && r.root) {
    double x = *r.root;            // または r.root.value()
    // ... x を使う
} else {
    // 失敗 (区間に符号変化なし / max iter 到達 / 数値特異 等)
    std::cerr << "no root: iter=" << r.iterations
              << " err=" << r.error_estimate << '\n';
}

ConvergenceCriteria<T>

収束判定の閾値を保持する struct。 全求根関数の最後の引数として渡せる (省略時はデフォルト)。

template<concepts::OrderedField T>
struct ConvergenceCriteria {
    T      abs_ftol = std::numeric_limits<T>::epsilon() * 100;   // |f(x)| 絶対許容誤差
    T      abs_xtol = std::numeric_limits<T>::epsilon() * 100;   // |Δx|   絶対許容誤差
    T      rel_ftol = std::numeric_limits<T>::epsilon() * 1000;  // |f|相対許容誤差
    T      rel_xtol = std::numeric_limits<T>::epsilon() * 1000;  // Δx 相対許容誤差
    size_t max_iterations = 100;                                  // 反復上限
};

使用例:

// 緩い許容誤差で素早く打ち切る
ConvergenceCriteria<double> loose;
loose.abs_ftol = 1e-6;
loose.max_iterations = 50;
auto r = newton_raphson(f, df, x0, loose);

// 多倍長 (Float) で高精度
ConvergenceCriteria<Float> tight;
tight.abs_ftol = Float("1e-50");
tight.max_iterations = 200;
auto r2 = brent_method<Float>(f, a, b, tight);

1D: 囲い込み法 (Bracketing)

区間 $[a, b]$ 内で $f(a)$ と $f(b)$ が異符号であることを前提とし、確実に収束する。

関数詳細

bisection

template<concepts::OrderedField T>
RootFindingResult<T> bisection(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 区間 $[a, b]$ を毎ステップ半分に絞り、$f$ の符号が変わる側を残す。 線形収束 (毎反復で誤差が 1/2)。

パラメータ:

返り値: RootFindingResult<T>。 区間に符号変化がない場合は converged = false、 root = std::nullopt。

用途: ロバスト性が最優先で速度は二の次のとき。 区間幅が確実に半減するため最悪反復回数を予測可能 (log2((b-a)/tol))。

auto r = bisection<double>([](double x){ return x*x - 2.0; }, 1.0, 2.0);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=42  root=1.41421356237311
//   ※ 真値 sqrt(2) = 1.41421356237310 と 12 桁一致

regula_falsi

template<concepts::OrderedField T>
RootFindingResult<T> regula_falsi(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 偽位置法 (Regula Falsi) の基本形。 区間端点を結ぶ直線と $x$ 軸の交点 $c$ を次の試行値とし、 $f(c)$ と同符号の端点を $c$ で置換して符号変化の括り区間を保つ。 常に根を挟むため収束は保証されるが、 $f$ が凸 (または凹) の場合は一方の端点が固定され続ける停滞 (stagnation) が起こり、 収束が線形に劣化しうる。 停滞を回避した実用版は illinois_method を参照。

パラメータ:

用途: 教育・参照用の素の Regula Falsi。 実用では停滞を回避する illinois_method (Illinois)、 king_method (Pegasus 系)、 anderson_bjork_method を推奨。

illinois_method

template<concepts::OrderedField T>
RootFindingResult<T> illinois_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 偽位置法 (基本形は regula_falsi) に Illinois 修正を加えた版。 純粋な Regula Falsi は凸/凹な関数で一方の端点が停滞し線形収束に劣化するが、 同側の端点が 2 回連続で更新されなかった時に重み $1/2$ で fa または fb をスケールすることで停滞を回避し、 超線形収束を確保する。 さらに高速な改良版として anderson_bjork_method および king_method (Pegasus 法) を参照。

パラメータ:

ridders_method

template<concepts::OrderedField T>
RootFindingResult<T> ridders_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 区間の中点 $m$ と $e^{Q(x-m)}$ で重み付けした false position を組合せる。 1 反復で $f$ を 2 回評価する代わりに $\sqrt{2} \approx 1.414$ 次の超線形収束を得る (実効収束次数は約 1.84)。

パラメータ:

用途: bisection より速く、 Brent より実装が単純。 関数評価が安価な場合の中庸選択。

brent_method

template<concepts::OrderedField T>
RootFindingResult<T> brent_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 逆 2 次補間 / 割線法 / 二分法を内部状態に応じて切替える Brent (1973) のハイブリッド。 区間幅が縮まないステップは bisection で強制縮退するため、 速度と保証収束を両立する。

パラメータ:

用途: 1D 求根のデフォルト推奨。 SciPy brentq 等の標準アルゴリズム。 区間端点で符号が異なれば常に確実に根を返す。

auto r = brent_method<double>([](double x){ return std::cos(x) - x; }, 0.0, 1.0);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=26  root=0.739085133215161
//   ※ Dottie 数 = 0.7390851332151607 と 14 桁一致

alefeld_potra_shi

template<concepts::OrderedField T>
RootFindingResult<T> alefeld_potra_shi(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Alefeld-Potra-Shi (1995) の TOMS Algorithm 748。 逆 3 次補間と double-length step を組合せて Brent より高い実効収束次数 ($\approx 1.84$ for smooth $f$) を持つ。

パラメータ:

用途: Brent と同じく区間保証付きで、 滑らかな $f$ では Brent より少ない反復で収束。 SciPy scipy.optimize.brentq の代替に推奨。

anderson_bjork_method

template<concepts::OrderedField T>
RootFindingResult<T> anderson_bjork_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Regula Falsi の改良版。 同側端点が 2 回連続で更新されなかったとき、 重み $m = 1 - f_c / f_b$ (もし $m \le 0$ なら $m = 1/2$ にフォールバック) を停滞側の関数値にかけて再補間する。 Illinois 法 (固定重み $1/2$) よりも収束が速く、 King-Pegasus 系と並ぶ実用的な超線形収束を達成する。

パラメータ:

用途: regula_falsi より高速な区間保証付き超線形法が欲しい場合の標準選択。 出典: Anderson & Björck, "A new high order method of regula falsi type for computing a root of an equation", BIT (1973)。

1D: 非区間法 (Non-bracketing)

初期値 $x_0$ (と任意で導関数) から反復する。収束は速いが初期値依存。

関数詳細

newton_raphson

template<concepts::OrderedField T>
RootFindingResult<T> newton_raphson(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 各反復 $x_{n+1} = x_n - f(x_n) / f'(x_n)$。 単根近傍で 2 次収束 ($|x_{n+1} - x^*| \le C |x_n - x^*|^2$)。 重根では線形に劣化。

パラメータ:

注意: $|f'(x)| < \epsilon \cdot 10$ の near-stationary 点で発散判定 → converged=false。 初期値が悪いと別の根に収束、 または発散することがある。 区間保証がほしい場合は brent_method を推奨。

// 例: f(x) = x^2 - 2 の根を Newton 法で求める
auto r = newton_raphson<double>(
    [](double x) { return x*x - 2.0; },
    [](double x) { return 2.0 * x; },
    1.0);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=5  root=1.41421356237310
//   ※ 真値 sqrt(2) = 1.41421356237310 と 15 桁一致 (5 反復で収束、 2 次収束の威力)

secant_method

template<concepts::OrderedField T>
RootFindingResult<T> secant_method(
    const std::function<T(T)>& f,
    T x0, T x1,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 2 点 $(x_{n-1}, x_n)$ を結ぶ割線 (secant) と $x$ 軸の交点を次の試行値: $x_{n+1} = x_n - f(x_n) \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$。 導関数を 1 次差分で代用するため Newton よりわずかに遅い ($\approx 1.618$ 次収束、 黄金比)。

パラメータ:

注意: $|f(x_n) - f(x_{n-1})| < \epsilon$ (傾きほぼゼロ) で停滞検出 → converged=false。

// 例: f(x) = x - cos(x) の根 (= Dottie 数)
auto r = secant_method<double>(
    [](double x) { return x - std::cos(x); }, 0.0, 1.0);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=6  root=0.739085133215161
//   ※ Dottie 数 0.7390851332151607 と 14 桁一致

fixed_point_iteration

template<concepts::OrderedField T>
RootFindingResult<T> fixed_point_iteration(
    const std::function<T(T)>& g,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 反復写像 $g$ に対して $x_{n+1} = g(x_n)$ を繰り返し、 不動点 $x^* = g(x^*)$ (すなわち $x = g(x)$ の解) を求める。 引数 $g$ は反復関数そのものであり、 $f(x)=0$ の $f$ ではない点に注意。 $g$ が $x^*$ 近傍で縮小写像 ($|g'(x^*)| < 1$) なら線形収束し、 特に $g'(x^*) = 0$ なら 2 次以上で収束する。 $|g'(x^*)| > 1$ では発散する。

パラメータ:

注意: 縮小写像でない ($|g'(x^*)| \ge 1$) 場合は発散し、 最大反復回数に達して converged=false となる (最良近似値は返る)。 newton_raphson は $g(x) = x - f(x)/f'(x)$ とおいた不動点反復の特殊形である。

// 例: g(x) = cos(x) の不動点 (= Dottie 数) を不動点反復で求める
auto r = fixed_point_iteration<double>(
    [](double x) { return std::cos(x); }, 0.0);
if (r.converged) std::cout << *r.root;
// g(x) = cos x は |g'(x)| = |sin x| < 1 ゆえ収束、 x* = 0.7390851332151607 (Dottie 数)

steffensen_method

template<concepts::OrderedField T>
RootFindingResult<T> steffensen_method(
    const std::function<T(T)>& f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 更新式 $x_{n+1} = x_n - \frac{f(x_n)^2}{f(x_n + f(x_n)) - f(x_n)}$。 分母 $f(x+f)-f(x)$ が $f(x)\,f'(x)$ を近似するため、 導関数なしで Newton と同じ 2 次収束を得る (Aitken の $\Delta^2$ 加速を不動点反復に適用)。 割線法と違い保持する点は 1 つだけ。 1 反復に $f$ を 2 回評価する。

パラメータ:

注意: 分母がほぼゼロで停滞検出 → converged=false。 初期値が根から遠いと発散する (大域収束保証はない)。

// 例: f(x) = x^2 - 2 の根 (= √2)、 導関数不要
auto r = steffensen_method<double>(
    [](double x) { return x * x - 2.0; }, 1.5);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=4  root=1.414213562373095

sidi_method

template<concepts::OrderedField T>
RootFindingResult<T> sidi_method(
    const std::function<T(T)>& f,
    T x0, T x1,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>(),
    size_t memory = 3);

動作: 直近 $K+1$ 点を通る Newton 補間多項式 $p$ の最新点での微分係数を使って $x_{n+1} = x_n - f(x_n)/p'(x_n)$ と更新する。 memory $=2$ で割線法 ($\varphi \approx 1.618$) に一致し、 増やすほど収束次数が上がる ($\to 2$)。 導関数不要。

パラメータ:

// 例: f(x) = x^2 - 2 の根 (= √2)、 2 次補間 (memory=3)
auto r = sidi_method<double>(
    [](double x) { return x * x - 2.0; }, 1.0, 2.0,
    ConvergenceCriteria<double>(), 3);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=6  root=1.414213562373095

inverse_quadratic_interpolation

template<concepts::OrderedField T>
RootFindingResult<T> inverse_quadratic_interpolation(
    const std::function<T(T)>& f,
    T x0, T x1, T x2,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 3 点を通る逆関数 $x = g(f)$ の二次補間多項式を $f = 0$ で評価して次点を得る ($x = \sum_i x_i \prod_{j \neq i} f_j/(f_i - f_j)$)。 Brent 法の内部要素を独立に使えるようにしたもの。 導関数不要・超 1 次収束 ($\approx 1.84$ 次) だが、 関数値が接近すると不安定なので本来は安全網 (Brent 法) と併用する。

パラメータ:

注意: 3 点の $f$ 値が接近すると補間が破綻 → converged=false。

// 例: f(x) = x^2 - 2 の根 (= √2)
auto r = inverse_quadratic_interpolation<double>(
    [](double x) { return x * x - 2.0; }, 1.0, 1.3, 1.7);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=5  root=1.414213562373095

muller_method

template<concepts::OrderedField T>
RootFindingResult<Complex<T>> muller_method(
    const std::function<Complex<T>(Complex<T>)>& f,
    Complex<T> z0, Complex<T> z1, Complex<T> z2,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 3 点を通る放物線 (2 次多項式) の根のうち最新点に近い方へ進む。 判別式が負でも複素数で進むため、 実数寄りの初期点からでも複素根に到達できる (多項式の複素根探索に有用)。 収束次数 $\approx 1.84$、 導関数不要。 $f$ は複素引数で評価できる必要がある。

パラメータ:

注意: 戻り値型は RootFindingResult<Complex<T>> (根が複素数)。 どちらの複素根に収束するかは初期点に依存する。

// 例: z^2 + 1 = 0 → ±i （実数寄りの初期点から複素根へ）
std::function<Complex<double>(Complex<double>)> g =
    [](Complex<double> z) { return z * z + Complex<double>(1.0); };
auto r = muller_method<double>(g, {1.0}, {0.5}, {0.0, 0.3});
if (r.converged) std::cout << r.root->re << " + " << r.root->im << "i";
// 実行結果: iter=2  root=0 + 1i  (= i)

1D: 高次収束法

関数詳細

halley_method

template<concepts::OrderedField T>
RootFindingResult<T> halley_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 更新式 $x_{n+1} = x_n - \frac{2 f f'}{2 f'^2 - f f''}$。 Newton の 2 次収束を 3 次に引き上げる Householder 1 次。 2 階導関数が必要、 計算コストは Newton の約 2 倍だが反復回数は減る。

パラメータ:

注意: 分母 $2 f'^2 - f f''$ がほぼゼロで発散判定。

// 例: f(x) = x^3 - 2 の根 (= 2^(1/3))
auto r = halley_method<double>(
    [](double x) { return x*x*x - 2.0; },
    [](double x) { return 3.0 * x * x; },
    [](double x) { return 6.0 * x; },
    1.0);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=3  root=1.25992104989487
//   ※ 真値 2^(1/3) = 1.2599210498948732 と 14 桁一致 (3 反復、 3 次収束の威力)

schroder_method

template<concepts::OrderedField T>
RootFindingResult<T> schroder_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 更新式 $x_{n+1} = x_n - \frac{f f'}{f'^2 - f f''}$。 重根に対して通常の Newton が 1 次に劣化するのに対し、 Schröder (1870) は 重根でも 2 次収束を維持する設計。

パラメータ:

用途: 多項式の重根を扱う場合、 Halley より安定。 単根では Halley (3 次) の方が速い。

king_method

template<concepts::OrderedField T>
RootFindingResult<T> king_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Regula Falsi で補間点を求めた後、 同じ端点が連続停滞した場合に固定端点の関数値を $f_1 \leftarrow f_1 \cdot f_2 / (f_2 + f)$ で縮小して再補間する。 系譜は Illinois 法 (固定重み $1/2$) → Pegasus 法 (Dowell & Jarratt, 1972; 重み $f_2/(f_2+f)$) → King 法 (1973)。 King 法は Pegasus の遅い副ステップを 1〜2 個取り除いて漸近的な収束効率を高めた改良版であり (1 次差分商のみで収束次数 $\approx 1.839$)、 king_method という命名は誤りではなく Pegasus を改良した King の方法を指す。 Illinois 法より速い超線形収束。

出典: R. F. King, "An improved Pegasus method for root finding", BIT 13, 423–427 (1973); M. Dowell & P. Jarratt, "The Pegasus method for computing the root of an equation", BIT 12, 503–508 (1972)。

パラメータ:

popovski_method

template<concepts::OrderedField T>
RootFindingResult<T> popovski_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 2 ステップ更新 (Newton ステップ + 補正項) で 1 階導関数のみから 4 次収束を実現。 King-Werner / Ostrowski 系の改良。 2 階導関数を計算しない代わり 1 反復あたり $f$ を 2 回評価する。

パラメータ:

用途: 2 階導関数が高コスト (e.g. AD で計算) なケースで、 反復回数より関数評価回数の方が問題な時。

householder_method

template<concepts::OrderedField T>
RootFindingResult<T> householder_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    const std::function<T(T)>& d3f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Householder 法系列の階数 $d = 3$ のメンバ。 Newton 法 ($d{=}1$)・Halley 法 ($d{=}2$) の延長で 4 次収束。 更新式 $x_{n+1} = x_n - \frac{6 f f'^2 - 3 f^2 f''}{6 f'^3 - 6 f f' f'' + f^2 f'''}$。 1〜3 階導関数を要する。

パラメータ:

注意: 分母 $6 f'^3 - 6 f f' f'' + f^2 f'''$ がほぼゼロで発散判定。 3 階導関数が安く得られる (記号微分・AD) 場合に Newton/Halley より少ない反復で収束する。

// 例: f(x) = x^3 - 2 の根 (= 2^(1/3))
auto r = householder_method<double>(
    [](double x) { return x*x*x - 2.0; },   // f
    [](double x) { return 3.0*x*x; },        // f'
    [](double x) { return 6.0*x; },          // f''
    [](double)   { return 6.0; },            // f'''
    1.5);
if (r.converged) std::cout << *r.root;
// 実行結果: iter=3  root=1.259921049894873  (= 2^(1/3))

多項式求根

全関数は std::vector<Complex<T>> を返す (実根は虚部 = 0)。 根はベクトル空間の元ではなく単なる複素数の集合なので、 sangi::Vector (線形代数用) ではなく STL の std::vector を使う。 入力 Polynomial<T> は係数を 昇べき順で保持 (p[0] = 定数項)。

低次 (公式解)

高次 (反復法)

関数詳細

solveLinear / solveQuadratic / solveCubic / solveQuartic

template<typename T> std::vector<Complex<T>> solveLinear(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveQuadratic(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveCubic(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveQuartic(const Polynomial<T>& p);

動作: それぞれ 1/2/3/4 次方程式に対する閉形式解。 次数不一致は assert で検出。

パラメータ:

桁落ち対策: solveQuadratic は判別式が小さい場合に大きい絶対値の根を直接計算し、 もう一方を Vieta の公式で導出 (Numerical Recipes 推奨)。 出典: 奥村「アルゴリズム事典」p205。

// 例: x^2 - 5x + 6 = (x-2)(x-3) = 0
Polynomial<double> p({6.0, -5.0, 1.0});  // 昇べき: c0=6, c1=-5, c2=1
auto roots = solveQuadratic(p);
for (auto& r : roots) std::cout << r.re << ' ';
// 実行結果: "3 2" (真値 {2, 3} と完全一致)

jenkinsTraub

template<typename T>
std::vector<Complex<T>> jenkinsTraub(const Polynomial<T>& poly);

動作: Jenkins-Traub (1970) の RPOLY アルゴリズム (3 段階 shift-shift-shift)。 各根を 1 つずつ求め、 デフレーション (因数除去) で次数を下げる。 同分野で数十年にわたり標準的な「最も信頼できる」多項式求根法。

パラメータ:

用途: 中次数 (5-19 次) のデフォルト推奨。 solvePolynomial も 5-19 次でこれを内部呼出。

// 例: x^6 - 1 = 0 (1 の 6 乗根)
Polynomial<double> p({-1.0, 0, 0, 0, 0, 0, 1.0});
auto roots = jenkinsTraub(p);
// 実行結果: 6 根 = {±1, ±exp(±iπ/3)}
//   |p(root)| 残差最大 4.97e-16 (16 桁精度)

laguerre

template<typename T>
std::vector<Complex<T>> laguerre(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: Laguerre 法 + デフレーション。 各反復で $L(z) = -n / \left( G \pm \sqrt{(n-1)(n H - G^2)} \right)$ 更新 ($G = p'/p$, $H = G^2 - p''/p$)。 任意の初期値から全根に大域収束する珍しい性質、 3 次収束。

パラメータ:

durandKernerAberth

template<typename T>
std::vector<Complex<T>> durandKernerAberth(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: Aberth-Ehrlich の全根同時反復 ($z_i \leftarrow z_i - p(z_i)/p'(z_i) / (1 - p(z_i)/p'(z_i) \cdot \sum_{j \ne i} 1/(z_i - z_j))$)。 収束後に元の多項式で各根を Newton polish 精密化する。 全根同時のため デフレーション誤差が蓄積しない。

パラメータ:

用途: 高次 (20 次以上) のデフォルト推奨。 solvePolynomial も 20 次以上で内部呼出。

bairstow

template<typename T>
std::vector<Complex<T>> bairstow(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: 試行 2 次因子 $x^2 + rx + s$ で多項式を割り、 Newton 法で $r, s$ を修正して 2 次因子を抽出。 実数演算のみで複素根 (共役ペア) を得る。 出典: Press et al., Numerical Recipes §9.5。

パラメータ:

用途: 複素数演算を避けたい組込み環境、 実係数のみで複素根の組をペア単位で扱いたい場合。

weierstrass

template<typename T>
std::vector<Complex<T>> weierstrass(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: 全根を同時反復で更新する古典的な Durand-Kerner (Weierstrass) 法。 durandKernerAberth と異なり Aberth 補正を持たない素朴な同時反復 (2 次収束)。 更新式 $z_i \leftarrow z_i - \frac{P(z_i)}{\mathrm{lc}\cdot\prod_{j \neq i}(z_i - z_j)}$。 複素初期値 $(0.4+0.9i)^i$ から始め、 重根が無ければ全根 (複素共役対を含む) へ収束する。

パラメータ:

注意: 重根・近接根では分母 $\prod_{j \neq i}(z_i - z_j)$ がほぼゼロになり収束が遅い。 数値的頑健性を重視するなら durandKernerAberth (Aberth 補正版) を使う。

// 例: (x-2)(x^2+1) = x^3 - 2x^2 + x - 2 → 2, ±i
Polynomial<double> p({-2, 1, -2, 1});
auto roots = weierstrass(p);
// 実行結果: roots = { 2, i, -i }

graeffe

template<typename T>
std::vector<Complex<T>> graeffe(
    const Polynomial<T>& poly,
    int maxSquarings = 6,
    T eps = std::numeric_limits<T>::epsilon() * T(100));

動作: Graeffe 根二乗法。 $p(x)\,p(-x) = s(x^2)$ を繰り返し作って根を $r_i \to r_i^{2^m}$ と二乗し、 モジュラスが分離した根のモニック係数比から $|r_i| \approx \left|\frac{a_{n-i}}{a_{n-i+1}}\right|^{1/2^m}$ でモジュラスを読む。 実根は符号を $P(\pm|r_i|)$ で復元し、 最後に Newton 法で数ステップ磨く。 係数は二重指数的に増大するため、 各二乗でモニック化し overflow 直前で二乗を打ち切る。

パラメータ:

注意 (適用範囲): 実根・モジュラスが分離した場合向け。 等モジュラス根・複素共役対・重根は基本 Graeffe では分離できない (歴史的手法)。 一般の複素根には jenkinsTraub / durandKernerAberth / weierstrass を使う。

// 例: (x-1)(x-2)(x-4) = x^3 - 7x^2 + 14x - 8 → 4, 2, 1
Polynomial<double> p({-8, 14, -7, 1});
auto roots = graeffe(p);
// 実行結果 (モジュラス降順): roots = { 4, 2, 1 }

lehmerSchur

template<typename T>
std::vector<Complex<T>> lehmerSchur(
    const Polynomial<T>& poly,
    T tol = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxDescend = 200);

動作: 全根を含む円板 (Cauchy 限界) を半径半分の小円板 9 枚 (中心 + 周囲 8) で被覆し、 各小円板に Schur-Cohn 判定 (Lehmer の単位円板テスト) を掛けて根を含む小円板へ降下する。 半径が十分小さくなったら中心を種に Newton 磨き → デフレーションし、 次数回繰り返して全根を得る。 判定は厳密な Schur 変換 $T q = \overline{a_0}\,q - a_n\,q^*$ を次数が下がるまで反復する boolean テスト (途中で $|a_0|^2 - |a_n|^2 < 0$ なら内部に根あり)。 種の精度は分割が、 根の精度は Newton 磨きが担うため、 複素根 (単位円上を含む) も求まる。

パラメータ:

用途: 根の局在化を幾何的に行いたい場合、 単位円上・近接する複素根を確実に分離したい場合。 数値的に最も安定な汎用解には jenkinsTraub を使う。

// 例: x^4 - 1 = 0 → 1, -1, i, -i （単位円上の複素根）
Polynomial<double> p({-1, 0, 0, 0, 1});
auto roots = lehmerSchur(p);
// 実行結果: roots = { 1, -1, i, -i }

統合ディスパッチャ

template<typename T>
std::vector<Complex<T>> solvePolynomial(
    const Polynomial<T>& p,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

solvePolynomial は次数に応じて自動的に最適なアルゴリズムを選択する:

• 1 次: solveLinear (閉形式)
• 2 次: solveQuadratic (閉形式、 桁落ち対策付)
• 3 次: solveCubic (Cardano + 三角関数)
• 4 次: solveQuartic (Ferrari)
• 5-19 次: jenkinsTraub (1 根ずつ + デフレーション、 高精度)
• 20 次以上: durandKernerAberth (全根同時 + Newton polish、 デフレーション誤差回避)

N-D: 多変数求根

$\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$ の零点を求める。 戻り値は RootFindingResult<V> (1D 版とテンプレートを共有する統合型、 V はベクトル型)。

関数詳細

newton_raphson_nd

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
RootFindingResult<V> newton_raphson_nd(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 各反復で $\mathbf{x}_{n+1} = \mathbf{x}_n - J(\mathbf{x}_n)^{-1} F(\mathbf{x}_n)$。 ヤコビアン特異時は対角要素を abs_xtol * 10 で正則化して solve をリトライ。

パラメータ:

返り値: RootFindingResult<V>。 root は std::optional<V>。

収束判定: $\|F(\mathbf{x})\|_2 < \mathrm{abs\_ftol}$ または $\|\Delta\mathbf{x}\|_2$ が許容値以下。

例: 単位円 $x^2 + y^2 = 1$ と直線 $y = x$ の交点を求める。 sangi::Matrix<T>::solve() を経由して concept 要件を満たすため <math/linalg/solvers.hpp> も include する。

#include <math/roots/root_finding.hpp>
#include <math/linalg/solvers.hpp>   // Matrix::solve() の本体定義はここ

using namespace sangi;

auto F = [](const Vector<double>& v) {
    Vector<double> r(2);
    r[0] = v[0]*v[0] + v[1]*v[1] - 1.0;   // x^2 + y^2 - 1
    r[1] = v[0] - v[1];                   // x - y
    return r;
};
auto J = [](const Vector<double>& v) {
    Matrix<double> m(2, 2);
    m(0,0) = 2*v[0]; m(0,1) = 2*v[1];
    m(1,0) = 1.0;    m(1,1) = -1.0;
    return m;
};
Vector<double> x0(2); x0[0] = 0.5; x0[1] = 0.5;
auto r = newton_raphson_nd<double, Vector<double>, Matrix<double>>(F, J, x0);
if (r.converged && r.root) {
    const auto& v = *r.root;
    std::cout << v[0] << ", " << v[1];
}
// 実行結果: iter=5  root=(0.707106781186548, 0.707106781186548)
//   ※ 真値 1/√2 = 0.7071067811865476 と 15 桁一致

broyden_method

template<concepts::OrderedField T, typename V, typename M>
RootFindingResult<V> broyden_method(
    const std::function<V(const V&)>& F,
    const V& x0,
    std::optional<M> J0 = std::nullopt,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 「Good Broyden」式の準 Newton 法。 ヤコビアンを毎回計算せず、 ランク 1 更新で近似 ($B_{n+1} = B_n + \frac{(\Delta F - B_n \Delta \mathbf{x}) \Delta \mathbf{x}^\top}{\|\Delta \mathbf{x}\|^2}$)。

パラメータ:

用途: ヤコビアンが解析的に書けない / 計算コストが高い場合。 初期 $B_0$ は単位行列、 後は更新で改善。

auto r = broyden_method<double>(F, x0);  // J 不要、 F のみ与える

newton_krylov

template<concepts::OrderedField T>
RootFindingResult<Vector<T>> newton_krylov(
    const std::function<Vector<T>(const Vector<T>&)>& F,
    const Vector<T>& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>(),
    size_t krylov_restart = 30,
    T inner_tol = T(1e-3));

動作: Jacobian-free Newton-Krylov (JFNK)。 ユーザは $F$ のみを渡す (解析 Jacobian も数値 Jacobian の陽な構築も不要)。 各 Newton ステップで $\mathbf{J}\,\boldsymbol\delta = -\mathbf{F}$ を GMRES (Krylov) で解くが、 Jacobian を行列として作らず方向微分の有限差分 $\mathbf{J}\,\mathbf{v} \approx (\mathbf{F}(\mathbf{x}+\varepsilon\mathbf{v}) - \mathbf{F}(\mathbf{x}))/\varepsilon$ だけで進める (matrix-free)。 バックトラッキング減衰で大域収束性を補う。

パラメータ:

用途: Jacobian を陽に組めない / 大きすぎる大規模非線形系。 broyden_method も $F$ のみで済むが、 こちらは Krylov 部分空間で各ステップの線形系を解くため、 高次元でメモリ効率がよい。

// x² + y² = 1, x - y = 0 → (1/√2, 1/√2)   ※ F だけ渡せばよい
using V = Vector<double>;
auto F = [](const V& v) -> V {
    return V({v[0]*v[0] + v[1]*v[1] - 1.0, v[0] - v[1]});
};
auto r = newton_krylov<double>(F, V({0.5, 0.5}));
// 実行結果: x = 0.7071067811865476, y = 0.7071067811865476

newton_fd

template<concepts::OrderedField T>
RootFindingResult<Vector<T>> newton_fd(
    const std::function<Vector<T>(const Vector<T>&)>& F,
    const Vector<T>& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>(),
    T fd_step = std::sqrt(std::numeric_limits<T>::epsilon()));

動作: 数値ヤコビアン Newton 法。 ユーザは $F$ のみを渡す。 各 Newton ステップでヤコビアンを前進差分 $J[:,j] = (F(x + h_j e_j) - F(x))/h_j$ で 陽な n×n 行列 (sangi の Matrix) として構築し、 $\mathbf{J}\,\boldsymbol\delta = -\mathbf{F}$ を直接解法 (LU) で厳密に解く。 中身は newton_raphson_nd に数値ヤコビアンを与える薄いラッパ。

パラメータ:

newton_krylov との違い: newton_fd は $J$ を陽に持ち (メモリ $O(n^2)$) 線形系を直接解くため小〜中規模で頑健。 大規模で $J$ を組めない場合は newton_krylov (matrix-free) を使う。

// x² + y² = 1, x - y = 0 → (1/√2, 1/√2)   ※ F だけ渡せばよい (数値 J を内部構築)
using V = Vector<double>;
auto F = [](const V& v) -> V {
    return V({v[0]*v[0] + v[1]*v[1] - 1.0, v[0] - v[1]});
};
auto r = newton_fd<double>(F, V({0.5, 0.5}));
// 実行結果: x = 0.7071067811865476, y = 0.7071067811865476

levenberg_marquardt

template<concepts::OrderedField T, typename V, typename M>
RootFindingResult<V> levenberg_marquardt(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: $(J^\top J + \lambda I) \Delta\mathbf{x} = -J^\top F$ を解く。 ダンピングパラメータ $\lambda$ を成功時減らし、 失敗時増やす (Marquardt-Levenberg 更新規則)。 Newton と最急降下の補間。

パラメータ:

用途: 初期値が真の根から遠い / ヤコビアンが特異近傍を通る / 最小二乗の感覚で根を探したい場合。 単なる方程式系解法というより、 残差最小化の混成手法。

powell_hybrid

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
RootFindingResult<V> powell_hybrid(
    const std::function<V(const V&)>& F,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Powell (1970) の dogleg-trust-region。 信頼領域内で Newton 方向と最急降下方向を組合せた dogleg 経路をたどる。 ヤコビアンは内部で差分近似するため ユーザは F のみ与えればよい。 MINPACK hybrd ルーチン (SciPy fsolve の標準) と同系統。

パラメータ:

用途: N-D 求根のデフォルト推奨。 局所収束性 (Newton 並み) と大域収束性 (信頼領域) を兼ね備える。

使用例

#include <math/roots/root_finding.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // 1D: f(x) = x^2 - 2 の根を Brent 法で求める
    auto result = brent_method<double>(
        [](double x) { return x * x - 2.0; },
        1.0, 2.0);
    if (result.converged && result.root)
        std::cout << "sqrt(2) = " << *result.root << std::endl;
    // 出力: sqrt(2) = 1.41421...

    // 多項式: x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3) の全根
    Polynomial<double> p({-6.0, 11.0, -6.0, 1.0});  // 昇べき
    auto roots = solvePolynomial(p);
    for (auto& r : roots)
        std::cout << r.re << std::endl;
    // 出力: 1, 2, 3 (順序は実装依存)

    // 5次以上: solvePolynomial が自動的に jenkinsTraub を呼ぶ
    Polynomial<double> q({1.0, 0.0, 0.0, 0.0, 0.0, 1.0});  // x^5 + 1
    auto roots5 = solvePolynomial(q);
    for (auto& r : roots5)
        std::cout << r.re << " + " << r.im << "i" << std::endl;
}

関連する数学的背景

以下の記事では、Roots モジュールの基盤となる数学的概念を解説している。

• 二分法 — 最も基本的な求根アルゴリズム
• ニュートン法 — 二次収束する反復法
• 割線法 — 導関数不要の超一次収束法
• 高度な多項式求根法 — Aberth-Ehrlich・companion matrix 法
• N 次根の精度倍増 Newton 法 — 多倍長 N 乗根
• Zimmermann 再帰平方根 — 高速多倍長平方根