イテレータとしての箱ティッシュとハンドソープ

箱ティッシュとハンドソープは、それぞれ「次の要素を取得する」という操作ができます。つまりイテレータですね。 しかし、一言でイテレータと言っても両者は少し性質が違います。

箱ティッシュは、要素の取得と同時に（ユーザから観測できる形で）次の要素がスタンバイされます。 あるいは、次の要素を取得しようとする前に次の要素の有無がわかります。

tissue_box = TissueBox()  # not empty
while True:
    tissue = tissue_box.next()
    # do some work
    if tissue_box.the_last_element_was_emitted():
        break

一方でハンドソープは、要素を取得しようとした時点でようやく次の要素の有無がわかります。

soap_dispenser = SoapDispenser()  # not empty
while True:
    soap = soap_dispenser.next()
    if soap is None:
        break
    # do some work

たとえば Rust でイテレータを使う際の構文は、ハンドソープの方の気持ちに近いような気がします。

soap_dispenser = SoapDispenser::new()
while let Some(soap) = soap_dispenser.next() {
    // do some work
}

内部的な実装はそうとは限らず、箱ティッシュのように次の要素をスタンバイさせている状況もあります。 <std::iter::successors as Iterator>::next() などは、次のような実装になっています。

fn next(&mut self) -> Option<Self::Item> {
    let item = self.next.take()?;
    self.next = (self.succ)(&item);
    Some(item)
}

返してしまう要素は所有権を渡してしまっているので、それを使って self.succ を呼ぶことはできないですからね。 これにより、「$n+1$ 個目の要素を取得しようとすると何らかの事情で panic するが、実際には $n$ 個目までの要素で終了する」というときに、期待に反して $n+1$ 個目の要素が計算されて panic したりします。

fn main() {
    for i in std::iter::successors(Some(1), |&x| Some(x * 1000)) {
        println!("{i}");
        if i == 10_u32.pow(9) {
            break;
        }
    }
}

これは、1 1000 1000000 を出力した後に panic します (playground)。

実際には箱ティッシュは exact な .the_last_element_was_emitted() は実装されているわけではなく、.probably_the_last_element_was_emitted() であることも多いですし、どちらのイテレータも .next() の呼び出し時に「もうすぐなくなりそうだな」感のある値も暗に返ってきている感はありますが。

買い替えなどのことを考えると、やっぱり箱ティッシュ型のイテレータの方が助かりますね。とはいえハンドソープ（やシャンプーなど、その他多くのハンドソープ型イテレータのインタフェースになっている製品）が箱ティッシュ型のイテレータになると、それはそれで不幸になりそうな気はします。

おわり

おわりです。

「○○○ー○としての○○○○○○と○○○○ー○」というタイトルでそのうち記事を書こうと思います

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025å¹´3月24æ—¥

correct rounding への道 (2) shift-and-add algorithm

前回 は、主に自分で考えた素朴なアルゴリズムを実装・証明しましょうというスタンスの回でした。 今回は（今後は？）既存の実装をなぞる形で勉強していきます。

紹介

実装

LLVM のものを参考にしています。

• libc/src/__support/FPUtil/generic/
• libc/src/math/generic/

資料

• Jean-Michel Muller. (2016). “Elementary Functions: Algorithms and Implementation.” Birkhäuser, Boston, MA. https://doi.org/10.1007/978-1-4899-7983-4
  • https://link.springer.com/book/10.1007/978-1-4899-7983-4
• Catherine Daramy-Loirat, David Defour, Florent de Dinechin, Matthieu Gallet, Nicolas Gast, et al. “CR-LIBM A library of correctly rounded elementary functions in double-precision.” [Research Report] LIP. 2006. ensl-01529804
  • https://ens-lyon.hal.science/ensl-01529804v1/document

ツール

• Gappa (Génération Automatique de Preuves de Propriétés Arithmétiques)
  • automatic generation of proofs of arithmetic properties
  • 誤差評価の証明などに用いる
• Sollya
  • ハードコード用の魔法の多項式を生成するのに用いる
  • 数値計算の用途でもしばしば使う
• SageMath
  • 諸々の数値計算などをしたいときに用いる
• SymPy
  • 数値以外計算をしたいときに用いる

ツールに関しては、今回はまだ使いません。

おわび

disclaimer: 魔法の多項式と Newton 法については一旦おあずけです。誤差評価がよくわかりません。$x, y\ge 1$ を踏まえると $y\oplus (x\oslash y)$ の部分の誤差は $y+\tfrac xy$ と比べて 1 ULP 程度に収まるのかな？

本題

さて、shift-and-add algorithm と呼ばれる典型テクをやっていきます。LLVM の sqrt でも使われているほか、Elementary Functions の 8 章でも紹介されています。

range reduction や reconstruction の部分は前回同様で、浮動小数点数 $x\in[1\lldot 4)$ について $\roundcirc{\sqrt x}$ を求めるとします。approximation の部分のアルゴリズムです。

考察

$$ \begin{aligned} y^{(n)} &= \sum_{i=0}^{n} y_i\times 2^{-i}, \\ r^{(n)} &= 2^n\cdot \left(x - (y^{(n)})^2\right) \end{aligned} $$ で定義します。ここで $y_i\in\{0, 1\}$ です。 初期値は $y^{(0)} = 1$ および $r^{(0)} = x-1$ です。

$n\ge 1$ とし、次のようにして漸化式を得ます。 $$ \begin{aligned} r^{(n)} &= 2^n\cdot\left(x - (y^{(n)})^2\right) \\ &= 2^n\cdot\left(x - (y^{(n-1)} + y_n\times 2^{-n})^2\right) \\ &= 2^n\cdot\left(x - \left( (y^{(n-1)})^2 + 2y^{(n-1)}\cdot y_n\times 2^{-n} + y_n^2\times 2^{-2n}\right)\right) \\ % &= 2^n\cdot\left(x - (y^{(n-1)})^2 - y^{(n-1)}\cdot y_n\cdot 2^{-(n-1)} - y_n^2\cdot 2^{-2n}\right) \\ &= 2\cdot 2^{n-1}\cdot\left(x - (y^{(n-1)})^2\right) - y_n\cdot\left(2y^{(n-1)} + y_n\cdot 2^{-n}\right) \\ &= 2r^{(n-1)} - y_n\cdot\left(2y^{(n-1)}+2^{-n}\right). \end{aligned} $$ 最後の等号は $y_n\in\{0, 1\}$ から従います*1。

ここで、$r^{(n)}\ge 0$ かつ $r^{(n)}$ が最小となるように $y^{(n)}$ を定めると、 $$ y_n = \begin{cases} 1, & \text{if}\:\: 2r^{(n-1)} \ge 2y^{(n-1)} + 2^{-n}; \\ 0, & \text{otherwise} \end{cases} $$ を得ます。

まず、$y^{(n)}$ や $r^{(n)}$ についての性質を示しておきましょう。

Claim 1: 任意の整数 $n\ge 0$ に対して、$y^{(n)}\le \sqrt{x}\lt y^{(n)}+2^{-n}$ が成り立つ。

Proof

数学的帰納法で示す。$P(n) \iff y^{(n)}\le \sqrt{x}\lt y^{(n)}+2^{-n}$ とする。

To-be-proved 1: $P(0)$

$y^{(0)} = 1$, $y^{(0)}+2^{-0} = 2$ および $1\le x\lt 4$ より従う。

To-be-proved 2: $P(n) \implies P(n+1)$

Case 1: $y_{n+1} = 0$

$y^{(n+1)} = y^{(n)} \le \sqrt x$ は明らか。$\sqrt x\lt y^{(n+1)}+2^{-(n+1)}$ を示す。

$y_{n+1}=0$ より、下記が成り立つ。

• $r^{(n+1)} = 2r^{(n)}$,
• $2r^{(n)}\lt 2y^{(n)}+2^{-(n+1)}$, and
• $y^{(n+1)} = y^{(n)}$.

すなわち、$r^{(n+1)}\lt 2y^{(n+1)}+2^{-(n+1)}$ が成り立つ。 $$ \begin{aligned} \left(y^{(n+1)}+2^{-(n+1)}\right)^2 &= (y^{(n+1)})^2+2^{-n}\cdot y^{(n+1)}+2^{-2(n+1)} \\ &= (y^{(n+1)})^2 + 2^{-(n+1)}\cdot(2y^{(n+1)}+2^{-(n+1)}) \\ &\gt (y^{(n+1)})^2 + 2^{-(n+1)}\cdot r^{(n+1)} \\ &= (y^{(n+1)})^2 + \left(x-(y^{(n+1)})^2\right) \\ &= x. \end{aligned} $$

Case 2: $y_{n+1} = 1$

$\sqrt x\lt y^{(n+1)}+2^{-(n+1)} = y^{(n)}+2^{-n}$ は明らか。$y^{(n+1)}\le \sqrt x$ を示す。

$y_{n+1} = 1$ より、下記が成り立つ。

• $r^{(n+1)} = 2r^{(n)}-(2y^{(n)}+2^{-(n+1)})$,
• $2r^{(n)}\ge 2y^{(n)}+2^{-(n+1)}$, and
• $y^{(n+1)} = y^{(n)}+2^{-(n+1)}$.

すなわち、$r^{(n+1)}\ge (2y^{(n)}+2^{-(n+1)})-(2y^{(n)}+2^{-(n+1)}) = 0$ が成り立つ。 よって、$r^{(n+1)} = 2^n\cdot(x-(y^{(n+1)})^2)$ より従う。$\qed$

Claim 2: 任意の整数 $n\ge 0$ に対して、$r^{(n)}\lt 4$ が成り立つ。

Proof

Claim 1 より、$\sqrt x-y^{(n)}\lt 2^{-n}$ が成り立つ。 よって、 $$ \begin{aligned} r^{(n)} &= 2^n\cdot(x-(y^{(n)})^2) \\ &= 2^n\cdot(\sqrt x-y^{(n)})(\sqrt x+y^{(n)}) \\ &\lt 2^n\cdot 2^{-n}\cdot(\sqrt 4+2) \\ &= 4.\quad\qed \end{aligned} $$

さて、$2r^{(n-1)} \ge 2y^{(n-1)} + 2^{-n}$ の判定について考えます。

$2y^{(n-1)}\in[2\lldot 4)$ が常に成り立つことから、$2^{-n}\ge 2^{-51}$、すなわち $n\lt 52$ であれば $2r^{(n-1)}\ge 2y^{(n-1)}\oplus 2^{-n}$ として計算できます。 よって、$n = 52$ についてのみ別途考えます。

Claim 3: 任意の浮動小数点数 $r\ge 0$ および $y\in[1\lldot 2)$ に対して、$r\ge y+2^{-53} \iff r\gt y$ が成り立つ。

Proof

($\implies$): 明らか。

($\impliedby$): 対偶 $r\lt y+2^{-53} \implies r\le y$ を示す。

Case 1: $r\lt 1$

$y\ge 1$ より従う。

Case 2: $r\ge 1$

