correct rounding への道 (5) table-maker’s dilemma

やりましょう。

前回、double-double や triple-double の話をしましたが、「triple-double は必要なのか？」とか「quadruple-double は不要なのか？」とかについて触れた方がいいかなと思ったので、先にそういう話をします。重要なお話です。

disclaimer: 重要なお話ではあるものの、計算するにあたってのアルゴリズムなどの話ではないです。

• 記法
• 導入
  • 例
  • 用語について
• 本題？ 雑談？
• 参考文献
• おわり

記法

多くの記法については前回定義したものを使います。特に、double で表せる浮動小数点数全体からなる集合を $F$ とします。

今回以降、実数の表記として主に指数表記つきの 16 進法や 2 進法を使っていきます。たとえば $$ \pi = {\small 3.14159265358979}{\footnotesize 323846264338327}{\scriptsize 950288419716939}{\tiny 937510582097494{\ldots}} $$ ではなく $$ \begin{aligned} \pi &= \texttt{1}.\texttt{{\small 921FB54442D18}{\footnotesize 469898CC51701}{\scriptsize B839A252049C1}{\tiny 114CF98E80417{\ldots}}}_{(16)}\times 2^{1} \\ &= 1. \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1011}_{\texttt{B}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1101}_{\texttt{D}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0110}_{\texttt{6}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 1011}_{\texttt{B}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1110}_{\texttt{E}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0001}_{\texttt{1}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} {\footnotesize \,\dots\,}_{(2)} \times 2^1 \end{aligned} $$ のように書きます。指数の底は $2$ です。2 進法で書く場合は 4 桁ごとに区切って 16 進法での値を添えるようにします。

言語によっては、浮動小数点型の 16 進法リテラルがサポートされていて、0x1.921FB54442D18p+1 のように書くこともできます。

C ではリテラルとして可能で、Python でも float.fromhex() を使うことで可能です。Rust では標準ではできず、hexfloat2 などのクレートを使う必要があります*1。また、出力に関しても、C では printf("%a", x) としたり、Python では x.hex() とすることで可能です。Rust はお察しです。

導入

例

$\sin(x)$ を計算したくなったことにしましょう。

ということで、次の手続き $\textsc{Fr54-Sin}$ を考えます*2。

• 任意の $x\in F\smallsetminus\{-\infty, 0, +\infty\}$ を入力とする。
• 下記を満たす $y\in \Q$ を出力する。
  • ある組 $(s, e, m) \in \{0, 1\}\times\Z\times([2^{53}\lldot 2^{54})\cap\N)$ が一意に存在して ${y = (-1)^s\cdot m\cdot 2^{e-53}}$ が成り立つ。
  • 上記の $e$ に対して $|y-\sin(x)|\le 2^{e-53}$ が成り立つ。

$\roundcirc{\sin(x)} = \roundcirc{\textsc{Fr54-Sin}(x)}$ として求めようというわけです。

たとえば $x = \texttt{1}.\texttt{D037CB27EE6DF}_{(16)}\times 2^{-3}\in F$ を考えます。 実際の計算結果は下記の通りです。 $$ \begin{aligned} &\phantom{{}={}} \textsc{Fr54-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 1\hphantom{000}}_{\texttt{8}} {}_{(2)} \times 2^{-3} \end{aligned} $$ 最後の部分について、$\underbrace{1010}_{\texttt{A}}\,100\ldots01\ldots$ と $\underbrace{1010}_{\texttt{A}}\,011\ldots10\ldots$ のどちらを丸めた結果なのかを判別できません。前者であれば、$\roundcirc{\sin(x)} = \texttt{1}.\texttt{CC40C3805229B}_{(16)}\times 2^{-3}$ ですし、後者であれば $\roundcirc{\sin(x)} = \texttt{1}.\texttt{CC40C3805229A}_{(16)}\times 2^{-3}$ なので、$\textsc{Fr54-Sin}(x)$ を用いて求めることはできません。

note: $x\in\Q\smallsetminus\{0\}$ より $\sin(x)\notin\Q$ なので、$\sin(x) = \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3}$ となることはありません。

ということで、もう少し高い精度で計算してみましょう。同様に $\textsc{Fr55-Sin}(x)$ を定義して計算します。 $$ \begin{aligned} &\phantom{{}={}} \textsc{Fr55-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 10\hphantom{00}}_{\texttt{8}} {}_{(2)} \times 2^{-3}. \end{aligned} $$ まだ同様の理由でだめですね。さらに高い精度でやりましょう。

〜しばらく繰り返し〜

$$ \begin{aligned} &\phantom{{}={}} \textsc{Fr105-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A8000000000000}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} \\ & \qquad \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} {}_{(2)} \times 2^{-3}. \end{aligned} $$ ま、まだだめです...。

