いろいろな数学関数たちの 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_-$ としておきます。
下記の記事でも多少触れています。
誤差評価の過程ではゼロの符号は問題にならないので、数式中では単に $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 な実装の典型的なパターンとして、次のような流れがあります。
- range reduction
- 入力を一般のケースから特殊なケースに帰着させる
- approximation
- 所望の関数に応じた手法を用い、特殊なケースについての答えを求める
- 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$ で更新する。
- $\roundcirc{(y\oplus\Delta y)\times(y\oplus\Delta y)+(-x_m)}\le 0$ であれば下記を行う。
- $y$ を出力する
note: 各ループの先頭において、$\Delta y = 0.5^i$ が成り立ちます。
Claim 1: 任意の浮動小数点数 $x\in[1\lldot 2)$, $y\in[1\lldot 2]$, $z\in[1\lldot 4)$ に対し、 $$ \begin{aligned} x\times y\le z &\iff \roundcirc{x\times y+(-z)}\le 0, \\ x\times y\lt z &\iff \roundcirc{x\times y+(-z)}\lt 0 \end{aligned} $$ が成り立つ。
Proof
任意の実数 $w$ に対し、$\roundcirc{w}\le 0 \iff w\le 2^{-1075}$ および $\roundcirc{w}\lt 0 \iff w\lt -2^{-1075}$ である。 よって、 $$ \begin{aligned} x\times y\le z &\iff x\times y-z \le 2^{-1075}, \\ x\times y\lt z &\iff x\times y-z \lt -2^{-1075}, \\ \end{aligned} $$ を示す。
( $\le$, $\implies$): 明らか。
( $\le$, $\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}.$$
( $\lt$, $\implies$): 対偶 $x\times y-z \ge -2^{-1075} \implies x\times y\ge z$ を示す。
$x\times y-z$ は $2^{-104}$ の倍数であるから、 $$x\times y-z \ge -2^{-1075}\implies x\times y-z\ge 0.$$
( $\lt$, $\impliedby$): 明らか。$\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 を用いることで、もっと簡単に求めることができます。
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:二分探索で済むなら他のことを学ぶ理由がないと思っている人々が多数派だから（偏見）。