整数 $m_r \ge 2^{52}$ および $m_y\ge 2^{52}$ を用いて $r-y = (m_r-m_y)\times 2^{-52}$ と表せる。 よって、$r-y\lt 2^{-53} \implies r-y\le 0$ となる。$\qed$

さて、$y_n = 1$ だった場合、$r^{(n)} = 2r^{(n-1)} - (2y^{(n-1)}+2^{-n})$ なので、この計算についても実行可能である必要があります。

Lemma 4: 任意の浮動小数点数 $x\in[1\lldot 4)$, $y\in[1\lldot 2]$ に対し、$x-y\in[0\lldot 2)$ ならば $x\ominus y = x-y$ となる。

Proof

Case 1: $x\in[1\lldot 2)$

$2y-x\ge 2\cdot 1-2=0$ より $x\le 2y$。$0\lt \tfrac y2\lt y\le x$ は明らかであり、Sterbenz lemma より従う。

Case 2: $x\in[2\lldot 4)$

note: たとえば $x=3-2^{-51}$, $y=1$ のとき、Sterbenz lemma の前提が成り立たない。

整数 $m_x\in[2^{52}\lldot 2^{53})$, $m_y\in[2^{52}\lldot 2^{53}]$ を用いて $x=m_x\times 2^{-51}$, $y=m_y\times 2^{-52}$ と表せる。

$$ \begin{aligned} x-y &= m_x\times 2^{-51} - m_y\times 2^{-52} \\ &= (2m_x-m_y)\times 2^{-52}. \end{aligned} $$ $x-y\in [0\lldot 2)$ より、$2m_x-m_y \in [0\lldot 2^{53})$ が成り立つ。$\qed$

Claim 5: 任意の整数 $0\le n\lt 52$ に対して、$y_n=1$ のとき $2r^{(n-1)} \ominus (2y^{(n-1)}\oplus 2^{-n}) = 2r^{(n-1)} - (2y^{(n-1)}+2^{-n})$ が成り立つ。

Proof

Claim 2 より $r^{(n)}\in [0\lldot 4)$ なので、$\tfrac12 r^{(n)}=r^{(n-1)}-(y^{(n-1)}+2^{-(n+1)}) \in[0\lldot 2)$ が成り立つ。 $y^{(n-1)}+2^{-(n+1)}\in[1\lldot 2]$ なので、$r^{(n-1)}\in[1\lldot 4)$ となる。

よって、Lemma 4 より従う。$\qed$

以上により、$y^{(0)}, \dots, y^{(52)}, y_0, \dots, y_{52}, r^{(0)}, \dots, r^{(51)}$ が得られます。 $r^{(52)}$ は得られていません。なんとかして $y_{53}$ を求めて、$y^{(52)}+2^{-52}\cdot y_{53}$ を返す必要があります。

note: 前回 の Claim 2 を踏まえつつ、$y_{53} = 1$ なら $y^{(52)} + 2^{-52}$、そうでなければ $y^{(52)}$ が答えです。

Claim 6: 任意の整数 $0\le n\le 52$ に対して、$r^{(n)}\equiv 0\pmod{2^{-52}}$ が成り立つ。

Proof

数学的帰納法で示す。$P(n) \iff r^{(n)} \equiv 0\pmod{2^{-52}}$ とする。

To-be-proved 1: $P(0)$

$x\in[1\lldot 4)$ より $x\equiv 0\pmod{2^{-52}}$ なので、$r^{(0)} = x-1\equiv 0\pmod{2^{-52}}$。

To-be-proved 2: $n\lt 52 \wedge P(n) \implies P(n+1)$

Case 1: $y_{n+1} = 0$

$$r^{(n+1)} = 2r^{(n)} \equiv 0\pmod{2^{-52}}.$$

Case 2: $y_{n+1} = 1$

$y^{(n)}\in[1\lldot 2)$ より $y^{(n)}\equiv 0\pmod{2^{-52}}$ なので、 $$ \begin{aligned} r^{(n+1)} &= 2r^{(n)}-(2y^{(n)}+2^{-(n+1)}) \\ &\equiv 0 \pmod{2^{-52}}. \quad\qed \end{aligned} $$

さて、$y_{53}$ を求めていきます。まだ計算できていない値たちの定義は下記です。 $$ \begin{aligned} r^{(52)} &= 2r^{(51)} - y_{52}\cdot(2y^{(51)}+2^{-52}), \\ y_{53} &= \begin{cases} 1, & \text{if}\:\: 2r^{(52)} \ge 2y^{(52)} + 2^{-53}; \\ 0, & \text{otherwise}. \end{cases} \end{aligned} $$ $y_{52} = 0$ のとき、$r^{(52)} = 2r^{(51)}$ および $y^{(52)} = y^{(51)}$ より、Claim 6 を踏まえて $$ \begin{aligned} y_{53} = 1 &\iff 2r^{(51)} - y^{(51)} \ge 2^{-54} \\ &\iff 2r^{(51)} - y^{(51)} \gt 0 \\ &\iff 2r^{(51)} \gt y^{(51)} \end{aligned} $$ です。$y_{52} = 1$ のときは、 $$ \begin{aligned} 2r^{(51)} &\ge 2y^{(51)}+2^{-52}, \\ r^{(52)} &= 2r^{(51)} - (2y^{(51)}+2^{-52}) \end{aligned} $$ であり、 $$ 2\cdot\left(2r^{(51)} - (2y^{(51)}+2^{-52})\right) \ge 2y^{(52)}+2^{-53} $$ すなわち $$ r^{(51)}-y^{(51)} \ge 2^{-1}\cdot y^{(52)} + 2^{-53} + 2^{-55} $$ を判定すればよいです。

Lemma 7: $y_{52} = 1$ のとき、$r^{(51)}\ominus y^{(51)} = r^{(51)}-y^{(51)}$ が成り立つ。

Proof

Claim 5 とほぼ同様にして示す。

$r^{(51)}-(y^{(51)}+2^{-53})\in[0\lldot 2)$ と $y^{(51)}+2^{-53}\in[1\lldot 2)$ より、$r^{(51)}\in [1\lldot 4)$ が成り立つ。 また、$r^{(52)}\lt 4$ より $r^{(52)}\le 4-2^{-51}$ なので、$r^{(51)}-y^{(51)} = \tfrac12 r^{(52)}+2^{-53}\in [0\lldot 2)$ である。 よって、Lemma 4 より $r^{(51)}\ominus y^{(51)} = r^{(51)}-y^{(51)}$ となる。$\qed$

Claim 8: $y_{52} = 1$ のとき、下記が成り立つ。 $$ y_{53} = 1 \iff r^{(51)}\ominus y^{(51)} \gt 0.5\otimes(y^{(52)} \oplus 2^{-52}). $$

Proof

$y^{(52)}\in[1\lldot 2)$ より $y^{(52)}\oplus 2^{-52} = y^{(52)}+2^{-52}\in[1\lldot 2]$ が成り立つ。

$z = (r^{(51)}-y^{(51)}) - (2^{-1}\cdot(y^{(52)} + 2^{-52}))$ として $z \equiv 0 \pmod{2^{-53}}$ が成り立つので、$z\ge 2^{-55} \iff z\gt 0$ となる。あとは Lemma 7 より従う。$\qed$

以上により、$\roundcirc{y^{(53)}}$ が求められました。めでたしめでたし。

手順

手順として整理すると次の通りです。

便宜上、Iverson bracket $[\bullet]$ を用います。すなわち、条件 $x$ が成り立つとき $[x]=1$、そうでないとき $[x]=0$ です。

• 入力: $x\in[1\lldot 4)$
• 出力: $\roundcirc{\sqrt x}$
• $y\gets 1$ で初期化する。
• $r\gets x\ominus 1$ で初期化する。
• $\Delta y \gets 0.5$ で初期化する。
• 各 $i\gets\angled{1, 2, \dots, 51}$ について下記を行う。
  • $r\xgets{\otimes}2$ で更新する。
  • $z \gets (2\otimes y)\oplus\Delta y$ で定義する*2。
  • $r\ge z$ であれば下記を行う。
    • $r\xgets{\ominus}z$ で更新する。
    • $y\xgets{\oplus}\Delta y$ で更新する。
  • $\Delta\xgets{\otimes}0.5$ で更新する。
• $(y^{(51)}, r^{(51)}) \gets (y, r)$ で定義する。
• $y_{52} \gets [r^{(51)} \gt y^{(51)}]$ で定義する。
• $y_{52} = 1$ であれば下記を行う。
  • $y\xgets{\oplus}\Delta y$ で更新する。
• $y_{53} \gets [y_{52} = 0][2\otimes r^{(51)}\ge y^{(51)}] + [y_{52} = 1][r^{(51)}\ominus y^{(51)}\gt 0.5\otimes(y\oplus \Delta y)]$ で定義する。
• $y_{53} = 1$ であれば下記を行う。
  • $y\xgets{\oplus} \Delta y$ で更新する。
• $y$ を出力する。

実装

fn sqrt(x: f64) -> f64 {
    if x < 0.0 {
        return f64::NAN;
    }
    if x == 0.0 || x.is_infinite() || x.is_nan() {
        return x;
    }

    let (x_m, x_e) = match frexp(x) {
        (m, e) if e % 2 == 0 => (4.0 * m, e - 2),
        (m, e) => (2.0 * m, e - 1),
    };
    assert_eq!(x_e % 2, 0);
    assert!(1.0 <= x_m && x_m < 4.0);

    let mut y = 1.0;
    let mut r = x_m - 1.0;
    let mut two_pmn = 0.5;
    for _ in 1..=51 {
        r *= 2.0;
        let tmp = 2.0 * y + two_pmn;
        if r >= tmp {
            r -= tmp;
            y += two_pmn;
        }
        two_pmn *= 0.5;
    }

    let (y51, r51) = (y, r);

    let y_52 = if r51 > y51 { 1 } else { 0 };
    if y_52 == 1 {
        y += two_pmn;
    }

    let y_53 = if y_52 == 0 {
        if 2.0 * r51 >= y51 { 1 } else { 0 }
    } else {
        if r51 - y51 > 0.5 * (y + two_pmn) { 1 } else { 0 }
    };
    if y_53 == 1 {
        y += two_pmn;
    }
    ldexp(y, x_e / 2)
}

おまけ

各整数 $0\le n\le 52$ に対して $y^{(n)}\in[1\lldot 2)$, $r^{(n)}\in[0\lldot 4)$, $r^{(n)}\equiv 0\pmod{2^{-52}}$ が成り立つことから、それぞれ $2^{52}$ 倍したものを考えて u64 で計算することも可能です。実際、LLVM の実装ではそのようになっています。 前回の二分探索の方針では u128 が必要な上に乗算が必要でしたが、今回の方針では u64 の加減算とシフト演算のみで実装可能ですね。

所感

なんでこんなつらい思いをしながら f64 で計算可能なことを示したがってるんだ？ 素直に u64 を使いたい気持ちでいっぱい

— えびちゃん🍑🍝🦃 (@rsk0315_h4x) 2025å¹´3月23æ—¥