$$ \begin{aligned} &\phantom{{}={}} \textsc{Fr106-Sin}(x) \\ &= \texttt{1}.\texttt{CC40C3805229A7FFFFFFFFFFFF8}_{(16)}\times 2^{-3} \\ &= 1. \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0100}_{\texttt{4}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0011}_{\texttt{3}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1010}_{\texttt{A}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1\hphantom{000}}_{\texttt{8}} {\footnotesize\,} {}_{(2)} \times 2^{-3}. \end{aligned} $$ わ、ようやく区別がつきました。 これにより $$ \sin(x) \lt \texttt{1}.\texttt{CC40C3805229A8}_{(16)}\times 2^{-3} $$ が言えるので、 $$ \roundcirc{\sin(x)} = \texttt{1}.\texttt{CC40C3805229A}_{(16)}\times 2^{-3} $$ であるとわかったことになります。めでたしめでたし*3。実際のところ、 $$ \sin(x) = \texttt{1}.\texttt{{\small CC40C3805229A}{\footnotesize 7FFFFFFFFFFFF}{\scriptsize 83E76BBF0FCBE}{\tiny 19AFB5AC38007{\ldots}}} {\,}_{(16)} \times 2^{-3} $$ です。

用語について

さて、上記の例では 106-bit での計算で十分でしたが、これが「106-bit で十分ですよ」ということを予め知っておくことは可能でしょうか？ 平方根など一部の例では可能ですが、三角関数を始めとする多くの関数ではこれは難しいということになっています。

William M. Kahan*4は次のように記しています^[6]。

Why can’t Y^W be rounded within half an ulp like SQRT ? Because nobody knows how much computation it would cost to resolve what I long ago christened “The Table-Maker’s Dilemma”.

拙訳：なぜ $Y^W$ は、平方根のように 0.5 ULP 以内（の誤差で収まるよう）に丸められないのか？ それは、以前私が 数表作成者のジレンマ (Table-Maker’s Dilemma) と名づけた事象を解決するためにどれだけの計算が必要かを誰も知らないからだ。

note: table-maker’s dilemma は、しばしば TMD と略されます。

この後、$\tfrac18\cosh(\pi\sqrt{163}) - (2^{53}-1)$ を 53 ビットに正確に丸めた値を計算しようとして、精度を上げつつ ${\small 7401389035307055}.{\footnotesize 49999978}$ や ${\small 7401389035307055}.{\footnotesize 5000000000015}$ を得る話が続きます。先ほどの $\textsc{Fr106-Sin}(x)$ の例と同じような感じです。

No general way exists to predict how many extra digits will have to be carried to compute a transcendental expression and round it correctly to some preassigned number of digits.

拙訳：超越式（超越関数を含む式？）を計算して予め決められた桁数に正確に丸めるにあたり、どれだけの桁数を余分に持っておく必要があるかを予測する一般的な方法は存在しない。

Table と言っているのは、たぶん教科書の巻末に載っている三角関数表のようなもののことだと思います。下記のようなイメージですね。

$x$	$\roundcirc{\sin(x)}$
$\texttt{1}.\texttt{D037CB2700000}_{(16)}\times 2^{-3}$	$\texttt{1}.\texttt{CC40C37F69D51}_{(16)}\times 2^{-3}$
$\vdots$	$\vdots$
$\texttt{1}.\texttt{D037CB27EE6DF}_{(16)}\times 2^{-3}$	$\texttt{1}.\texttt{CC40C3805229A}_{(16)}\times 2^{-3}$
$\vdots$	$\vdots$
$\texttt{1}.\texttt{D037CB2800000}_{(16)}\times 2^{-3}$	$\texttt{1}.\texttt{CC40C3806348B}_{(16)}\times 2^{-3}$

これの各行について、正確に丸めた値を計算するにあたってどれだけの精度が必要なのかを予測できず、徐々に桁数を上げながら繰り返す*5にしてもどの程度待てばいいかわからないねえという感じでしょう。

たとえば、たかが 53-bit での correctly-rounded な値を得るための途中計算に $10^{3000}$ bits の精度が必要ないとどうして言い切れるでしょうか？ 実際、先の $\sin(x)$ のような例が存在することも予想していなかった人が多いのではないかと思います。 下記のような例もあります。 $$ \begin{aligned} &\phantom{{}={}} \tan(\texttt{1}.\texttt{50486B2F87014}_{(16)}\times 2^{-5}) \\ &= 1. \underbrace{\footnotesize 0101}_{\texttt{5}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 1110}_{\texttt{E}} {\footnotesize\,} \underbrace{\footnotesize 1011}_{\texttt{B}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1111}_{\texttt{F}} {\footnotesize\,} \underbrace{\footnotesize 1001}_{\texttt{9}} {\footnotesize\,} \underbrace{\footnotesize 1100}_{\texttt{C}} {\footnotesize\,} \underbrace{\footnotesize 0111}_{\texttt{7}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} \underbrace{\footnotesize 1000}_{\texttt{8}} {\footnotesize\,} \\ & \qquad \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0000}_{\texttt{0}} {\footnotesize\,} \underbrace{\footnotesize 0010}_{\texttt{2}} {\footnotesize\,} {\footnotesize \, \ldots} {\,}_{(2)} \times 2^{-5}. \end{aligned} $$

