【Rustで作る電卓】Remezのアルゴリズムでlog計算を最適化 #17

これまでにTUI電卓開発では、350年前の数学的知恵を借り、エストリン法でモダンなCPUの力を使えるように改良検討などしてきました。しかし、まだ解決していない「最適化問題」があります。

それは、「理論上の係数（1/3, 1/5…）が、浮動小数点の世界でも最良である」という思い込みです。

今回は、数値計算のプロたちが必ず通る道、Remez（レメッツ）のアルゴリズムを導入し、標準ライブラリに肉薄する「最強の対数関数」を完成させます。

1. Remezのアルゴリズム

Remezのアルゴリズムが目指すのは、数学的には「Minimax（最小最大）問題」の解決です。
具体的には、以下の式を満たすような多項式 \(P(x)\) の係数を見つけ出す作業を指します。
\[\min_{P \in \Pi_n} \left( \max_{a \le x \le b} |f(x) – P(x)| \right)\]

1. 1 グレゴリー級数で展開した式を見直す

ここで、前回までに使用してきたグレゴリー級数の式を振り返ってみましょう。
\[\ln(m) = 2 \left( z + \frac{z^3}{3} + \frac{z^5}{5} + \dots \right) ・・・[1]\]
この式の係数（\(1, 1/3, 1/5 \dots\)）は、数学的には極めて美しく、中心点（\(z=0\)）において完璧な精度を誇ります。しかし、計算機で実行する際には一つ問題があります。私たちが知りたいのは「\(z=0\) での完璧な値」ではなく、「対象とする全範囲（例：\(1 \le m < 2\)）において、どこを切り取っても十分な精度が出ること」だからです。

ここで先ほどの Remezのアルゴリズム（Minimax問題） が真価を発揮します。

そもそも なぜ lnのテイラー展開が上記の式[1]になっているのかは以下の記事で解説しています。
tui関数電卓 #15 : 【Rustで作る電卓】浮動小数点数の構造をハックして対数関数(log)を実装する #15
解説記事　なぜわざわざ ｚに置換したの？： ln自然対数の計算プログラム実装でzに置換する理由【ホーナー法・グレゴリー級数】

1. マクローリン展開の「甘え」を断つ

マクローリン展開（およびグレゴリー級数）は、特定の点（\(x=0\)など）の周りでは極めて正確ですが、そこから離れるにつれて急速に精度が悪化します。

これを「誤差の偏り」と呼びます。中心部は過剰に正確なのに、端の方では計算機が許容できる誤差を超えてしまう。限られた「次数」というリソースを、中心点だけに浪費している状態です。

[Figure 1] : テイラー展開の誤差

2. Remezのアルゴリズム：最小最大近似（Minimax）

ここで登場するのが、1.項で紹介したソビエト連邦の数学者エフゲニー・レメッツが提唱したRemezのアルゴリズムです。
このアルゴリズムの目的は、指定した範囲内において、「最大誤差を最小にする（Minimax）」ような多項式の係数を見つけることです。

(アルゴリズムを調べてみると、ソ連の数学者による突飛かつ最適化された実装（三進法コンピュータの「トリット」）など、今ではそこまで活発じゃないけど面白い分野が見つかり、非常に興味深いです。ロシア語はまだ何となくしか分かりませんが、いつか原典をちゃんと読めるようになったらもっと面白そうだなと思っています。)

特徴：等波特性（Equiripple Property）

[Figure 2] : 純粋なテイラー展開とRemezの誤差比較

Remezによって導かれた多項式の誤差グラフを見ると、まるで波のように、誤差が範囲内で均等にプラスマイナスへ振れます。係数をあらかじめ計算しておくだけでこの精度を出せるのがRemezのアルゴリズムの強みです。

「特定の値だけの最適化」をやめ、「全体が平均以上の成績」を目指す。

これが、現代の math.h や Rust の標準ライブラリが採用している戦略になります。

3. 実践：係数を「マジックナンバー」に書き換える

3.1 係数の「再チューニング」

Remezを適用するということは、式の形（奇数次の項の積み重ね）は維持したまま、その係数を理論値から「あえて僅かにずらす」ことを意味します。
\[P(z) = c_0 z + c_1 z^3 + c_2 z^5 + \dots\]
この\(c_0, c_1, c_2 \dots\)を、理論上の

\[1, 1/3, 1/5,\dots\]

ではなく、「指定した範囲内での最大誤差を最小化するマジックナンバー」へと最適化するのです。

なぜ誤差が減るのか？