示せるとうれしいという気持ちはあります。

オーバーフローしたときに上の桁を捨てるのが u64 で、下の桁を捨てる（+ どのくらいオーバーフローしたかを別途持っておく）のが f64 だという感覚が身についてきました*3。 そういう意味で、「誤差を出さない」と「オーバーフローさせない」がある種の似た見方だと思うようになってきました。

Elementary Functions では、shift-and-add algorithm の使用例として log や exp を挙げており、こちらも後々紹介することになるかな？と思っています。 1 ビットずつ決めていくため仮数部の長さに対して線形時間はかかってしまいますが、値が $0$ または $1$ であることを利用した式変形が使えがちで面白いですね。 今回の例で言うところの $$ y_n\cdot(2y^{(n-1)}+y_n\cdot 2^{-n}) = y_n\cdot(2y^{(n-1)}+2^{-n}) $$ のようなものです。

気分としては、次は sin, cos をしようかなと思っていますが、どうなるかはわかりません。 魔法の多項式や Newton 法とも仲よくなりたい気持ちはありますからね。

おわり

おわりです。

*1:LLVM のコメント上では最後の項が $2^{-n-1}$ になっていますが、その後の実装を見ても $2^{-n}$ 相当になっているように見えますし、これで正しいと思っています。

*2:immutable 前提の変数の導入をしたいとき、「初期化」は違うかもと思って用語を分けました。「束縛」？

*3:浮動小数点数で言うところのオーバーフローは $\pm\infty$ になることだと思いますが、ここでは仮数部に対して四則演算をしたときの桁あふれくらいの意味で解釈してください。

画像のてすと

ここにテスト用の画像を貼る

ここにテスト用の画像を貼ります。

correct rounding への道 (1) ç´ æœ´ãª sqrt

いろいろな数学関数たちの correct rounding な実装をしていこうという遊びです。 丸め方向については一旦は tiesToEven のみを考えています。

$\sin$ や $\log$ のような超越関数 (transcendental functions) は大変ではありますが、さまざまな典型テクがあるので、そのうち紹介することになると思います。 今回は、特に洗練されていない方法を用いて $\sqrt{x}$ の correctly-rounded な値を得ましょうという回です。

前提

ç«‹å ´

ここでは correct rounding であることを一番に重要視します。すなわち、真の値が $f(x)$ のとき、それを正確に丸めた値が欲しいです。 いくらか調べたところ、この分野では $\roundcirc{f(x)}$ と書くことが多そうなので、ここでもそうします*1。今までの記事で $\roundp{f(x)}$ のように書いていたものと同じです。

初めての人からすると「正確に丸める」というのが慣れない概念かもしれないので、一応具体例を出しておきます。$\roundcirc{\sqrt2}$ を考えます。 $$ \sqrt{2} = {\small 1.}{\footnotesize 414213562373095}{\scriptsize 048801688724209}{\tiny 698078569671875}{\tiny 376948073176679{\dots}} $$ で、binary64 で表せる数のうち、$\sqrt{2}$ 以下で最大のものと $\sqrt{2}$ 以上で最小のものは、それぞれ下記の通りです。 $$ \begin{aligned} \texttt{16A09E667F3BCC}_{(16)}\times 2^{-52} &= {\small 1.}{\footnotesize 414213562373094}{\scriptsize 923430016933707}{\tiny 520365715026855}{\tiny 46875}, \\ \texttt{16A09E667F3BCD}_{(16)}\times 2^{-52} &= {\small 1.}{\footnotesize 414213562373095}{\scriptsize 145474621858738}{\tiny 828450441360473}{\tiny 6328125}. \end{aligned} $$

前者との差は $1.253\times 10^{-16}$ 程度、後者との差は $9.667\times 10^{-17}$ 程度であり、後者の方が近いため、後者に丸められることになります*2。 ここで前者に丸めてしまうものは correct とは認められません。なお、負方向または正方向で最も近いもののどちらか（方向は引数によって変わってよい）に丸めるものは faithful rounding と呼ばれます。

note: $\sqrt{2}$ や $\tfrac13$ など、（二進法で）無限小数となるものばかりが意識されがちですが、$(1+2^{-52})^2$ のような有限小数でも丸め誤差は生じます。

応用先によっては faithful rounding でよいでしょうし、あるいはもっとラフに「相対誤差 $10^{-14}$」のようなので満足できる状況も多々あるでしょうが、ここではそうではないということです。誘惑に負けずにいきましょう。

また、正当性の証明を重要視します。これは浮動小数点型に限った話ではないですが、未証明のアルゴリズムは使う気がしないですからね。 人間が数式で行うのが主になると思いますが、Gappa や Coq などの証明支援システムを使うことも視野に入れています。

速度も速いに越したことはないですが、二の次です。たとえ「faithfully-rounded な結果は 0.01 ms で得られるが、correctly-rounded な結果は 30 s かかる」としても、我々は correct rounding に固執します。 正当性を保ったままで高速化をがんばればよいのであって、correct rounding を捨てるのは論外です。 correct rounding の狂信者になったと思ってください。

åž‹

浮動小数点型として binary64 (double, f64) を考えます。丸めモードとしては tiesToEven だけを一旦想定します。

プリミティブ演算

• 定数リテラル
  • NaN ã‚„ $\pm\infty$ を含む定数の生成
• ビット表現同士の変換 (transmute)
  • f64 から u64 へのビット表現を保った変換
  • u64 から f64 へのビット表現を保った変換
• 128-bit 以下の整数型（符号なし、符号つき）の各演算
  • 四則演算 (+, -, *, /, %)
    • / および % は、商を $0$ 方向に丸める切り捨て除算とする*3
  • ビット演算 (&, |, ^, !, <<, >>)
  • 比較演算 (==, !=, <, >, <=, >=)
• 浮動小数点型の各演算
  • 四則演算 (+, -, *, /)
    • 数式中では $\oplus$, $\ominus$, $\otimes$, $\oslash$ と表記する
  • fused multiply-add (FMA)
    • $\roundcirc{x\times y + z}$ のこと
  • 比較演算 (==, !=, <, >, <=, >=)
  • NaN の判定

浮動小数点型の各演算は correct rounding とします。後々になってから整数型などでエミュレートする方法を書くかもしれませんが、一旦は与えられているものとします。

しばしば暗黙に行いますが、$2^{-1022}$ や $2^{54}$ のような定数は、指数関数のようなものを意図しているのではなく、あくまで定数リテラルに相当するものとして使っています。

上記の演算のみで計算できるアルゴリズムを構成し、その正当性の証明の際に必要であれば無限精度の実数の演算も用います。数式中の計算は、丸め $\roundcirc{\bullet}$ を明示しない限りは無限精度で（というかただの実数の演算として）行います*4。

note: $x \oplus y$ のような浮動小数点数用の演算子は、オペランドの $x$ や $y$ が浮動小数点数として表せる値である場合のみ使う想定ですが、そういう前提で $\roundcirc{x+y}$ の略記だと思って差し支えないです。

本題

やりましょう。

frexp

まずは frexp と呼ばれる関数を用意しておきましょう。 与えられた浮動小数点型の値を、fraction と exponent に分ける関数です。 具体的には、非零な有限値 $x$ が与えられたとき、下記を満たすような組 $(x_m, x_e)$ を返します。$x_m$ は浮動小数点型で、$x_e$ は符号つき整数型です。

• $x_m\times 2^{x_e} = x$, and
• $|x_m|\in[0.5\lldot 1)$.

$x$ がゼロ ($0_{\pm}$) または NaN、$\pm\infty$ のときは $(x, 0)$ を返すことにしておきます。

ゼロの符号について

正負のゼロ +0.0 と -0.0 をそれぞれ $0_+$, $0_-$ と書いています。$\roundcirc{0} = 0_+$ ã‚„ $-0_+ = 0_-$ などが成り立ちます。 +0.0 == -0.0 ですが、記号の比較という意味で $0_+ \ne 0_-$ としておきます。

下記の記事でも多少触れています。

rsk0315.hatenablog.com

誤差評価の過程ではゼロの符号は問題にならないので、数式中では単に $0$ と書いておくことにします。$\eod$

signaling NaN や quiet NaN の区別などに関しては、一旦気にしないことにしておきます*5。

ビット表現を得てから指数部を調整します。非正規数に注意しましょう。glibc の実装では $2^{54}\cdot x$ のケースに帰着させていますね*6。

実装

const TWO_P54: f64 = 18014398509481984.0;

fn frexp(x: f64) -> (f64, i32) {
    let mut ix = x.to_bits();
    let mut ex = (ix >> 52 & 0x7FF) as i32;

    if ex != 0x7FF && x != 0.0 {
        // Not zero and finite.
        let mut e = ex - 1022;
        if ex == 0 {
            // subnormal.
            let x = x * TWO_P54;
            ix = x.to_bits();
            ex = 0x7FF & (ix >> 52) as i32;
            e = ex - 1022 - 54;
        }
        ix = (ix & 0x800FFFFFFFFFFFFF_u64) | 0x3FE0000000000000_u64;
        (f64::from_bits(ix), e)
    } else {
        (x, 0)
    }
}

ldexp

続いて frexp の逆操作 (load exponent) を用意しておきましょう。$(x_m, x_e)$ に対して $\roundcirc{x_m\times 2^{x_e}}$ を返します。 $x_e$ の範囲次第では overflow ã‚„ underflow が起きるので $\circ$ を明示していますが、それ以外のケースでは誤差は生じません。 引数は $x_m\in[0.5\lldot 1)$ である必要はありません。

実装

const TWO_P54: f64 = 18014398509481984.0;
const TWO_PM54: f64 = 5.5511151231257827021181583404541015625e-17;

fn ldexp(mut x: f64, e: i32) -> f64 {
    if !x.is_finite() || x == 0.0 {
        return x;
    }
    let mut ix = x.to_bits();
    let mut k = (ix >> 52 & 0x7FF) as i32;
    if k == 0 {
        // subnormal.
        x *= TWO_P54;
        ix = x.to_bits();
        k = (ix >> 52 & 0x7FF) as i32 - 54;
    }
    if e < -50000 {
        return 0.0_f64.copysign(x);
    }
    if e > 50000 || k + e > 0x7FE {
        return f64::INFINITY.copysign(x);
    }
    k += e;
    if k > 0 {
        let bits = (ix & 0x800FFFFFFFFFFFFF_u64) | (k as u64) << 52;
        return f64::from_bits(bits);
    }
    if k <= -54 {
        return 0.0_f64.copysign(x);
    }
    k += 54;
    let bits = (ix & 0x800FFFFFFFFFFFFF_u64) | (k as u64) << 52;
    f64::from_bits(bits) * TWO_PM54
}

sqrt

さて、いよいよ本題です。$\roundcirc{\sqrt{x}}$ を求めます。

まず correct rounding な実装の典型的なパターンとして、次のような流れがあります。

1. range reduction
   • 入力を一般のケースから特殊なケースに帰着させる