本題？ 雑談？

double-double だけだと 106-bit 弱になってしまうので、先のような例では負けてしまいます。triple-double ではどうでしょうか？ 足りますか？ わかりません... いかがでしたか？

「100 台程度のワークステーションを使って 4 年かけて最悪ケースの探索をしました」という論文^[2]が 2000 年に出ていますが、下記のように書かれています。

The results [...] give the worst cases for functions $\sin$, $\arcsin$, $\cos$, $\arccos$, $\tan$ and $\arctan$. For those functions, we have worst cases in some bounded domain only, because trigonometric functions are more difficult to handle than the other functions.

拙訳・要約：三角関数は他の関数たちより扱いが難しく、いくつかの区間における最悪ケースしか得られておりません。

そんなあ。

ということで、次からの記事は log や exp の correct rounding な実装の話になるかもしれないし、三角関数関連の話をするかもしれないし、サーベイの旅に出たっきり帰ってこないかもしれません。

Kahan も [6] で触れていますが、Bessel 関数 $J_n(x)$ のような名前のついた非初等関数や、$y^w$ のように複数の引数がある関数についての最悪ケースの探索はほぼ望み薄でしょう。身近なところだと atan2 とかはどうなるんでしょうね。hypot は超越関数ではないのでそこまで悲惨なことにはならないかもしれません。

これはあくまで感情側の話なのですが、標準ライブラリの実装が「ま〜 1–10 ULP くらいの誤差に収まっている気がするしいいでしょ（未証明）」という感じだったり規格側が「別に精度については規定しないよ」という感じだったりするのは困ります。一方で、規格側が「（たとえば）3 ULP までの誤差は許容する」と言ってきても「その値の根拠はどこから...？」となりますし、そもそも「可搬性は...？」となります。値自体の可搬性はないけど事後条件に関しては担保するよ*6みたいな？ 可搬性を第一目標にするなら「この数学関数はこういうアルゴリズムで計算する」とでもすればよいですが、よりよい精度・速度のアルゴリズムが提案されたときにどうするの？という話にもなってきます。

というような旨のことが [7] の文献で書かれていて、次のような主張*7の節があります。

3.2 $\quad$ The only specification that makes sense is correct rounding

なかなか強めの思想な気もしないではないですが、まぁ否定もできないかなという感じです。

といった感じで、三角関数に関してはちゃんとした最悪ケースがまだわかっていないと思っているのですが、correct rounding を謳うライブラリたちはどうやっているのでしょう。Ziv’s rounding test をしながら最大 768-bit くらいの精度まで試すみたいな話はどこかで見たような気はしつつ、それで十分である証明なしには correct と主張できないはずですし。

次数 $d\ge 2$ の代数的数 $\alpha$ に対しては、ある定数 $C_{\alpha}$ が存在して任意の整数 $p$, $q$ ($q\ge 1$) に対して $|\alpha-\tfrac pq| \ge \tfrac{C_{\alpha}}{q^d}$ が成り立つことが Liouville によって示されているらしいですが、超越数は...？という感じで、参っています。どうしましょう？ やはり枝刈り全探索するしか？

たとえば何らかのアルゴリズムによって高々 $2^{-768}$ 程度の相対誤差しかない値を得られるのであれば 53-bit の correct rounding にこだわる必要はないのでは？という主張は真っ当な気がしつつ、まぁそれはそれかなという気持ちでやっています。

参考文献

• [1] Vincent Lefèvre (1997). “An Algorithm that Computes a Lower Bound on the Distance Between a Segment and $\Z^2$.” ftp://ftp.lip.ens-lyon.fr/pub/LIP/Rapports/RR/RR97/RR97-18.ps.Z
• [2] Vincent Lefèvre, Jean-Michel Muller (2000). “Worst Cases for Correct Rounding of the Elementary Functions in Double Precision.” https://inria.hal.science/inria-00072594/document
• [3] Vincent Lefèvre (2000). “The Table Maker’s Dilemma.” http://www.vinc17.org/research/slides/aoc2000-11.pdf
• [4] Vincent Lefèvre, Jean-Michel Muller (2003). “The Table Maker's Dilemma: our search for worst cases.” https://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm
• [5] Florent de Dinechin, David Defour, Christophe Lauter (2004). “Fast correct rounding of elementary functions in double precision using double-extended arithmetic.” https://inria.hal.science/inria-00071446/file/RR2004-10.pdf
• [6] William M. Kahan (2004). “A Logarithm Too Clever by Half.” https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
• [7] Nicolas Brisebarre, Guillaume Hanrot, Jean-Michel Muller, Paul Zimmermann (2025). “Correctly-rounded evaluation of a function: why, how, and at what cost?” https://hal.science/hal-04474530/document

Vincent Lefèvre が大活躍すぎ。また、最後のサーベイ論文に関しては 2025 年（初版は 2024 年）で、びっくりです。

おわり

ということでえびちゃんはこれからお勉強です。フォロヮもごきげんよう。

おわりです。