テイラー展開が「中心点付近の精度を極限まで高め、遠くの精度を犠牲にする」のに対し、Remezは「中心の精度を少しだけ削り、その分を端の精度の補填に回す」という全体最適を行います。
その結果、次数（計算コスト）を増やすことなく、実効的な精度（最大誤差の少なさ）を劇的に向上させることが可能になります。まさに「理論」を「実装」へとコンパイルするための、数値計算における究極の調整作業といえると思います。

計算機にとっての「1/3」は 1/3 ではない？

例えば、理論上の係数が 0.3333333333333333（1/3）だとしても、Remezで最適化すると 0.33333331… のような、微妙にズレた値になります。

// 従来の理論値ではなく、Remezで算出した「最良」の係数を採用
const COEFFS_REMEZ: [f64; 6] = [
    1.0000000000000000,   // c0
    0.3333333333333123,   // c1 (理論値は 1/3)
    0.2000000000021456,   // c2 (理論値は 1/5)
    // ... 
];

4. 12次のln近似のRemez係数を求めてみた

Remezの係数計算をフルスクラッチするとそれだけで１シリーズ書けそうなくらい奥深いので、今回はpythonのライブラリの恩恵を享受しようと思います。

// 算出したRemez係数
const REMEZ_COEFFS: [f64; 12] = [
    0.9999999999999999, // c0
    0.3333333333333317, // c1
    0.2000000000001053, // c2
    0.1428571428522026, // c3
    0.1111111112446774, // c4
    0.0909090886161822, // c5
    0.0769231268962635, // c6
    0.0666663244247854, // c7
    0.0588258284534882, // c8
    0.0526189578278274, // c9
    0.0477174665487779, // c10
    0.0428456801269550, // c11
];

5. 実験：自作 ln vs std::f64::ln

ホーナー法による高速な計算に、Remezによる最強の係数を搭載。これで,フルスクラッチで作れるln()としての理論値に近いものが完成しました。

結果の考察

実験結果より、Remez法を用いることで、0.5ms程度高速化し、遅くとも約4ms程度で１万回の計算ができるようになった。Horner法だけの場合と比べて優位に高速化したといえる。また、Std::ln()と比べると、遅くなることもあるが、安定して同じ速度で計算できるというのは計算機のアルゴリズムとしては良いといえる。
最終的な実装は horner法+ Remezのアルゴリズムでの最適化を用いた方法を選択することにした。

結び：理論と実装の間に

全三回にわたって「対数関数の再実装」をおこなってきました。いかがでしたでしょうか？

1. 理論編： グレゴリー級数という350年前の知恵を知る。
2. 構造編： エストリン法でCPUの並列性を引き出す。
3. 近似編： Remez法で計算機にとっての「最良」を定義する。

私たちが何気なく使っている関数一つ一つに、これだけの歴史と工夫が詰まっています。「ただ動く」のその先にある、数値計算の深淵を少しでも楽しんでいただけたなら幸いです。

次は、さらなる困難が予想される（？）——sinの実装三角関数の範囲還元(周期性を生かす)でお会いしましょう。

三角関数の範囲還元 :　関連記事は、2026年5月26日に公開予定 (あと2時間)

ここまで読んでいただき、ありがとうございました。

関連記事

第11回　平方根の計算実装：【Rust自作TUI電卓】mod.rsを活用したスマートなモジュール管理と基本数学関数の実装 #11
第13回　n乗根の計算実装：【Rust自作TUI電卓】ニュートン法×マジックナンバーで実装する汎用n乗根 #13
第14回　ハレー法による立方根：【Rust自作TUI電卓】ニュートン法を超えろ！ハレー法による立方根(cbrt)の極限最適化 #14
第15回　logの計算基本編：【Rustで作る電卓】浮動小数点数の構造をハックして対数関数(log)を実装する #15
バビロニア法・ニュートン法：平方根の導出展開、手動でも解ける差分アプローチ【バビロニア法・ニュートン法】
理論編：lnの近似 グレゴリー級数：【グレゴリー・ライプニッツ級数】微積分学が確立される前の無限への挑戦を紐解く
解説 lnの近似、なぜわざわざ ｚに置換したの？： ln自然対数の計算プログラム実装でzに置換する理由【ホーナー法・グレゴリー級数】
構造編：解説 ホーナー法：ホーナー法とは？多項式計算の効率化手法そのメリットとデメリット
構造編：解説 エストリン法：【多項式計算の最適化】ホーナー法と、並列処理に強い「エストリン法」の比較検証