2. approximation
   • 所望の関数に応じた手法を用い、特殊なケースについての答えを求める
3. reconstruction
   • 2. の結果を用い、元々の入力に対する答えを求める

今回もそれに従って行います。負の数や NaN などのケースについては予め処理しておくことにします。 また、$\sqrt{0_{\pm}\vphantom{0^0}} = 0_{\pm}$（複号同順）と定義されているので、それも先に処理しておきます。

range reduction

$\sqrt{2^{2 x_e}\cdot x_m} = 2^{x_e}\cdot\sqrt{x_m\vphantom{2^2}}$ であることから、$x_m\in[1\lldot 4)$ のケースについて考えます。 $\texttt{frexp}(x) = (x_m', x_e')$ に対して下記を行うことで帰着できます。 $$ (x_m, x_e) = \begin{cases} (4x_m', x_e'-2), & \text{if }x_e' \equiv 0\pmod 2; \\ (2x_m', x_e'-1), & \text{if }x_e' \equiv 1\pmod 2. \end{cases} $$

approximation

最初に、$y\in[1\lldot 2)$ であって、$y\le \sqrt{x_m\vphantom{2^2}}\lt y+2^{-52}$ となるものを求めます。 これは、二分探索の要領で、以下の手続きによって可能です。

• 入力: $x_m\in[1\lldot 4)$
• $y \gets 1$ で初期化する。
• $\Delta y \gets 0.5$ で初期化する。
• 各 $i \gets \angled{1, 2, \dots, 52}$ について下記を行う。
  • $\roundcirc{(y\oplus\Delta y)\times(y\oplus\Delta y)+(-x_m)}\le 0$ であれば下記を行う。
    • $y\xgets{\oplus}\Delta y$ で更新する。
  • $\Delta y \xgets{\otimes} 0.5$ で更新する。
• $y$ を出力する

note: 各ループの先頭において、$\Delta y = 0.5^i$ が成り立ちます。

Claim 1: 任意の浮動小数点数 $x\in[1\lldot 2)$, $y\in[1\lldot 2]$, $z\in[1\lldot 4)$ に対し、$$x\times y\le z \iff \roundcirc{x\times y+(-z)}\le 0$$が成り立つ。

Proof

任意の実数 $w$ に対し、$\roundcirc{w}\le 0 \iff w\le 2^{-1075}$ である。 よって、 $$ x\times y\le z \iff x\times y-z \le 2^{-1075} $$ を示す。

($\implies$): 明らか。

($\impliedby$): 対偶 $x\times y\gt z \implies x\times y-z \gt 2^{-1075}$ を示す。

ある整数 $m_x\in[2^{52}\lldot 2^{54})$, $m_y\in[2^{52}\lldot 2^{53}]$, $m_z\in[2^{52}\lldot 2^{54})$ を用いて $x = m_x\times 2^{-52}$, $y = m_y\times 2^{-52}$, $z = m_z\times 2^{-52}$ と表せる。 $$ \begin{aligned} x\times y - z &= (m_x\times 2^{-52})\times (m_y\times 2^{-52})-(m_z\times 2^{-52}) \\ &= (m_x\cdot m_y - m_z\cdot 2^{52})\times 2^{-104} \end{aligned} $$ より $x\times y-z$ は $2^{-104}$ の倍数であるから、 $$x\times y\gt z\implies x\times y-z\ge 2^{-104}\gt 2^{-1075}. \quad\qed$$

すなわち、$\roundcirc{(y\oplus\Delta y)\times(y\oplus\Delta y)+(-x_m)}\le 0$ の部分は $(y\oplus\Delta y)^2\le x_m$ と同値です*7。 $y\oplus\Delta y = y+\Delta y\in[1\lldot 2)$ は常に成り立ちます。 また、$y$ や $\Delta y$ の更新に際して生じる誤差は $0$ です。

次に、丸めの境界に真の値が現れることはないことを示します。

Claim 2: 任意の $x_m$ に対して $y+2^{-53} \ne \sqrt{x_m\vphantom{2^2}}$ が成り立つ。

Proof

背理法による。

上記の手続きの結果 $y+2^{-53} = \sqrt{x_m\vphantom{2^2}}$ なるような浮動小数点数の組 $(x_m, y)$ が存在したとする。 このとき、$x_m = (y+2^{-53})^2 = y^2+2^{-52}\cdot y+2^{-106}$ となる。

$y\in[1\lldot 2)$ より $y$ は $2^{-52}$ の倍数であり、 $$ y^2+2^{-52}\cdot y + 2^{-106} \equiv 2^{-106} \not\equiv 0 \pmod{2^{-104}} $$ となる。一方、$x_m\in[1\lldot 4)$ であったから、$x_m$ は $2^{-52}$ の倍数であり、$x_m\equiv 0\pmod{2^{-104}}$ となる。$\qed$

よって、$y+2^{-53}\lt \sqrt{x_m\vphantom{2^2}}$ であれば $y\oplus 2^{-52}$ が、そうでなければ $y$ が approximation step の答えとなります。

非負性は明らかなので、$(y+2^{-53})^2\lt x_m$ と同値です*8。 すなわち、$y^2+2^{-52}\cdot y+2^{-106}-x_m\lt 0$ を判定できればよいです。

Claim 3: 任意の浮動小数点数 $x\in[1\lldot 4)$ および $y\in[1\lldot 2)$ に対し、$$y^2+2^{-52}\cdot y-x\lt -2^{-106} \iff y^2+2^{-52}\cdot y-x\lt 0$$が成り立つ。

Proof

($\implies$): 明らか。

($\impliedby$): 対偶 $y^2+2^{-52}\cdot y-x\ge -2^{-106} \implies y^2+2^{-52}\cdot y-x\ge 0$ を示す。

ある整数 $m_x\in[2^{52}\lldot 2^{54})$, $m_y\in[2^{52}\lldot 2^{53})$ を用いて $x = m_x\times 2^{-52}$, $y = m_y\times 2^{-52}$ と表せる。 $$ \begin{aligned} y^2+2^{-52}\cdot y-x &= (m_y\times 2^{-52})^2 + 2^{-52}\cdot(m_y\times 2^{-52}) - (m_x\times 2^{-52}) \\ &= (m_y^2+m_y-m_x\cdot 2^{52})\times 2^{-104} \end{aligned} $$ より $y^2+2^{-52}\cdot y-x$ は $2^{-104}$ の倍数であるから、 $$y^2+2^{-52}\cdot y-x\ge -2^{-106} \implies y^2+2^{-52}\cdot y-x\ge 0.\quad\qed$$

Claim 3 より $y^2+2^{-52}\cdot y-x_m\lt 0$ を判定すればよいことになりました。 すなわち、$y\times(y+2^{-52})\lt x_m$ です。これは Claim 1 から FMA で計算可能です。

reconstruction

さて、ここまでで $\roundcirc{\sqrt{x_m\vphantom{2^2}}} = y$ が得られました。$\texttt{ldexp}(y, \tfrac{x_e}2)$ が最終的な答えです。

実装

const TWO_PM52: f64 = 2.220446049250313080847263336181640625e-16;

fn sqrt(x: f64) -> f64 {
    if x < 0.0 {
        return f64::NAN;
    }
    if x == 0.0 || x.is_infinite() || x.is_nan() {
        return x;
    }

    let (x_m, x_e) = match frexp(x) {
        (m, e) if e % 2 == 0 => (4.0 * m, e - 2),
        (m, e) => (2.0 * m, e - 1),
    };
    assert_eq!(x_e % 2, 0);
    assert!(1.0 <= x_m && x_m < 4.0);

    let mut y = 1.0_f64;
    let mut dy = 0.5;
    for _ in 0..52 {
        if (y + dy).mul_add(y + dy, -x_m) <= 0.0 {
            y += dy;
        }
        dy *= 0.5;
    }

    if y.mul_add(y + TWO_PM52, -x_m) < 0.0 {
        y += TWO_PM52;
    }

    let y = ldexp(y, x_e / 2);
    assert_eq!(y, x.sqrt());
    y
}

次回予告

とりあえず $\roundcirc{\sqrt{x}}$ が計算可能であることは示しましたが、ビット長 $p$ に対して $\Theta(p)$ 時間というのはちょっとね... というのが正直なところです。 たとえば $\Theta(\log(p))$ 時間くらいになってくれたらうれしいですよね。

ということで一旦実験しました。

まず $3$ 次くらいの魔法の多項式でざっくり近似値（誤差 $10^{-4}$ くらい）を求めて、Newton 法で $2$ 回くらい反復させると、1 ULP くらいの誤差に収まってくれそう？ 1 ULP（というか定数 ULP）に収まることが示せるのであれば、今回やったように不等式評価して補正すればいいですからね。

魔法の多項式は次のような感じです。 $$ \begin{aligned} f_0(x) &= \left(\begin{matrix} \roundcirc{0.371351660146978} \\ \roundcirc{0.784942635287931} \\ \roundcirc{-0.180689144911217} \\ \roundcirc{0.0244769088312593} \end{matrix}\right)^{\top}\cdot\left(\begin{matrix} x^0 \\ x^1 \\ x^2 \\ x^3 \end{matrix}\right), \\ f_1(x) &= \left(\begin{matrix} \roundcirc{0.525170554189621} \\ \roundcirc{0.555038260254535} \\ \roundcirc{-0.0638832598267602} \\ \roundcirc{0.00432694705426709} \end{matrix}\right)^{\top}\cdot\left(\begin{matrix} x^0 \\ x^1 \\ x^2 \\ x^3 \end{matrix}\right). \\ \end{aligned} $$ $f_0$ は $[1\lldot 2)$ 用、$f_1$ は $[2\lldot 4)$ 用です。 反復は $y \xgets{\oplus} (x_m\oslash y)$ と $y\xgets{\otimes} 0.5$ です。 下記は $\sqrt{2}$ を求めようとしているところのイメージです。

fn main() {
    let poly = |x: f64| {
        let mut y = 0.00432694705426709_f64;
        y = y.mul_add(x, -0.0638832598267602);
        y = y.mul_add(x, 0.555038260254535);
        y.mul_add(x, 0.525170554189621)
    };
    let iter = |x: f64, y: f64| (y + x / y) * 0.5;
    
    let x = 2.0;
    let y = poly(x);    // 0x1.6A1181648E5E8p0, 1.4143296118257869
    let y = iter(x, y); // 0x1.6A09E67C6699Fp0, 1.414213567134176
    let y = iter(x, y); // 0x1.6A09E667F3BCCp0, 1.414213562373095
    let y = iter(x, y); // 0x1.6A09E667F3BCCp0, 1.414213562373095
}

初期値を魔法で決めずに $1$ とかから始めても $6$ 回程度で収束してくれていそうでした。

disclaimer: えびちゃんは気分屋さんなので、この通りに進まないかもしれません。

とりあえず hypot, cbrt, log, exp, sin, cos あたりをやるまでに飽きないといいですね。順番は未定です。

おまけ

二分探索するのであれば、f64 を使わない方法もあって、たとえば u128 などを使いながら整数として二分探索してもよいですよね。

整数 $m_x\in[2^{52}\lldot 2^{54})$, $m_y\in[2^{52}\lldot 2^{53})$ であって $$ (m_y-\tfrac12)\times 2^{-52} \lt \sqrt{m_x\times 2^{-52}} $$ なる $y$、すなわち $$ \begin{aligned} m_y-\tfrac12 &\lt 2^{52-26}\cdot\sqrt{m_x\mathstrut} \\ \Floor{(2m_y-1)^2\times 2^{-54}} &\lt m_x \end{aligned} $$ を満たすような $m_y$ の最大値を求めればよいです。

fn sqrt(x: f64) -> f64 {
    // ... Same as the preceding one.
    assert_eq!(x_e % 2, 0);
    assert!(1.0 <= x_m && x_m < 4.0);

    let m_x = {
        let mant = (1 << 52 | x_m.to_bits() & !(!0 << 52)) as u128;
        if x_m >= 2.0 { 2 * mant } else { mant }
    };
    let m_y = {
        let mut lo = 1_u128 << 52;
        let mut hi = 2 * lo;
        while hi - lo > 1 {
            let mid = lo + (hi - lo) / 2;
            let too_lo = (2 * mid - 1) * (2 * mid - 1) >> 54 < m_x;
            *(if too_lo { &mut lo } else { &mut hi }) = mid;
        }
        lo as u64
    };
    let y = f64::from_bits(0x3FF << 52 | (m_y & !(!0 << 52)));
    ldexp(y, x_e / 2)
}

また、下記の記事の Corollary 10 と u128::isqrt を用いることで、もっと簡単に求めることができます。

rsk0315.hatenablog.com

fn sqrt(x: f64) -> f64 {
    // ... Same as the preceding one.

    let m_x_p52 = m_x << 52;
    let m_y = m_x_p52.isqrt();
    let m_y = (if m_y * (m_y + 1) < m_x_p52 { m_y + 1 } else { m_y }) as u64;

    let y = f64::from_bits(0x3FF << 52 | (m_y & !(!0 << 52)));
    ldexp(y, x_e / 2)
}

こうして、u128 で $\floor{\sqrt{x_m\mathstrut}\times 2^{104}}$ を求めることができれば、それを使った定数回の演算で f64 での $\roundcirc{\sqrt{x_m\mathstrut}}$ も求められることがわかってしまいました。 Rust の isqrt では、Karatsuba square root algorithm を用いているよう ([src]) です。競プロ er にはあまり馴染みのないアルゴリズムなのではないでしょうか*9。

所感

う〜〜む、結局「整数型でなんとかなるじゃん」という例が出てしまい、「これでよかったのか...？」という気持ちです。 もちろん、こんなのはおそらくsqrt のような簡単な algebraic な関数だからできることで、今後の transcendental な関数ではこうはいかないと思っています。 次回は、それに備えて諸々の典型テクを導入するはずです。

おわり

おわりです。

*1:おそらく round-to-the-nearest から $\mathrm{RN}(f(x) )$ と書いている文献もあります。また、任意の丸めモードを考えるときは $\operatorname{\diamond}{(f(x) )}$ と書かれることが多そうでした。

*2:$\sqrt2$ を正確に表せないこと自体には特に関心がなく、$\roundcirc{\sqrt2}$ を求めることに興味があります。

*3:$-\infty$ 方向丸めの方がやりやすいが...

*4:$\roundcirc{\bullet}$ の $\circ$ は丸めの関数で、$\bullet$ は（$\sqrt{\bullet\mathstrut}$ のように使う）プレースホルダです。形が似ていてわかりにくいですね。

*5:glibc の実装が規格に反しているような気がしつつ、あまり自信がないですね。

*6:$2^{52}\cdot x$ でいい気がするのですが、意図はよくわかりません。

*7:$x^2$ は通常の $x\times x$ のことであり、$x\otimes x$ のことではありません。

*8:$y\oplus 2^{-53} \ne y+2^{-53}$ なので、この時点では Claim 1 は使えません。

*9:二分探索で済むなら他のことを学ぶ理由がないと思っている人々が多数派だから（偏見）。

浮動小数点型に関するポエム20250320

えびちゃんです。お気持ち表明記事です。最近「お気持ち表明」というフレーズをあまり聞かなくなった気がしますね。

各章はそれぞれ基本的に独立です。特に結論めいたものはありません。

出会い

えびちゃんが競プロを始めたのは 2016 年頃で、それ以来なんやかんやで浮動小数点型には「ブラックボックス」感を抱いていた記憶があります。 競プロ界隈は別に浮動小数点型に詳しい界隈ではないので、（#define EPS 1e-6 などに代表されるような）怪しい民間療法のようなものが当時から信仰されていました。

「できる限り整数型で処理できるように式変形・考察しましょう」という定石は全く正しいと思っていますが、浮動小数点型への忌避意識のようなものは強まるばかりです。 そうして避けたまま、いざ避けられなくなったときに「基礎が身についていないため、典型的なミスをする」というのがありがちなパターンな気がします。

一生そうした感覚で生き続けることに耐えられなくなったため、2023 年頃にちゃんと勉強を始めました。 もっとも、浮動小数点型に限らず、そうしたブラックボックスめいたものの内部を知りたいがちな性格ではあったので、ようやくかという感じではあったかもしれません。

整数型との対比

これはぱっと思いついた例ですが、JPEG は「よくわからんがぐちゃぐちゃになっちゃうもの」という印象があり、PNG は「きれいなもの」という印象がありますよね。浮動小数点型と整数型の対比もそれと似ている気がします。

整数型はきっちりとした計算を行うことができるもので、たとえば 0x29AEDBEFB219581E18E1217B2715654F ã‚’ 0x89FC309DF824DC7C で割ったあまり (%) などもやろうと思えば手計算で模倣することができると思われていそうです。たいへん面倒なのでやりたがる人はいないと思いますが、やり方自体がわからないということはないと思います。仮にお金がもらえるならやりますよね。

一方で、浮動小数点型においては、「コンピュータのお気持ち次第で誤差が出るような予測不能な変なもの」と思っている人がそれなりにいるのではないでしょうか。 たとえば 0x1.26ECE273B8E6Ap-11 と 0x1.EC8689D97BFB8p+3 の積 (*) を手計算で模倣することはできますか？

このへんは人によるかもで、「あくまで float -> float -> float のような関数の挙動が実質的には予測不能なだけで、値を固定すれば手計算（あるいは整数型の演算に帰着させるなどして）でも模倣可能」くらいに思っている人も多いのかもしれません？ 多くの人がどう思っているのかもうわかりません。

とはいえ、整数型での演算でもブラックボックスなものはいくらでもあって、そうしたものについては「ほんとにどんなときでもこわれないの？」という気持ちは生じそうなので、浮動小数点型ばかりが除け者というわけではないのかもしれません。AtCoder Library の math にあるようなものたちを想像して話しています。

rsk0315.hatenablog.com

固定小数点型や十進浮動小数点型など

浮動小数点型はカスだが固定小数点型はカスではないとか、浮動小数点型の NaN はカスだとか、あれこれ理由をつけて（二進）浮動小数点型を貶めようとする言説はしばしば見かけます。

下記などが参考になります。

qiita.com

浮動小数点型はそんなに目の敵にされるべき発明ですか？ むしろ概念自体は直感的まであるような気がします。 批判の要旨が「直感的・素朴すぎて、もっと洗練されるべきだ」のような感じなら納得しなくもないですが、もっと素朴で扱いにくいものが代替案として出されがちなので微妙な気持ちになります。

多倍長整数が定数時間で計算できると思って使っていそうな人々もいますし、自分の関心のない部分は「全部勝手にうまくいってくれる」と（無意識に？）思って代替案を出しているんだろうと思うと、理解はできます。

誤差

よくある 0.1 + 0.2 != 0.3 のような話題はもう見飽きましたが、「「$0.1$ を浮動小数点型で表せるように丸めたやつ」と「$0.2$ を浮動小数点型で表せるように丸めたやつ」の和を浮動小数点型で表せるように丸めたやつ」が「$0.3$ を浮動小数点型で表せるように丸めたやつ」と異なるという事実自体は、別にそこまで驚くべきことでもないのではないか？という気持ちはあります。

そもそも（オペランド自体の丸めもそうですが）+ 自体が単なる $+$ とは異なるものなので、「浮動小数点型では $0.1 + 0.2 = 0.3$ は成り立たない」のような言い方はちょっとうれしくないなぁと思ってしまいます。

よくある REPL などで「$0.1$ を浮動小数点型で表せるように丸めたやつ」が「0.1」と表示されることは非常に confusing で、変な誤解を生む要因の一つなんじゃないかなぁという思いはあります。下記みたいなのが人類にとって幸福かと言われると、それはそれで悩むところではあります。

>>> 0.1
0.1000000000000000055511151231257827021181583404541015625

これだと結局は「わけのわからないぐちゃぐちゃな値を $0.1$ に足したもの」だと認識するしかなく、「浮動小数点型は変なやつ」として避けられそうな気がします。 初心者に「実数を模倣できるやつ」だと誤解されて乱用されるよりはマシなのかもしれない気もします。

えびちゃんは、上記の値が $\tfrac1{10}(1+2^{-54})$ と正確に等しい値だと知ってから、それだけでお友だち度が少し増したような気持ちになりました。

比較

よく言われている「セオリー」の一つに、浮動小数点型の比較は x == y で行うのではなく、1e-15 などの十分小さい値を使って (x - y).abs() < 1e-15 のようにして行いましょうというものがあります。

やっぱりこれはどうにも怪しいという感想になってしまいます。1e-15 のような値が適切であるということはどうやって証明したのですか？という話です。 それが「なんとなく」なのであれば、怪しいのは浮動小数点型ではなく実装者です*1。

この小さい値の決め方（あるいはこうした比較の仕方が可能かどうか自体も）は状況次第であり、都度証明する必要がある*2ものです。 そもそも一般に、アルゴリズムの正当性というものに興味がない人ももしかして多いですか？ 困ったものですね。

仮数部のビット長などの意味での精度 (precision) と、その型を用いて何らかの式を計算した際の精度 (accuracy) を区別せずに認識していて、「precision が 53 bits であれば、（ある程度複雑な式の）計算結果も $2^{-53} \approx 10^{-16}$ 程度くらいになるだろう」と思っている人もいるのではないでしょうか。さすがにそんなことはないですか？

数式での証明

IEEE 754 の規格によって、浮動小数点型の計算は（IEEE 754 に準拠したものであれば）環境によらずに再現性のあるものになっています。 この規格の登場前はプロセッサなどごとにさまざまな挙動をしていてカオスだったと聞いています。

再現性があるということは数式などによる証明も有効で、これは非常にうれしいことです。 にもかかわらず、やはりわざわざ「float を使うたびに数式（など）で証明を書く」という営みをしている人はそう多くはないような気がします。 なにをどう証明すればいいかもわからないし、そもそも証明できる（あるいは証明するべき）ことだと思っていない人が多数派なのかな？と思っています。

最近では AI によるコード生成が流行っていますが、（Gappa なり Coq なりを駆使しつつ）証明もやってくれたらうれしいものでしょうか。できるのかな？ できないとはあまり思ってないですけど。

数学関数たちについて

数学関数たちの実装は必ずしも correctly-rounded になっていないことも多く、嫌な気持ちになります。LLVM は correct rounding な実装を 目指していそう なので、うれしい気持ちになります。がんばえーっっ

実務？

長めの保守を前提とする実務のコードで使いたいか？となってくると、これは話が別です。

リファクタの名目で括弧のつけ方を勝手に変えられるだけでめちゃくちゃになったり、そもそも証明を理解してメンテできる人の教育コストはたくさんだったり、単体テストでなにかを保証するのが大変だったり、本番環境でバグが出たときの再現に困ったりなど、ぱっと思いつくだけでも嫌な要素ばかりです。う〜〜ん参った。

でも現実はたぶんもっと酷くて、普通に桁落ちでめちゃくちゃになりうるようなコードが世の中の本番環境にたくさんデプロイされてたりするんじゃないですか？とも思ったりで、嫌〜〜という気持ちになります。

というわけなので、えびちゃんは浮動小数点型が好きではあるものの、「世の中の人がもっと使ってくれ〜〜」と思っているわけでもないです。 「世の中の人が浮動小数点型関連の証明をできるようになってくれ〜〜」とも別に思ってはいないですが、できるようになった世界も見てみたくはありますね。

浮動小数点型を使ったコードを書いているとレビュー時に反例を挙げて reject してくるえびちゃん概念、怖くて嫌ですね*3。あるいは「証明教えてください」とか言われても嫌だな。

所感

今までのも全部所感だろというのは正しいですね。メタ所感？

ネット上で見かけた怪しい言説や、浅い理解の解説ブログなどの印象に引っ張られすぎて、「大半の人類はまともに理解してないんじゃない？」と思っているだけかもしれません。書きながらそういう自覚は芽生えてきました。錯覚だったらうれしいです。

なにかしらを見るたびにもやもやするのに飽きたので、何度も思ったようなことをばーっと書き殴ったという感じです。

数値計算のエキスパートみたいな人自体はごろごろいるはずですから、えびちゃんもそういう人になりたいですね。 過去何度か挫折していますが、最近は数学関数の correct rounding な実装について調べています。 correct rounding な数学関数のライブラリに関しては、すでにメンテが終了したプロジェクトが複数あってかなしい気持ちです。 サーベイ論文 が 2025 年 2 月に更新されていたりして、まだまだ動きのある分野なのかな？と思ったりしています。 それなりに古くから（たとえば競プロ界隈ではあまり有名でないような）さまざまなことが研究されていそうなので、がんばらなきゃというところです。

「競プロ界隈で有名」ってなんだよという指摘はありますね。競プロ界隈で有名なものはセグ木と union-find しかないかもしれません。

おわり

明日からがんばります。

*1:過激発言すぎる。

*2:毎回厳密な証明をしましょうというよりは、「こういう根拠があるからいけるよね」くらいの考慮すらしない態度はだめだよね、くらいの気持ちに近いかもしれません。

*3:でもレビュアーとしてするべきではありそう。実装の代替案は出さずに反例だけ挙げられたら嫌。

u64 の平方数判定を f64 の sqrt でやるやつの正当性

Rust で言うところの下記のような関数です。

fn is_square(nn: u64) -> bool {
    let n = (nn as f64).sqrt() as u64;
    n.wrapping_mul(n) == nn
}

fn is_square(nn: u64) -> bool {
    ((nn as f64).sqrt() as u64).wrapping_pow(2) == nn
}

これが問題ないということを示そうねというのが今回のお話です。

力技

u64 の範囲の平方数なんてのは $2^{32}$ 個しかないので、それらとその off-by-one を全部テストすることは数秒で可能です。

#[test]
fn boundaries() {
    assert!(is_square(0));
    assert!(is_square(1));
    assert!(!is_square(2));
    for i in 2..=u32::MAX as u64 {
        let ii = i * i;
        assert!(is_square(ii));
        assert!(!is_square(ii - 1));
        assert!(!is_square(ii + 1));
    }
    assert!(!is_square(u64::MAX));
}

æ•°å­¦

数式での証明もしましょう。

IEEE 754 を前提として、sqrt は correctly-rounded で計算できるとします*1。 なので sqrt は問題ないのですが、整数型からの変換 (as f64) の部分で誤差が生じます。 それを考慮した上で、問題なく計算できるというのが面白いところですね。

double (binary64) を前提として、$p = 53$ とします。

Observation 1: $\roundp{\sqrt{94906267^2}} = 94906267$.

ここで、$94906265^2\lt 2^p\lt 94906267^2$ です。

>>> 94906267**2
9007199515875289
>>> 9007199515875289.0
9007199515875288.0
>>> math.sqrt(9007199515875289.0)
94906267.0

$$ \begin{aligned} \roundp{94906267^2} &= 94906267^2-1 \\ \sqrt{94906267^2-1} &= 94906266.{\footnotesize 999999994731644}{\scriptsize 012507624877103}{\tiny 348660069207956{\dots}} \\ \roundp{\sqrt{94906267^2-1}} &= 94906267. \end{aligned} $$

Lemma 2: 実数 $x\gt 0$ に対して、$1+\tfrac12x \gt \sqrt{1+x}$ が成り立つ。

Proof

$f(x) = (1+\tfrac12x)-(1+x)^{1/2}$ とする。$f(0) = 0$ なので、$x\gt 0 \implies f'(x)\gt 0$ を示せば十分。

両辺の差を微分して、 $$ \begin{aligned} f'(x) &= \tfrac12-\tfrac12(1+x)^{-1/2} \\ &= \tfrac12(1-(1+x)^{-1/2}). \end{aligned} $$ $x\gt 0$ のとき $(1+x)^{-1/2} \lt 1$ より従う。$\qed$

note (fact): $\sqrt{1+x} = 1 + {\small 0.5}x - {\small 0.125}x^2 + {\small 0.0625}x^3 - {\small 0.0390625}x^4 + O(x^5)$.

Lemma 3: 実数 $0\lt x\lt 1$ に対して、$1-x\lt\sqrt{1-x}$ が成り立つ。

Proof

$0\lt 1-x \lt 1$ であるから、$(1-x)^1 \lt (1-x)^{1/2}$ は明らか。$\qed$

Lemma 4: $f(x)$ を次で定義する。 $$ f(x) = \frac{\tfrac x2}{(1-\tfrac x2)-\sqrt{1-x}}. $$ このとき、任意の $0\lt x\lt 1$ に対して $f(x) \gt \tfrac 4x-3$ が成り立つ。

Proof

$$ \begin{aligned} f(x) &= \frac{\tfrac x2}{(1-\tfrac x2)-\sqrt{1-x}} \\ &= \frac{\tfrac x2}{(1-\tfrac x2)^2-(1-x)}\cdot( (1-\tfrac x2)+\sqrt{1-x}) \\ &= \frac{\tfrac x2}{\tfrac{x^2}4}\cdot( (1-\tfrac x2)+\sqrt{1-x}) \\ &\gt \tfrac2x\cdot( (1-\tfrac x2)+(1-x) ) \\ &= \tfrac2x\cdot(2-\tfrac32x) \\ &= \tfrac4x-3.\quad\qed \end{aligned} $$

note (fact): $f(x) = \tfrac4x-2-{\small 0.25}x-{\small 0.125}x^2-{\small 0.078125}x^3-{\small 0.0546875}x^4+O(x^5)$.

Therorem 5: 整数 $2^{p/2}\le n\le 2^p$ に対し、$\roundp{\sqrt{\roundp{n^2}}} = n$ が成り立つ。

Proof

この範囲において $\roundp{n}=n$ であることに注意する。 $\roundp{\sqrt{\roundp{n^2}}}\le n$ かつ $\roundp{\sqrt{\roundp{n^2}}}\ge n$ を示せばよい。

( $\le$ ): $\sqrt{\roundp{n^2}} \lt n + \hfloor{n}\cdot 2^{-p}$ を示す。

$$ \begin{aligned} \sqrt{\roundp{n^2}} &\le \sqrt{n^2(1+2^{-p})} \\ &= n\sqrt{1+2^{-p}} \\ &\lt n\cdot(1+2^{-(p+1)}) \\ &\lt n+\hfloor{n}\cdot 2^{-p}. \end{aligned} $$

( $\ge$ ):

Case 1: $\hfloor n=n \implies \sqrt{\roundp{n^2}} \ge n - \hfloor{n}\cdot 2^{-(p+1)}$ を示す。

$\hfloor n=n$ より、ある整数 $k$ に対して $n=2^k$ となる。 すなわち、$\roundp{n^2} = \roundp{2^{2k}} = 2^{2k}$ であり、$\sqrt{\roundp{n^2}} = n$ から、明らかに従う。

Case 2: $\hfloor n\lt n \implies \sqrt{\roundp{n^2}} \gt n - \hfloor{n}\cdot 2^{-p}$ を示す。

$$ \begin{aligned} \sqrt{\roundp{n^2}} &\ge \sqrt{n^2(1-2^{-p})} \\ &= n\sqrt{1-2^{-p}}. \end{aligned} $$

$n\le 2\hfloor{n}-1$ より $\hfloor{n}\ge \tfrac{n+1}2$ なので、 $$ \begin{aligned} n-\hfloor n\cdot2^{-p} &\le n-\tfrac{n+1}2\cdot 2^{-p} \\ &= n\cdot(1-2^{-(p+1)})-2^{-(p+1)}. \end{aligned} $$ よって、下記を示せば十分。 $$ n\sqrt{1-2^{-p}} - \left(n\cdot(1-2^{-(p+1)})-2^{-(p+1)}\right) \gt 0. $$

これは、 $$ n \le 2^p \lt 2^{p+2}-3 \lt f(2^{-p}) = \frac{2^{-(p+1)}}{(1-2^{-(p+1)})-\sqrt{1-2^{-p}}} $$ より従う。$\qed$

Corollary 6: 整数 $0\le n\le 2^p$ に対し、$\roundp{\sqrt{\roundp{n^2}}} = n$ が成り立つ。

Proof

$0\le n\lt 2^{p/2}$ のときは、$\roundp{n^2} = n^2$ より明らか。Theorem 5 より従う。

note: $n=2^p+1$ のときに成り立たないのは明らか。

>>> 2**53+1
9007199254740993
>>> sqrt((2**53+1)**2)
9007199254740992.0

wrapping に関して

ここ追記です。

n.wrapping_mul(n) の部分で変なことにならないことを示しておきます。 n.wrapping_mul(n) が wrap するのは n が$2^{32}$ 以上のときです。一方、入力 nn は $2^{64}-1$ 以下なので、n は $2^{32}$ 以下になります（$\sqrt{2^{64}-1}\le 2^{32}$ かつ $\roundp{2^{32}} = 2^{32}$ より）。

そのため、n が $2^{32}$ となる場合のみ考えればよいです。

n が $2^{32}$ のときは n.wrapping_mul(n) == 0 ですが、nn == 0 と両立することはないので、誤って true を返すことはありません。 また、nn が $2^{64}$ となることは（u64 の範囲から）ないので、誤って false を返すこともありません。$\qed$

あとがき

これで安心して sqrt を使って平方数判定をして生きていけます。 結果が所望の整数になってくれるため、「念のため round しておく」のような処理も不要であることがわかります。

？？「u64::isqrt は stabilize されてますよ笑」

ところで、Theorem 5 の ( $\ge$ ) の Case 2 を示すのに苦戦しました。以前自分で書いた記事を見て「なるほどな〜」となりました。衰えを感じます。悔しかったので少し違う方針で示しました。

rsk0315.hatenablog.com

fact の部分の Taylor 展開は sympy.series などで求めています。「こんな不等式成り立つの？」という気持ちになったときに、とりあえずグラフや Taylor 展開を見て「たしかに成り立ちそう・示せそうか〜」とやっていました。甘えかもしれません。

Gappa や Sollya や SageMath ともお友だちになろうとしているところです。

そのうちまた初等関数についての話を書くと思います。

おわり

おわりです。

*1:IEEE 754 では、四則演算や平方根などをはじめとする各演算は correctly-rounded であることが要求されて (shall) います。一方、C の規格では四則演算を含む各演算の精度は処理系定義ということになっています。もしかして解散ですか？ GCC や LLVM においては問題ないような気もしますが。

cbrt を求めよう (with glibc)

「そろそろ四則演算の誤差評価も飽きたな〜」「数学関数の計算とかもやってみてよ」という声が聞こえてきた*1ので、今日はそれをします。

å°Žå…¥

よく言われているような勘違いについて、一旦触れておきます。

cbrt(x) って pow(x, 1.0/3.0) でよくない？

IEEE Std 754-2019 の 9.2.1 Special values には、下記の記述があります。

$\mathrm{pow}(x, y)$ signals the invalid operation exception for finite $x\lt 0$ and finite non-integer $y$.

実際、pow(-1.0, 1.0/3.0) を計算してみると、C では NaN になり、Python では ValueError: math domain error が出ました。

また、精度においても違いがあります。 1.0/3.0 というのは、$\tfrac13$ とは異なる値です。 $$ \begin{aligned} 1\oslash3 &= 0.{\small 333333333333333}{\footnotesize 314829616256247}{\scriptsize 390992939472198}{\tiny 486328125} \\ &= \tfrac13\cdot(1-2^{-54}). \end{aligned} $$ $\tfrac13\cdot 2^{-54}$ の誤差のせいで、（cbrt が correctly-rounded な値を返すとして実装されていたとしても）所望の値になってくれないケースがしばしば発生します。 $$ \begin{aligned} 64^{\tfrac13\cdot(1-2^{-54})} &= 3.{\small 999999999999999}{\footnotesize 692180816275335}{\scriptsize 216604455276753}{\tiny 748400773143551{\dots}}, \\ \roundp{64^{\tfrac13\cdot(1-2^{-54})}} &= 3.{\small 999999999999999}{\footnotesize 555910790149937}{\scriptsize 383830547332763}{\tiny 671875} \\ &= 4-2^{-51}. \end{aligned} $$

そういうわけで、専用の関数を使う必要があるわけですね。

note: 1.0 / 3.0 を渡された pow 関数側に「オッ cbrt を計算したいんやな、気ィ利かせて補正しといたろｗ」とされると、それはそれでめちゃくちゃになりますね*2。

glibc 実装

精度？

cbrt(3375.0) を計算してみると、14.9999999999999982236431605997495353221893310546875 になりました。ウ〜ンム。 $3375 = 15^3$ ですが、$15-2^{-49}$ が返ってきていて嫌な気持ちになりますね。

note: GCC で単に cbrt(3375.0) と書くと 15.0 が出力される場合がありますが、これはコンパイル時に計算しているためです（アセンブリを見たりするとわかります）。C の規格的には、コンパイル時の計算と実行時の計算で精度などに違いがあることは許容されています。cbrt(strtod("15", NULL)) などとすると、そういう最適化を防げそうです。それでも防げなかった場合は、実行時に与えるようにするとよいですね。

N2347 からの引用

6.6/5

An expression that evaluates to a constant is required in several contexts. If a floating expression is evaluated in the translation environment, the arithmetic range and precision shall be at least as great as if the expression were being evaluated in the execution environment.

5.2.4.2.2/8

The accuracy of the floating-point operations (+, -, *, /) and of the library functions in <math.h> and <complex.h> that return floating-point resuls is implementation-defined, as is the accuracy of the conversion between floating-point internal representations and string representations performed by the library functions in <stdio.h>, <stdlib.h>, and <wchar.h>. The implementation may state that the accuracy is unknown.

そんな無責任なあ。$\eod$

19.7 Known Maximum Errors in Math Functions には、glibc が correct rounding な結果を目指しているわけではない旨や、crsin のように correct rounding 版を用意するかもしれない旨や、既知の誤差の表などが記載されています。

note: GCC (the GNU Compiler Collection) と glibc (the GNU C Library) は別物で、こうした標準関数は GCC ではなく glibc が提供しているものです。GCC はあくまで、（別途用意されているところの）glibc の関数を呼び出す実行ファイルを作っているだけですね (cf. 4.14 Library Functions)。

まえおき

さて、correctly-rounded ではないと知ってやる気を失くしている人もいるかもしれませんが、一応実装を見に行きましょう。 読者の人々（や、昔のえびちゃん）は数値計算のアルゴリズムに明るくないので、「どうせ Taylor 展開とかをしてるんでしょ？」のようなぼやぼやっとしたイメージをしていることと思いますが、そうではないです。

sysdeps/ieee754/dbl-64/s_cbrt.c (glibc-2.41) を読んでいきます。

disclaimer: sysdeps 配下の構造をいまいちわかっていないです。また、数学関数のライブラリのファイル名には e_ k_ s_ t_ v_ w_ などの prefix がついているのですが、これがなにを意味するのかもわかっていないです。glibc 以外の数学関数のライブラリでもこうした名前になっており、そうした慣習がある（or 昔のえらいプロジェクトがそうしていた？）とかなのかなぁと思っています。w_ は wrapper なのかなという雰囲気は感じています。k_ は kernel を意味していそうなのを見た記憶はあります*3。

方針

大まかな方針は次の通りです。

$\gdef\quot{\operatorname{quot}}$ $\gdef\rem{\operatorname{rem}}$

• ${|x|} = x_m\times 2^{x_e}$ として分解する。ここで、$x_m\in[0.5\lldot 1)$ で、$x_e$ は整数である。
• $\sqrt[3]{|x|} = x_m^{1/3}\times 2^{x_e/3} = (x_m^{1/3}\cdot 2^{(x_e\rem 3)/3})\times 2^{x_e \quot 3}$ と変形し、それぞれの部分を求める。
  • $\quot$ と $\rem$ はそれぞれ C の / と % に対応する除算と剰余算。商を $0$ 方向に丸めるため、あまりが負になることもある。

$(x_m, x_e)$ を求める部分は、frexp という標準関数によって行います。これは、内部のビット表現の操作によるものなので、誤差なしで可能です。 この際に、$0_{\pm}$, $\pm\infty$ ã‚„ NaN であることの判定も可能なので、それらの場合はそのまま返してしまいます。

定数の計算

さて、本題となる下記について考えていきます。 $$ \sqrt[3]{|x|} = (x_m^{1/3}\cdot 2^{(x_e\rem 3)/3})\times 2^{x_e \quot 3}. $$ まず、$2^{(x_e\rem 3)/3}$ の部分は $2^{-2/3}$, $2^{-1/3}$, $2^{0/3}$, $2^{1/3}$, $2^{2/3}$ の 5 通りしかないため、定数として用意しておきます。 $$ \begin{aligned} r_{-2} &= 0.{\small 629960524947436}{\footnotesize 484310969717625}{\scriptsize 994235277175903}{\tiny 3203125}, \\ r_{-1} &= 0.{\small 793700525984099}{\footnotesize 680698591328109}{\scriptsize 614551067352294}{\tiny 921875}, \\ r_0 &= 1, \\ r_1 &= 1.{\small 259921049894873}{\footnotesize 190666544360283}{\scriptsize 296555280685424}{\tiny 8046875}, \\ r_2 &= 1.{\small 587401051968199}{\footnotesize 583441787581250}{\scriptsize 537186861038208}{\tiny 0078125}. \end{aligned} $$

$\angled{r_1, r_2} = \angled{\roundp{\sqrt[3]{2}}, \roundp{\sqrt[3]{4}}}$ ですが、$\angled{r_{-2}, r_{-1}} = \angled{1\oslash r_2, 1\oslash r_1} \ne \angled{\roundp{\sqrt[3]{0.25}}, \roundp{\sqrt[3]{0.5}}}$ となっています。

これらの値に関して

sysdeps/i386/fpu/s_cbrt.S でも同様のアルゴリズムがアセンブリで書かれていますが、こちらでは $\angled{r_{-2}, r_{-1}} = \angled{\roundp{\sqrt[3]{0.25}}, \roundp{\sqrt[3]{0.5}}}$ となっており、意図した定義になっているのかは微妙な気がします。コンパイル時計算などで実は correctly-rounded の値になる可能性も否定はできませんが、手元の環境ではそうはなっていないようでした。

volatile double x = 0.125; // コンパイル時に計算されるのを防ぐ
cbrt(x);

として、$0.5$ になれば $r_{-2} = \roundp{\sqrt[3]{0.25}}$ となっています。

volatile double x = 0.25;
cbrt(x);

こちらは、$0.{\small 629960524947436}{\footnotesize 706355574642657}{\scriptsize 302320003509521}{\tiny 484375}$ になれば $r_{-1} = \roundp{\sqrt[3]{0.5}}$ となっています。なお、このときの cbrt(0.25) の値は $\roundp{\sqrt[3]{0.25}}$ とは異なっていることに注意してください 🥲。

$$ \begin{aligned} \roundp{\sqrt[3]{0.5}} &= 0.{\small 793700525984099}{\footnotesize 791720893790625}{\scriptsize 268593430519104}{\tiny 00390625}, \\ \roundp{\sqrt[3]{0.25}} &= 0.{\small 629960524947436}{\footnotesize 595333272180141}{\scriptsize 648277640342712}{\tiny 40234375}. \quad\eod \end{aligned} $$

近似多項式

さて、残りは $x_m^{1/3}$ を求めるパートです。 定義から $x_m\in[0.5\lldot1)$ です。 このように、入力が特定の区間に含まれるケースに帰着する手法は、これ以外の数学関数においてもしばしば見られます。

以下で定義される多項式 $f$ を考えます。 $$ \begin{aligned} f(x) &= 0.{\small 354895765043919}{\footnotesize 841907893442112}{\scriptsize 253978848457336}{\tiny 42578125} \\ &\qquad {}+ 1.{\small 508191937815849}{\footnotesize 037463067361386}{\scriptsize 492848396301269}{\tiny 53125}\,x \\ &\qquad {}- 2.{\small 114994941673713}{\footnotesize 046981220031739}{\scriptsize 212572574615478}{\tiny 515625}\,x^2 \\ &\qquad {}+ 2.{\small 446931225635344}{\footnotesize 375746171863283}{\scriptsize 962011337280273}{\tiny 4375}\,x^3 \\ &\qquad {}- 1.{\small 834692774836130}{\footnotesize 801929812150774}{\scriptsize 523615837097167}{\tiny 96875}\,x^4 \\ &\qquad {}+ 0.{\small 784932344976639}{\footnotesize 218003811038215}{\scriptsize 644657611846923}{\tiny 828125}\,x^5 \\ &\qquad {}- 0.{\small 145263899385486}{\footnotesize 366933051272098}{\scriptsize 964545875787734}{\tiny 9853515625}\,x^6. \end{aligned} $$

この $f$ に対して $x_m^{1/3} \approx f(x_m)$ が成り立ちます。 $f(x_m)$ の計算の際には Horner’s method を用います*4。$f$ の計算にも誤差が出るため、Horner’s method による $f(x_m)$ の値は $\hat f(x_m)$ と表記することにしましょう。

この多項式は、おそらくは Remez algorithm などで求めたのだろうと推測しています。特定の区間における誤差の最大値を最小化するような（予め決めた次数以下の）近似多項式を求めるアルゴリズムです。 $x_m\in[0.5\lldot 1)$ としているので、$[0.5\lldot 1)$ でのよい近似が得られていればよいですね。

f の近似の様子のグラフ

点線がその境界、破線が近似されるべき真の値 $y=x^{1/3}$ で、赤の実線が上記の多項式 $y=\hat f(x)$ です。よさそうな見た目をしていますね。

念のためもう少し近づけておく

$\gdef\signum{\operatorname{sgn}}$

さて、もう少しだけ続きます。$\signum(x)\times(\hat f(x_m)\otimes r_{x_e\rem 3})\times 2^{x_e\quot 3}$ を返すわけではないです。

返すのは、$u = \hat f(x_m)$ および $t_2 = u\otimes u\otimes u$ として、次の値です。 $$ \signum(x)\times (u \otimes (t_2 \oplus 2x_m)\oslash(2t_2\oplus x_m)\otimes r_{x_e\rem 3})\times 2^{x_e\quot 3} $$

やっぱり $\circledast$ の表記だと分数が使えないので読みにくい 🥲。 計算せんとしている式の絶対値は次の通りです（$\signum(x)$ の部分は煩雑になるので省きたいです）。 $$ u\cdot\frac{u^3+2x_m}{2u^3+x_m}\cdot r_{x_e\rem 3}\times 2^{x_e\quot 3} $$

$u \xgets{\times} \tfrac{u^3+2x_m}{2u^3+x_m}$ で更新し、$u\cdot r_{x_e\rem 3}\times 2^{x_e\quot 3}$ を返そうとしているという見方もできます。というわけで、この更新について見ていきます。

これは、Newton’s method ではなく Halley’s method と呼ばれる反復アルゴリズムの更新ステップです。次のような更新をします。 $$ x_{n+1} = x_n - \frac{2g(x_n)g'(x_n)}{2g'(x_n)^2-g(x_n)g''(x_n)} $$ この反復により、$g(x) = 0$ の根を求めます。 たとえば $a$ の 3 乗根を求めたいときは $g(x) = x^3-a$ として行えばよいですね。 $g'(x) = 3x^2$, $g''(x) = 6x$ です。

$$ \begin{aligned} x_{n+1} &= x_n - \frac{2\cdot(x_n^3-a)\cdot 3x_n^2}{2\cdot(3x_n^2)^2-(x_n^3-a)\cdot 6x_n} \\ &= x_n - \frac{6x_n^5 - 6ax_n^2}{18x_n^4 - 6x_n^4+6ax_n} \\ &= x_n - \frac{x_n^4 - ax}{2x_n^3+a} \\ &= \frac{2x_n^4+ax_n-x^4+ax_n}{2x_n^3+a} \\ &= x_n\cdot \frac{x_n^3+2a}{2x_n^3+a} \end{aligned} $$ $(x_n, a) = (u, x_m)$ とすることで、上で見た式と一致します。

$u^3$ ã‚’ t2 という名前にした意図はよくわかりませんでした。何らかの一般的な手法で使われていた記号を流用したとかでしょうか。

完成

$\signum(x)\times u\cdot r_{x_e\rem 3}\times 2^{x_e\quot 3}$ を返します。

最後の $2^{x_e\quot 3}$ を掛けるパートは、ldexp という標準関数により行います。

余談

だったら $f$ なんて使わずに最初から Halley’s method を使えばいいんじゃんという声も聞こえます。 おそらくは除算が重いため避けているのではないかなと推測するところです。

(xm, xe) = math.frexp(x)
yi = math.ldexp(xm, xe // 3)
for _ in range(4):
    yi3 = yi * yi * yi
    yi *= (yi3 + 2 * xi) / (2 * yi3 + xi)

    # yi3 = yi * yi * yi
    # yi3_xi = yi3 + xi
    # yi *= (yi3_xi + xi) / (yi3_xi + yi3)

反復 4 回の収束の様子

$(0\lldot 2]$ の範囲で $10^{-5}$ 刻みで試したところ、4 回の反復でだいたい $5\times 10^{-16}$ 程度の相対誤差に収まっていそうでした。 反復が 3 回のときは $6\times 10^{-6}$ くらいでした。3 次収束らしいので、1 回の反復でだいぶ差が出ますね。

反復 1 回あたり +, *, / がそれぞれ 2, 5, 1 回（計算順序を変えて yi3 + xm 先に計算することにしても 3, 3, 1 回）です。 Horner’s method による $\hat f$ の計算回数はそれぞれ 6, 6, 0 回（- ã‚‚ + として計上）なので、Halley’s method 2 回分よりもお得ですね。適当な初期値から反復するよりも、多項式を用いて計算する方がお得です。

Halley’s method で収束した後の値に対して $r_{x_e\rem 3}$ を掛けるのは、損しているような気がしないでもないような気もします。

また、Newton’s method で yi *= (2 * yi3 + xi) / (3 * yi3) を用いた場合は、反復を 7 回しても $4\times 10^{-12}$ 程度の誤差がありました。

参考文献

• Muller, Jean-Michel, and Jean-Michael Muller. Elementary functions. BirkhÅ©ser Boston, 2006.

感想

今回は glibc の cbrt を見ましょうというだけの回でしたが、やっぱりそれだけじゃ満足できなくて、（使い道がどの程度あるかなんてのはさておき）correctly-rounded な値を計算できることに憧れますよね。 競プロ界隈においては（というか人間界においては）、「そんな標準関数の精度を気にするより、普通の四則演算で catastrophic cancellation を起こさないように気をつける方が先ですよ」というケースの方が多いので、correctly-rounded なんてのはまだ先の先の話かもしれないですが、人間界の事情なんてえびちゃんの知ったことではないですね。

近似多項式 $f$ も天下りで与えたものを使っただけなので、ちょっと物足りなさはあります。 Remez algorithm についても深掘りできたらいいなと思っています。

C や IEEE 754、POSIX に C++ と、さまざまな規格で思い思いの数学関数を定義しておられるので、めちゃくちゃな個数の数学関数があります。特に C++ はもうめちゃくちゃです。

数学関数たち

C の規格で定義されている数学関数には下記があります。ビット操作やそれに類するものなどは除外しています。

• trigonometric functions
  • cos $\cos(x)$, sin $\sin(x)$, tan $\tan(x)$
  • acos $\arccos(x)$, asin $\arcsin(x)$, atan $\arctan(x)$
  • atan2
• hyperbolic functions
  • acosh $\arcosh(x)$, asinh $\arsinh(x)$, atanh $\artanh(x)$
• exponential and logarithmic functions
  • exp $e^x$, exp2 $2^x$, expm1 $e^x-1$
  • log $\ln(x)$, log10 $\log_{10}(x)$, log1p $\ln(1+x)$, log2 $\log_2(x)$
• power and absolute-value functions
  • cbrt $x^{1/3}$
  • fabs $|x|$
  • hypot $\sqrt{x^2+y^2}$
  • pow $x^y$
  • sqrt $\sqrt x$
• error and gamma functions
  • erf $\operatorname{erf}(x) = \displaystyle\frac2{\sqrt\pi}\int_0^x e^{-t^2}\,\mathrm{d}t$
  • erfc $\operatorname{erfc}(x) = 1 - \operatorname{erf}(x)$
  • lgamma $\ln(|\Gamma(x)|)$
  • tgamma $\Gamma(x)$
• nearest integer functions
  • ceil $\ceil{x}$
  • floor $\floor{x}$
  • round, roundeven, trunc
• floating multiply-add
  • fma $(x\times y)+z$

IEEE Std 754-2019 の additional mathematical operations には、下記のような関数もあります。

• exp2m1 $2^x-1$, exp10 $10^x$, exp10m1 $10^x-1$
• log2p1 $\log_2(1+x)$, log10p1 $\log_10(1+x)$
• rSqrt $1/\sqrt{x}$
• compound $(1+x)^n$
• rootn $x^{1/n}$
• sinPi $\sin(\pi x)$, cosPi $\cos(\pi x)$, tanPi $\tan(\pi x)$
• asinPi $\arcsin(x)/\pi$, acosPi $\arccos(x)/\pi$, atanPi $\arctan(x)/\pi$
• atan2Pi

POSIX には次の関数もあります。

• j0, j1, jn — Bessel functions of the first kind
• y0, y1, yn — Bessel functions of the second kind

C++ にもたくさんあります (cf. [cmath.syn])。列挙を諦めました。ところでこれいつ使うんですか？ $\eod$

え、これ全部の correctly-rounded 版を実装させられるんですか...？ 😭（correctly-rounded の本質は、実装ではなく証明パートという話もあります。）

おわり

おわりです。

*1:自分の中から聞こえてきた声です。そもそも浮動小数点型に興味を持っているフォロヮが少ないため。

*2:現実としては、そういうめちゃくちゃな使い方をしたい人が多そうなので、そういう仕様の方が「幸福度」が高くなるのかもしれませんが...。そういう ad hoc な意味不明仕様にはなってほしくないですね。

*3:たぶん Linux kernel とかの kernel とは違う文脈で、「核となる部分だよ〜」くらいだと解釈しています。

*4:$a_0+a_1x+a_2x^2+\cdots$ を計算する際に、$a_0\oplus (a_1\otimes x)\oplus(a_2\otimes x\otimes x) \oplus \cdots$ とするのではなく、$( (\cdots\oplus a_2)\otimes x\oplus a_1)\otimes x\oplus a_0$ とする手法。