小清水さんとコンピューター数学

コンピューター・数学 に関することを書きます (特に丸め誤差の話が多いです。)

エラーフリー変換の紹介 および FastTwoSum アルゴリズム の紹介と証明 -- glibc のコードを読むための参考に --

小清水 (@curekoshimizu) です。

仕事が忙しくなかなか更新する余力がありませんでしたが一年半振りにブログを更新します!

エラーフリー変換とは?

浮動小数点数の計算というのは、このブログで何度も紹介していますが、 丸め誤差 を伴います。

簡単な  x+y という計算であっても  x + y を行ったあとに浮動小数点に変換されるために丸め誤差が発生する可能性があるのです。 その他にも

  •  xy
  •  xy + z
  •  xy + zw
  •  \sum_{n=1}^N a_n
  •  \sum_{n=1}^N a_n b_n

などなど各種計算も同様に計算過程で丸め誤差をともないます。

つまり

(求めたい計算式) = (愚直に浮動小数点体系で計算して得た結果) + (丸めによる誤差)

例えば

 x = 0.300000011920928955078125  = 1.20000004768371581677722\times 2^{-2}  y = 0.20000000298023223876953125 = 1.6000000238418585033165\times 2^{-3}

において  x + y を float32 で計算すると  0.5 = 1.0 \times 2^{-1}

となり、誤差  0.00000001490116119384765625 が生じます。

この 誤差項 も含めて結果を得たいというケースがあるのです。

これを行うのが エラーフリー変換 (Error Free Transformation) と呼ばれる計算です。

こんなこと本当にしたいのでしょうか?

一つの例としてみなさんがお世話になっている glibc でも使われているのだというのを見てみましょう。

glibc でエラーフリー変換が使われている例

たとえば glibc では ソフトウェアFMA を計算するために次のような処理をしています。

https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/s_fma.c の LN188 - LN205 あたり。

  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
  #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
  double x1 = x * C;
  double y1 = y * C;
  double m1 = x * y;
  x1 = (x - x1) + x1;
  y1 = (y - y1) + y1;
  double x2 = x - x1;
  double y2 = y - y1;
  double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  double a1 = z + m1;
  double t1 = a1 - z;
  double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  double a2 = t1 + t2;

github.com

ここで登場するコメントにも書かれている Dekker の乗算アルゴリズムKnuth の加算アルゴリズム

はエラーフリー変換の一つです。

  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
  #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
  double x1 = x * C;
  double y1 = y * C;
  double m1 = x * y;
  x1 = (x - x1) + x1;
  y1 = (y - y1) + y1;
  double x2 = x - x1;
  double y2 = y - y1;
  double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;

例えばこの箇所は x*y を計算する m1 = x * y という計算によってどの程度の誤差が生じたのかを知るために m2 を厳密に計算しています。

これにより

m1 + m2 = x * y

となることを数学的に保証するというものです。

こちらの加算も同様で

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  double a1 = z + m1;
  double t1 = a1 - z;
  double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  double a2 = t1 + t2;

a1 = z1 + m1 で求めたい加算結果を得て、それによって生じた誤差を a2 で手に入れています。

乗算に関するエラーフリー変換と FMA 命令

上の Dekker の乗算アルゴリズムFMA 命令がないアーキテクチャのためのアルゴリズムとして知られています。

もし FMA命令 (xy+zを計算する命令) を持っている場合には

  /* Multiplication m1 + m2 = x * y using FMA instruction.  */
  double m1 = x * y;
  double m2 = x * y - m1;

これでエラーフリー変換できてしまいます。 そのため、現代のアーキテクチャで、 この Dekker の乗算アルゴリズムを実際に使う機会は少ないかもしれません。

FMA についてはこちらで詳しく記述しました:

math-koshimizu.hatenablog.jp

加算に関するエラーフリー変換

一方で加算はそうはいきません。FMA命令では解決できないのです。

歴史 -- FastTwoSum(Fast2Sum)・TwoSum アルゴリズム

1971 年に Dekker に紹介された FastTwoSum アルゴリズム というアルゴリズムがあります。(T.J. Dekker: A floating-point technique for extending the available precision)

これが加算に関するエラーフリー変換として最も古いアルゴリズムです。 とはいえ、

1965 年に Kahan が総和演算に関する Compensated summation アルゴリズムに関する論文の中でこのアルゴリズムを暗に利用しているのも事実ですので、 最古といっていいか悩ましいところはあります…。(https://convexoptimization.com/TOOLS/Kahan.pdf)

この Dekker のアリゴリズムを FastTwoSum と名付けたのは Shewchuk の 1997年の論文です:

f:id:curekoshimizu:20190421191324p:plain
Shewchuk : Adaptive precision floating=point arithmetic and fast robust geometric predicates

ただし、この FastTwoSum は上記通り、条件がついており、引数a, bについて大小比較をする必要 があります。

このため、現在のアーキテクチャの高速化技法等を考えると、条件分岐が発生するこのアルゴリズムよりも、

KunuthとMøllerによる条件分岐の発生しない加算エラーフリー変換である、 TwoSum アルゴリズム のほうが一般的に利用されています。

上の glibc で書かれていた

  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
  double a1 = z + m1;
  double t1 = a1 - z;
  double t2 = a1 - t1;
  t1 = m1 - t1;
  t2 = z - t2;
  double a2 = t1 + t2;

この式はこのアルゴリズムです。

命名は先ほど登場した shewchuk の論文です。

論文中に FastTwoSum(a, b), Fast2Sum(a, b), TwoSum(a, b), 2Sum(a, b) などと登場したらこれらのアルゴリズムのことなんだなと思うと良いと思います。

FastTwoSum アルゴリズムの紹介と証明

今回のブログではこの古典的な 加算エラーフリー変換アルゴリズムである Dekker の FastTwoSum アルゴリズムについて紹介と証明を与えてみようと思います。

FastTwoSum アルゴリズム

次の計算により |a| \geq |b| の場合に  a+b = s + t = \mathrm{RN}(a+b) + b なる  s, t が得られます。ただしここでは2進数浮動小数点環境とする:

Dekker FastTwoSum
 s = \mathrm{RN}(a + b)
 z = \mathrm{RN}(s - a)
 t = \mathrm{RN}(b - z)

ただし  \mathrm{RN} は最近接偶数方向丸めを表します。この最近接偶数方向丸めについては次の記事を参照ください:

math-koshimizu.hatenablog.jp

FastTwoSum アルゴリズム証明

ここで

  •  a = M_a 2^{e_a}
  •  b = M_b 2^{e_b}
  •  s = M_s 2^{e_s}  = \mathrm{RN}(a+b)

とおく、ただし  M_a, M_b, M_s, e_a, e_b, e_s \in \mathbb{Z} かつ   |M_a|, |M_b|, |M_c| <  2^p ,  e_{min} \leq  e_a, e_b, e_s \leq e_{max} とする。

ここで、次のブログの Prop. (整数を用いた  \beta p 桁の浮動小数点数の表示) を使用した。 この証明ではこの性質を頻繁につかっている。

math-koshimizu.hatenablog.jp

また、

 |a + b| \leq |a| + |b| \leq 2|a| =  2 |M_a| 2^{e_a} =  |M_a|  2^{e_a+1}

ここから、 M_s e_s のとり方の自由度はあるものの、  e_s \leq e_a + 1 と制限した上で  M_s を定めても問題ないことがわかる。

ここで次の3つにわけて考える

  • Case1.  e_s > e_a (すなわち  e_s = e_a + 1)

  • Case2.  e_s \leq e_a かつ  e_s > e_b

  • Case3.  e_s \leq e_a かつ  e_s \leq e_b

このとき、

 s-a浮動小数点として正確に表現できることを示し、 s-a = (整数)\times 2^{e_b-p+1} と表せることを示す。

Case1.  e_s = e_a + 1 のとき

このとき

 a + b = M_a {2}^{e_a} + M_b {2}^{e_b}
= \left(  \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}}        \right) \times 2^{e_a+1}
= \left(  \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}}        \right) \times 2^{e_s}

となり、


\left|  \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}}  \right|
 \leq   \frac{|M_a|}{2}  + \frac{|M_b|}{2} 
 \leq | M_a | < 2^{p}

と、 \left(  \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}}        \right) 2^p で抑えられる実数として表現される。

 a +b 最近接偶数方向丸めが  s であるから、先の項を整数化することになる。よって

 
\left| M_s - \left(  \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}}  \right)\right| \leq \frac{1}{2}

と評価でき、

これより次の式が得られる:


| 2M_s - M_a | \leq 2 \left| M_s - \left(  \frac{M_a}{2} + \frac{M_b}{2^{e_a - e_b+1}}  \right)\right| + 2\left| \frac{M_b}{2^{e_a - e_b+1}} \right|
\leq 1 + | M_b | \leq 2^p

ゆえに


s - a = M_s 2^{e_a+1} - M_a 2^{e_a} = (2 M_s - M_a ) 2^{e_a}

から、この  2 M_s - M_a は整数であって、先の評価から  2^p 未満の場合、 s-a浮動小数点表示される。もちろん  2^{p+1}に一致する場合も、指数部を  2^{e_a+1} に修正しても、整数部は丸めによる誤差は生じないので、浮動小数点表示されていることがわかる。

そしてこの表示から  e_b \leq e_a より  s-a = (整数)\times 2^{e_b} と表されることもわかる。

Case2.  e_s \leq e_a かつ  e_s > e_b のとき

ここで

 a + b = M_a 2^{e_a} + M_b 2^{e_b} 
= \left(  2^{e_a - e_s} M_a + \frac{M_b}{ 2^{e_s-e_b}} \right) 2^{e_s}

最近接遇数方向の定義より

 \left|M_s- \left(  2^{e_a - e_s} M_a + \frac{M_b}{ 2^{e_s-e_b}} \right) \right | \leq \frac{1}{2}

ここから、

 |s-a| = | M_s 2^{e_s} - M_a2^{e_a}| 
\leq \left| M_s - \left(  2^{e_a - e_s} M_a + \frac{M_b}{ 2^{e_s-e_b}} \right) \right| 2^{e_s} +  \frac{|M_b|}{2^{e_s-e_b}} 2^{e_s}
\leq \left(  \frac{1}{2} + \frac{|M_b|}{2^{e_s-e_b}}  \right) 2^{e_s}

が得られる。  s-a = \left( M_s 2^{e_s} - M_a 2^{e_a-e_s} \right) 2^{e_s} =: K 2^{e_s} と整数部を  K と表すと

 |K| \leq    \frac{|M_b|}{2^{e_s-e_a}} + \frac{1}{2}  \leq \frac{2^p - 1}{2} + \frac{1}{2} < 2^p

となることから、 s-a浮動小数点表示されていることがわかる。また上評価より指数部は  2^{e_s} であることがわかっており、  e_s > e_b であるから  s-a = (整数)\times 2^{e_b} と表されることもわかる。

Case3.  e_s \leq e_a かつ  e_s \leq e_b のとき

 a + b = (2^{e_a-e_s} M_a + M_b 2^{e_b-e_s}) 2^s

ここで  e_s \leq e_a かつ  e_s \leq e_b であるから  (2^{e_a-e_s} M_a + M_b 2^{e_b-e_s}) は整数となる。

また、  a+b の丸めが  s = M_s 2^{s} であるから、

 (2^{e_a-e_s} M_a + M_b 2^{e_b-e_s}) が整数であることや、 2^{e_s} を比較して考えると、

丸めの定義を考えればそれによって値が変わることはなく、

 M_s = 2^{e_a-e_s} M_a + M_b 2^{e_b-e_s} となることがわかる。

ゆえに、  a+b = s が示され、

 s-a = b であり、 s-a浮動小数点表示されていることがわかる。また、  s-a = (整数)\times 2^{e_b} であることもわかる。

ここまでわかったことまとめ

 s-a浮動小数点表示でき、 s-a = (整数) \times 2^{e_b} となることがわかった。

ゆえに

Dekker FastTwoSum
 s = \mathrm{RN}(a + b)
 z = \mathrm{RN}(s - a)
 t = \mathrm{RN}(b - z)

 z = \mathrm{RN}(s-a) = s - a であることがわかった。

最後に

 t = \mathrm{RN}(b - z) = b - z であることを示せれば、  s + t = s + (b - (s - a) ) = a + b とエラーフリー変換であることが示される。

  t = \mathrm{RN}(b - z) = b - z の証明

 s = \mathrm{RN}(a + b) であるから

 |(a + b) - s| \leq | (a + b) - a| = | b |

である。これは、 s a+b に関して浮動小数点の中で最良近似であるからである。

ゆえに   |b-z| = | a + b -s| \leq |b| が示された。

ここで  s -a = (整数)\times 2^{e_b} と表すことができたので その丸めである  z = \mathrm{RN}(s-a) (整数)\times 2^{e_b} と表される。これは、丸めが入るかもしれないが、整数部が変化するのみであることからわかる。

よって  b -z についても同様であり、  b - z =: L 2^{e_b} と整数 L を定義すれば

 |b-z| \leq |b| であることから  |L| \leq |M_b| <   2^{p} であり、  b-z浮動小数点表示できる。

すなわち    t = \mathrm{RN}(b - z) = b - z が示され、証明完。

まとめ

今回のブログでは エラーフリー変換 について紹介をしました。

そのエラーフリー変換として有名で glibc でも使われている 加算と乗算のエラーフリー変換を紹介しその歴史を紹介し、

さらに、そのエラーフリー変換の証明方法の例として古典的な Dekker のFastTwoSum を紹介しました。

この記事を書くにあたって、 kashi さんのブログを拝見し、TwoSum アルゴリズムで 気をつけるべき例 (途中でオーバーフローが発生する場合にエラーフリーとならないこと) が紹介されていて面白いなと思いました。

verifiedby.me

いつか私も TwoSum アルゴリズムの証明を読んでみたいと思います。

(追記)

調べると TwoSum アルゴリズムは、アンダーフローやオーバーフローが途中で起こらない限りという条件がついていることがわかりました。

そのため、Dekker の FastTwoSum に比べて条件が厳しいことがわかりました。

参考:

f:id:curekoshimizu:20190423121211p:plain
http://perso.ens-lyon.fr/jean-michel.muller/Conf-Journees-Knuth.pdf

FMA (Fused Multiply-Add) について色んな観点でまとめてみた

小清水 (@curekoshimizu) です。

本日は FMA についてお話したいと思います。

FMA とは?

FMA とは Fused Multiply-Add ことで

 \circ (x\times y + z)

の演算のことです。 ここで  \circ (x) x \in \mathrm{R} の丸めを表しました。

本当にただこれだけの内容なのですが、 今回の記事は、この FMA について熱く書いてみたいと思います。

FMA について書くモチベーションなのですが、

本ブログは 精度に関する話題 を多く取り上げてきました。

その中で FMA と関係する話題が非常に多く登場し、

これからも登場予定 です。

そのたびに、 FMA について補足すべきことが多く、 ここでまとめておこうと思い立ちました。

例えばこの記事で FMA 命令と丸め誤差の話がすでに登場しています:

math-koshimizu.hatenablog.jp

この記事を読むと

1. FMA の凄さがわかる

精度や高速化と関わることがわかります。

2. 自分の環境で FMA 命令がサポートしているかがわかる

FMA 命令を積んだ HW の歴史についてある程度まとめてみました。

3. サポートがなくとも無理やり FMA 演算を実行する方法がわかる

FMA 命令を積んでいなくとも FMA 計算する方法がありますので紹介します。

4. コンパイルして FMA 命令を出す方法がわかる

FMA 命令をもっていても、実際に使われなくては意味がありません!

5. FMA の応用がわかる

これと関係して、 Horner 法の起源論文なども載せていたりしますので、 このあたりに興味ある方も是非ご覧ください。

具体的に何がすごいの?

2つのメリットがあります!

1. 丸め回数が減るので精度的に有利

2. HWサポートがある場合に速度的に有利

すごい点その1 – 丸め回数が減るので精度的に有利

上の  x \times y + zFMA 命令を使わずに計算すると

\circ \left( \circ (x\times y) + z \right)

と加算命令1回、乗算命令1回の演算で計算されることになります。

つまり、 2回の丸め が発生します。

それが FMA 命令だと

 \circ (x\times y + z)

1回の丸めしか入らないため、精度的な観点で有利 になります。

すごい点その2 – HWサポートがある場合に速度的に有利

この FMA は多くの最新 CPU で FMA 命令を HW サポートしています!

さらに CPU だけでなく GPU などの HW でもサポートしていることが多いです。

これにより

 \circ (x\times y + z)

という演算を 1回で計算できます。 もっていなければ 乗算命令と加算命令の 2回で 実行する必要があるわけです。 (丸めによる誤差は度外視すると)

これがどういった意味をもつのでしょうか?

この FMA の力をみるべく、 Geforce GTX 1080

という GPU を例に 理論ピーク性能を算出してみましょう。

Geforce GTX 1080 の能力は

  • 2560 CUDA コア
  • 1.733 GHz

で動作し、この CUDA コア1つで 単精度の FMA 命令を発行できます。

FMA は 乗算・加算の 2つの命令分なので 2回の浮動小数点が計算できるという定義になります。

これより 単精度理論ピーク性能

2560 (CUDA コア)  \times 1.733 (GHz)  \times 2 (FLOPS/(CLOCK  \times CUDAコア)) = 8872 GFLOPS

となります。

FMA 命令をもっていないとこの 「2」 の乗算がなくなりピーク性能指標は、半分になってしまいます。

余談になりますが、この GTX 1080 は単精度で 8.8TFLOPS の力を持っているということで、 2017年現在の GPU はこんなにも能力が高いのか…と思わされます。

このように、 性能指標にも活躍するのが FMA 命令です! すなわち、 この命令を使わないと事実上、ピーク性能の半分しか理論上だすことができません。

このため、FMA 命令をちゃんと使ってあげることが、 HW の力を引き出しているかどうかと密接な関係があります。

FMA 命令 HW サポートしてる?

ただしこの FMA、 多くの CPU でサポートと書きましたが、 比較的古い CPU を使っている方はサポートされていない可能性 があります。

私の デスクトップ・ノートPC は次のような環境です:

(デスクトップPC)

  • CPU : Core i7-3820 CPU @3.60GHz
  • メモリ : 64GB (DDR3-1600)

(ノートPC)

  • CPU : Core i5-4288U CPU @ 2.60GHz
  • メモリ : 8GB (DDR3-1600)

こちら、メモリだけでみると デスクトップの性能がすごくいいです!

しかしながら FMA 命令をサポートしているかどうかでいうと、

ノートPC は FMA命令をサポート しており、デスクトップは 非サポート であります。

そのため、FMA を試す環境としては、私のデスクトップ環境は非常によろしくないです。

このあたりの FMA が載っている載っていないの境界は 2011〜2013 年ごろの CPU にあり、 それについても後述します。

また FMA が載っていなくとも、 FMA の計算、すなわち

 \circ (x\times y + z)

の計算をできます!

これは HW サポートをしていないのでかなり遅いですが、 精度上の恩恵を受けるために使うことができる という意味です。

これについても 後述 します。

FMA とその歴史

FMA が載っているかどうか怪しいという話を上でしましたが、 何時頃登場したのでしょうか?

最初に FMA 命令を積んだのは 1990年 の IBM RS/6000 です。

IBM RS/6000 に関する論文 「Design of the IBM RISC System/6000 floating-point execution unit」 には次のように書かれています:

The IBM RISC System/6000® (RS/6000) floating-point unit (FPU) exemplifies a second-generation RISC CPU architecture and an implementation which greatly increases floating-point performance and accuracy. The key feature of the FPU is a unified floating-point multiply-add-fused unit (MAF) which performs the accumulate operation (A times B) + C as an indivisible operation.

このように、性能と精度を向上するべく搭載されたのが FMA です。 この当時、この命令は MAD (multiply-add-fused) と略されていたようです。

その後、浮動小数点の規格として採択されたのは 2008 年のことになります。 IEEE754-2008 に FMA について厳格な仕様が定められました。

その仕様が定まる前にも既に

などで実装されていました。

そのため、2008 年頃にはだいたい CPU は FMA 命令もっているのかな?と思います。

歴史的に気をつけたいのが、 我々のデスクトップ環境やノートPCなどで使われている x86_64Intel CPU や AMD CPU です。

対応されたのはまあまあ最近の話 だということがわかります。 (最近の感覚がずれていたらすみません)

FMAIntel CPU

これらは Intel CPU であれば Haswell になって、 ようやく FMA サポートしています。

具体的には Intel Core シリーズを最近までの表を下に載せてみますと

  • 第1世代: Nehalem
  • 第2世代: Sandy Bridge
  • 第3世代: Ivy Bridge
  • 第4世代: Haswell (2013年頃)
  • 第5世代: Broadwell
  • 第6世代: Skylake
  • 第7世代: Kaby Lake (2017年現在頃)
  • 第8世代: Coffee Lake

となっており 第4世代 Haswell から FMA 命令が HWサポートされました。

(デスクトップPC)

  • CPU : Core i7-3820 CPU @3.60GHz
  • メモリ : 64GB (DDR3-1600)

(ノートPC)

  • CPU : Core i5-4288U CPU @ 2.60GHz
  • メモリ : 8GB (DDR3-1600)

先ほどの環境の、 Core i7-3820」 は第2世代 Sandy Bridge (正確には Sandy Bridge-E) であり FMA 命令をもっていません。 一方でノートPCに載っていた Core i5-4288U」第4世代 Haswell のことであります。

FMAAMD CPU

一方、 AMD CPU であれば AMD FX、Bulldozer から FMA 命令がサポートされており、 最近の何かと話題な AMD Zen シリーズの Ryzen もサポート済みです。

このあたりは細かいことを言うと、 さらにもう少しややこしく、 レジスタをいくつ使うかという点で異なる FMA3・FMA4 の どちらを主流とするかの歴史バトルが IntelAMD でありました。

FMA3 とは例えば

 x = \circ (x\times y + z)
 x = \circ (y\times z + x)

という3レジスタの命令です。一方で FMA4

 w = \circ (x\times y + z)

という4レジスタの命令です。

Intel は FMA3 しか最初からサポートしていませんでしたが、 AMD の最初に FMA サポートした Bulldozer (2011年頃) は FMA4 のみをサポートしていました。

最近登場した Ryzen は FMA4 を非サポートとし、 FMA3 のみをサポート としていますから、

FMA3 のほうが主流となりそうです。

さて、上記のような FMA 命令の中のさらに分類の話はおいておくとして、

結局のところ 2011〜2013年付近の x86_64 頃から、大雑把に FMA 命令が使えるようになってきた ということになります。 (具体的には Intel なら HaswellAMD ならBulldozer)

ここで誤解しないで頂きたいのは、 Intel 初の FMA 搭載 CPU が Haswell というわけではありません。 例えば、 Itanium という IA-64 が 2001 年頃にありまして、 既に FMA 命令もっていました。

加えて、この年代を超えていたら必ず FMA を持っているというわけではありませんので、 自身の CPU を調べてみてください! 例えば Linux であれば /proc/cpuinfo の中をみてみるとサポートしているかわかります。

(FMAサポートしている:Core i5-4288Uでの結果)

f:id:curekoshimizu:20170806145146p:plain

話が脱線気味なので次の話題に移ります。

せっかくFMA命令をもっていても、 実際に使われていなければ意味がありません! プログラミング言語コンパイラの最近の状況はどうなのでしょうか?

FMAプログラミング言語

C言語C++言語だけの話に限定しますが、 C言語であれば C99FMA 計算用の関数をサポートしており、 C++言語であれば C++11 でサポートしています。

具体的には fmaf・fma・fmal という関数名で提供されています。 順に float・double・long double 用 になります。

C99 は名前の通り1999年の規格なわけですが、 FMA 演算

 \circ (x\times y + z)

を HW サポートしていたのは上でも書いたとおり最近の話になります。

実は HW サポートしていなくとも FMA 計算はできます!

ただしかなりの命令数を使うため HWサポートされていない場合は遅いです!!

あまりにも命令数がかかるため、ここでは「FMA計算」とあえて書きました。

ここに注意してください。

FMAコンパイラの最適化

例えば FMA を HW サポートしている場合、 コンパイラは自動的に FMA 命令に変更してくれたりするのでしょうか?

実はこのあたりが、さらにややこしい話をもたらします。

このあたり gcc と clang で話が少し違うのです。

コンパイラと最適化オプション

gcc も clang も最適化オプションをもっており、 デフォルトでは最適化を行わない設定になっています。

gcc・clnag の最適化オプションとして代表的なものは次です:

  • -O0 (default:最適化なし)
  • -O1
  • -O2
  • -O3
  • -Ofast

下にいくほどコンパイラにより最適化されます。

参考までに他にも次のような最適化オプションなどがあります:

  • -Og : Optimize debugging experience
  • -Os : Optimize for size

それぞれの最適化でどういったことをするかについて調査するには、 次を参考にするといいかもしれません:

(gcc)

(clang)

-march=native オプション

上の最適化オプションだけでは CPU による違いは気にしません。 具体的には

先ほどの Intel CPU の

  • 第1世代: Nehalem
  • 第2世代: Sandy Bridge
  • 第3世代: Ivy Bridge
  • 第4世代: Haswell (2013年頃)
  • 第5世代: Broadwell
  • 第6世代: Skylake
  • 第7世代: Kaby Lake (2017年現在頃)
  • 第8世代: Coffee Lake

といった違いは、 上の -O系のオプションだけをつけても吸収してくれません。

具体的には 上の最適化オプションだけをつけても、 FMA 命令をもつ Haswell と FMA 命令をもたない Sandy Bridge のような CPU のどちらでも 動作できるようにするため、

結論としては 「FMA命令をもたないコード」が生成される ことになります。

高精度・高速度な命令である FMA をもっていても、-O系のオプションをつけるだけでは使われない のです。

非常にもったいないです!

その最適化を支援するために、 この CPU だけでしか使わないよ? とか、 この演算ができるCPUでしか使わないよ? と伝えるオプションがあります。

それが例えば -march=native になります。

例えば gcc 5.4.0 を使って Core i7-3820 が乗った環境では

$ gcc -march=native

をつけるとと次のようなオプションが有効になります:

-march=sandybridge -mmmx -mno-3dnow -msse -msse2 -msse3 -mssse3 -mno-sse4a -mcx16 -msahf -mno-movbe -maes -mno-sha -mpclmul -mpopcnt -mno-abm -mno-lwp -mno-fma -mno-fma4 -mno-xop -mno-bmi -mno-bmi2 -mno-tbm -mavx -mno-avx2 -msse4.2 -msse4.1 -mno-lzcnt -mno-rtm -mno-hle -mno-rdrnd -mno-f16c -mno-fsgsbase -mno-rdseed -mno-prfchw -mno-adx -mfxsr -mxsave -mxsaveopt -mno-avx512f -mno-avx512er -mno-avx512cd -mno-avx512pf -mno-prefetchwt1 -mno-clflushopt -mno-xsavec -mno-xsaves -mno-avx512dq -mno-avx512bw -mno-avx512vl -mno-avx512ifma -mno-avx512vbmi -mno-clwb -mno-pcommit -mno-mwaitx --param l1-cache-size=32 --param l1-cache-line-size=64 --param l2-cache-size=10240 -mtune=sandybridge 

例えば キャッシュの大きさや SSE命令 (Streaming SIMD Extensions) のどれをサポートしているかなどを、 コンパイラに伝えていることがわかります。

ここでわかることは、

-mno-fma
-mno-fma4
-avx512ifma

というオプションがついており、FMA命令をもっていないんだろうな… ということが容易に想像つきます。

-march=native でどういったオプションが有効になるかについては次の記事が参考になります:

また、ここまで紹介した -O1/-O2/-O3/-Ofast や -march=native をつけるとどの程度速くなるかについては、 例えば次の記事が参考になるかもしれません:

http://www.phoronix.com/scan.php?page=news_item&px=GCC-Optimizations-E3V5-Levels

このあたりの最適化オプションを勉強するだけでも面白いです。

FMA が発行されるための最適化オプションは?

FMA 命令をもった Haswell 環境で次のコードをコンパイルしてみましょう。

float fma_test(float x, float y, float z)
{
    return x*y + z;
}
$ gcc -S hoge.c

コンパイルすると次のようになります:

        .....
    mulsd   -16(%rbp), xmm0
    addsd   -24(%rbp), %xmm0
        .....

つまりは FMA命令が使われず乗算と加算命令 で実行されることになります。

$ gcc -O3 -march=native

コンパイルするとどうでしょうか?

        .....
    vfmadd132sd     %xmm1, %xmm2, %xmm0
        .....

FMA 命令 (詳しくは上で説明したFMA3命令) が発行されていることがわかります。

一方で

$ clang -O3 -march=native

はどうでしょうか?

        .....
    vmulsd  %xmm1, %xmm0, %xmm0
    vaddsd  %xmm2, %xmm0, %xmm0
        .....

FMA 命令が使われません。

つまり、気をつけたいポイントは gcc・clang のオプションを揃えても FMA 化されるか違う場合がある ということです。

この違いを調査した表が次になります:

float fma_test(float x, float y, float z)
{
    return x*y + z;
}

上のコードは次のオプションで FMA命令 が呼ばれるかの表

(Core i5-4288U CPU (Haswell/ FMA をもつCPU))

オプション gcc 5.4.0 clang 3.8.0
-O0 × ×
-O1 × ×
-O2 × ×
-O3 × ×
-Ofast × ×
-march=native -O0 × ×
-march=native -O1 × ×
-march=native -O2 ×
-march=native -O3 ×
-march=native -Ofast

(Core i7-3820 (Haswell/ FMA をもたないCPU))

オプション gcc 5.4.0 clang 3.8.0
-O0 × ×
-O1 × ×
-O2 × ×
-O3 × ×
-Ofast × ×
-march=native -O0 × ×
-march=native -O1 × ×
-march=native -O2 × ×
-march=native -O3 × ×
-march=native -Ofast × ×

まず march=native の重要性がわかります。 これがないと FMA 命令をもたないかもしれませんので FMA 命令は発行されません。

また gcc と clang で FMA が有効になる最適化ポイントが違います。 gcc ならば -O2 から clang なら -Ofast から のようです。

-ffp-contract というオプションがあるのですが、 これをつけると clang でも -O1 から可能ならば fma 命令が発行されるようになります。

オプション gcc 5.4.0 clang 3.8.0
-ffp-contract=fast -march=native -O0 × ×
-ffp-contract=fast -march=native -O1 × ○(変更点)
-ffp-contract=fast -march=native -O2 ○(変更点)
-ffp-contract=fast -march=native -O3 ○(変更点)
-ffp-contract=fast -march=native -Ofast

この -ffp-contract は gcc のリファレンスにこのように記載されています:

-ffp-contract=off disables floating-point expression contraction. -ffp-contract=fast enables floating-point expression contraction such as forming of fused multiply-add operations if the target has native support for them. -ffp-contract=on enables floating-point expression contraction if allowed by the language standard. This is currently not implemented and treated equal to -ffp-contract=off.

The default is -ffp-contract=fast.

これをみると gcc は既に -ffp-contract オプションがついています。 一方で clang は挙動が変わったことからもデフォルトが fast になっていないようです。

このあたりに コンパイラの思想の違いがある ように見受けられます。

また、 gcc でも -O1 で最適化してもらうには -fexpensive-optimizations をつけると大丈夫です。

さらに、 -march=native は いくつかのオプションをつけてもらうためのものだったわけですが、 実際のところは -mfma が本質的です。

そのためまとめると

float fma_test(float x, float y, float z)
{
    return x*y + z;
}

というコードは

$ gcc   -march=native -O2
$ gcc   -march=native -O1 -fexpensive-optimizations

$ gcc   -mfma         -O2
$ gcc   -mfma         -O1 -fexpensive-optimizations

$ clang -march=native -Ofast
$ clang -march=native -O1 -ffp-contract=fast

$ clang -mfma         -Ofast
$ clang -mfma         -O1 -ffp-contract=fast

のいずれかで FMA 命令に展開されることがわかります。

#include <math.h>

float fma_test(float x, float y, float z)
{
    return fma(x, y, z);
}

としてみましょう。

この場合は

(Core i5-4288U CPU (Haswell/ FMA をもつCPU))

オプション gcc 5.4.0 clang 3.8.0
-O0 × ×
-O1 × ×
-O2 × ×
-O3 × ×
-Ofast × ×
-march=native -O0
-march=native -O1
-march=native -O2
-march=native -O3
-march=native -Ofast

(Core i7-3820 (Haswell/ FMA をもたないCPU))

オプション gcc 5.4.0 clang 3.8.0
-O0 × ×
…. …. ….
-march=native -Ofast × ×

となり、HW命令さえもっていれば FMA に展開されます。

そのため、どうしても FMA を発行させたい場合には有効な方法と言えます。

まとめますと

#include <math.h>

float fma_test(float x, float y, float z)
{
    return fma(x, y, z);
}

というコードは

$ gcc   -march=native
$ gcc   -mfma

$ clang -march=native
$ clang -mfma

で展開されます。こちらは -O0 でも大丈夫なところが大きな違いです。

また、 clang であっても -ffp-contract=fast も必要ありません。 それは既にユーザーが FMA 演算箇所を指定しているからです。

FMA の ハードウェアサポートがない場合に C99 FMA 関数 を呼ぶと?

この計算が

 \circ (x\times y + z)

きちんと発行できます!

ただし、問題は 複数命令で実行しているので遅い という点です。

それをみるべく次のことをみてみましょう。

HWサポートがない場合に 単精度浮動小数点用の C99 FMA 関数の fmaf を使った実行ファイルを nm コマンドにかけてみますと

$ nm a.out | grep fmaf
                 U fmaf@@GLIBC_2.2.5

となり glibc の fmaf 関数を呼んでいる ことがわかります。

具体的には glibc のコードの

sysdeps/ieee754/dbl-64/s_fmaf.c

が呼ばれており、ソフトウェアで実装されています。

この実装は次の論文が元になっています:

つまり、高速化という観点ではなく丸めの影響を抑えるため、 複数の計算で実現しているというわけです。

これについては丸め誤差の面白い技法が使われており、 詳しく紹介してみたいと思っております。

FMA 応用事例 – 多項式計算による FMA の利点

ここまでいろいろ知らないと、 FMA 命令の恩恵が得られないことがわかります。

そして、FMA 命令を使うことなんてそんなにあるの? という声があります。

非常にたくさんあります!

紹介する事例が多すぎるため、

ここでは 多項式計算 で考えてみたいと思います。

次の多項式を考えてみましょう

 2 {x}^{7} -7 {x}^{6} + 12{x}^{5} - 50{x}^{4} - 2{x}^{3} - 66 {x}^{2} - 60{x} + 54

こうした例えば 5次の方程式を計算することがあると思いますが、  a_{k} x^{k} という式を  k 回の乗算 で計算すると 0+1+2+3+4+5+6+7 = 28回の乗算と 7回の加減算で実行することになります。 つまり 35 回の丸めが発生します。

これを例えば次のような括弧の付け方で考えます

  2 \left( {x}^{7}  -7 \left( {x}^{6} + 12 \left( {x}^{5} -50 \left( {x}^{4} - 2 \left( {x}^{3} - 66 \left( {x}^{2} + ( -60 {x} + 54) \right) \right) \right) \right) \right) \right)

すると、  x^{k} という式を  (k-1) 回の乗算 で計算すると 0+0+1+2+3+4+5+6 = 21回の乗算と 7回のFMAで実行することになります。 つまり 28 回の丸めが発生します。

つまり 35 回の丸めが入る式から 28 回の丸めが入る式 に変わりました。

さらに、上の 5次方程式は次の Horner 法とよばれる形に変換してあげることができます:

 \left( \left( \left( \left( \left(  \left( 2{x} -7 \right){x} +12 \right){x} -50 \right) {x} -2 \right) {x} -66 \right) {x} -60 \right) {x} + 54

こうすると FMA 7 回で計算できてしまいます。

ちなみに Honer法 (ホーナー法) は非常に有名であり、 気になり起源論文も調べてみました。 1819年の Honer の論文になります:

もちろんこの Honer法 は有名なのですが、 高速化の観点にたてば、例えばこのような式変形も注目すべきです。

 ( {x}^{2} + 5)\left( ( {x}^{2} + 3) \left( ({x}^{2}-2)(2{x}-7)-8 \right) -9 \right) + 9

この方法にしますと 7回の FMA 演算で計算できてしまいます。

この特徴は 並列に計算可能 な 4つの FMA 演算が

  •  {x}^2 + 5
  •  {x}^2 + 3
  •  {x}^2 - 2
  •  2x-7

登場することであり、この点で Honer 法より有利に働きます。

こうした多項式計算技法はいろいろ知られており、どこかでまとめたいと思っています。

最後に

ここまで FMA についていろいろ紹介してみました。

FMA はかなり役に立つ命令だったり、 性能向上に役立つものですので、 是非ご活用ください!

また、もし家庭でPCの購入予算がおりない場合、

今のPCはFMAがサポートしてないけど最近のはFMAが載っているから、 だいたい2倍くらい性能いいらしいよ!っていえば 奥さんも説得できそうな気がします。

そんな時にも使ってみてください。

FMA丸め誤差の話とすごく関わるため、今後このブログでは頻繁に登場予定です!

よろしくお願いします。

こんなにも数式が出てこなかったのは、この記事が初めてかもしれません。

逆数の近似命令と精度補正について (その1)

小清水 (@curekoshimizu) です。

久しぶりの投稿になります。

長期にわたり転職活動をしており、 かなり投稿に時間が空いてしまいました。

今回の記事は

 \displaystyle \frac{1}{a} に関する内容です!

逆数の近似から精度を高めたい!!!

 \displaystyle \frac{1}{a} の近似値  {x}_{1} が与えられたときに

精度を高めたいということありませんでしょうか?

このモチベーションは色んな場面で登場します。

それはなぜか?

これは、 ハードウェア が 逆数の近似を行う命令 を持っていることが多く、 それを利用したいためです。

例えばこちらのブログでは

qiita.com

AVX-512 の vrcp28pd 命令 を使って  \displaystyle \frac{1}{a} の近似値 を得た際のお話が書かれています。

AVX-512 サポートしていれば vrcp28pd命令・vrcp14命令 などの逆数近似命令をサポートしています。

AVX-512 をサポートしているプロセッサはかなり限られていますが、 AVX をサポートしているだけでも rcpps という逆数近似命令を持っています。

また、もう今は死んでいると言っても過言ではない IA-64 にも frcpa 命令 (Floating-Point Reciprocal Approximation) という逆数近似命令があります。

ちなみに IA-64 は Hardware として普通の除算命令をもっておらず (IEEE-754準拠した除算の意)、 この近似逆数命令からうまく精度を高めて IEEE-754 に適合した 除算に変える必要がありました。

その点では IA-64 ユーザーからすれば、 モチベーションである「逆数の近似から精度を高めたい!」ということは普通に行われていました。

ここでは代表的な逆数近似命令を4つ挙げてみました。 あくまでも近似ではありますが、 きちんと誤差保証もついています。

その対応表は例えば次です:

要求 命令 最大相対誤差
IA-64 frcpa  2^{-8.886}
AVX rcpss  1.5 \times 2^{-12}
AVX-512 vrcp14pd  2^{-14}
AVX-512 vrcp28pd  2^{-28}

こうした話は CPU だけでなく、 NVIDIA GPU では fdivdef などで検索してみていただければと思います。

参考文献:

これらは ハードウェアの制約からきまった精度 であり、

IEEE754-2008 が規定するような精度は不要なものの、 もう少し精度を高めたい。そして それに保証をつけたい!

こういう時に役立つ記事を目指します。

精度を高める前に重要な事実

逆数の近似において

 e = 1 - a {x}_{1}

という式は重要な意味をもちます。

 e = 1 - a {x}_{1} \displaystyle \frac{1}{a} {x}_{1} の(符号付き)相対誤差を表す。

(証明)  \displaystyle \frac{1}{a} {x}_{1} の(符号付き)相対誤差を記し、そのまま計算すると、

(符号つき)相対誤差  \displaystyle = \frac{ \frac{1}{a}- {x}_{1}}{ \frac{1}{a} } = 1 - a {x}_{1} = e

となることから証明完了。


また、この  eFMA (Fused-Multiply-ADD)  (a \times b + c 型) の1命令 で計算できる点も 注目したい事実です。

FMA をご存知ない方はこちらの記事をご覧ください:

math-koshimizu.hatenablog.jp

精度を高めるアルゴリズム紹介

精度を高めるアルゴリズムはいくつか知られております。 そのうちの3つを紹介します。

精度向上アルゴリズムその1
 e = 1 - a {x}_{1}
 x_{2} = {x}_{1} + {x}_{1} e
FMA 2命令で得られる  x_{2} は相対誤差  e^2 となる

(証明)

(符号つき)相対誤差  \displaystyle = \frac{ \frac{1}{a}- {x}_{2}}{ \frac{1}{a} } = 1 - a ( {x}_{1} + {x}_{1} e ) = e - a{x}_{1}e = e^2
精度向上アルゴリズムその2
 e = 1 - a {x}_{1}
 x_{2} = {x}_{1} + {x}_{1} e
 x_{3} = {x}_{1} + {x}_{2} e
FMA 3命令で得られる  x_{3} は(符号つき)相対誤差  e^3 となる

(証明) 同様.

精度向上アルゴリズムその3
 e = 1 - a {x}_{1}
 x_{2} = {x}_{1} + {x}_{1} e, \quad E = e^{2}
 x_{4} = {x}_{2} + {x}_{2} E
FMA 3命令で得られる  x_{3} は相対誤差  e^4 となる。
ただし 2命令目は並列に実行すること。

(証明) 同様.

上の証明は証明になっているか?

一般的な書物には、上の内容だけが書かれているのですが、

証明としてこれでいいのでしょうか?

というのは、FMA であっても 丸め誤差が生じるからです。

実数世界であれば上の証明で正しいのですが、

2・3命令であっても丸め誤差が生じる演算が数回発生した結果、

得られた結果は  {e}^{2}, {e}^{3}, {e}^{4} などの結果と比べて近い結果になるのでしょうか。

本ブログではここについて踏み込みたいと思っています。

丸めの影響も含めた証明

ここで丸め関数として RN (x の 最近接偶数方向丸め) のみを取り上げることにします。

他の RU・RD・UZ についても同様に証明できますので割愛します。

また、丸めについては下の記事を参考ください:

math-koshimizu.hatenablog.jp


目的: 「精度向上アルゴリズムその1」を RN も含めて誤差評価を行う:
 e' = \mathrm{RN}(1 - a {x}_{1})
 x'_{2} = \mathrm{RN}({x}_{1} + {x}_{1} e')

これについて 2つの方針で考えていきたいと思います。

この評価は本などに書かれていなかったため、計算間違いをしている可能性もあり、 間違っている場合にはご指摘いただけると助かります。

前者が大変、後者が簡単な式変形になります。

アルゴリズムその1 – 丸め誤差評価気合編

気合で三角不等式でなんとか結果をだした評価式がこちらになります。

わかりやすいように、次のように記号を定義します:

  •  e = (1 - a {x}_{1})
  •  e' = \mathrm{RN}(e)
  •  y = ({x}_{1} + {x}_{1} e')
  •  y' = \mathrm{RN}(y) (= x'_{2})

(方針) y' の相対誤差を出す前に 先に  y の評価を試みる。

  •  y \displaystyle \frac{1}{a} の相対誤差

 = | 1 - a y |

 = | 1 - a ({x}_{1} + {x}_{1} e') |

 = |  e - a {x}_{1} e' |

 \leq | e - e' | + | e' - a {x}_{1} e'|

 =| e - e' | + | e' e |

 = | e - e' | + | e (e' - e) + {e}^{2}  |

 \leq |e| \left| \frac{e - e'}{e} \right| + {|e|}^{2} \left| \frac{e - e'}{e} \right|   + {|e|}^{2}

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

 \leq |e| {2}^{-p} + {|e|}^{2} {2}^{-p}   + {|e|}^{2}

と良い評価式を得る。

最後に

  •  y' \displaystyle \frac{1}{a} の相対誤差

 = | 1 - a y' |

 \leq | 1 - ay | + | ay - a y' |

 = | 1 - ay | + | a | | y - y' |

 = | 1 - ay | + |a y| \left| \frac{y - y'}{y} \right|

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

 = | 1 - ay | + |a  y | {2}^{-p}

 = | 1 - ay | + |a  ({x}_{1} + {x}_{1} e') | {2}^{-p}

 = | 1 - ay | + | (1-e)(1+e') | {2}^{-p}

 \leq | 1 - ay | + \left( | (1-e)(1+e') - (1-e)(1+e) | + |(1-e)(1+e)| \right) {2}^{-p}

 = | 1 - ay | + \left( |1-e||e'-e| + |(1-{e}^{2}| \right) {2}^{-p}

 \leq | 1 - ay | + \left( (1+|e|) |e| \left| \frac{e - e'}{e} \right| + 1 \right) {2}^{-p}

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

 \leq | 1 - ay | + \left( |e|{2}^{-p} + {|e|}^{2} {2}^{-p} + 1 \right) {2}^{-p}

この結果に 先ほど計算した  y の評価式を使って

 \leq  |e| {2}^{-p} + {|e|}^{2} {2}^{-p}   + {|e|}^{2}  + \left(  |e|{2}^{-p} + {|e|}^{2} {2}^{-p}   + 1 \right) {2}^{-p}

 =  {|e|}^{2} (1 +  {2}^{-p} + {2}^{-2p}) +  |e|  (  {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

という結果を得る。

この方針で得られた

丸め誤差を含めた アルゴリズムその1 の評価式は

   {|e|}^{2} (1 +  {2}^{-p} + {2}^{-2p}) +  |e|  (  {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

 {|e|}^{2} に比べて少し大きくなっています。

2命令しかないのにこんなに評価が長くなってしまいました…。

同じような話が アルゴリズムその2, その3 にも適用できるとは思うのですが…。正直やりたくないです…。


そこで次のような式変形を思いついてみました!

アルゴリズムその1 – 丸め誤差評価数式処理編

ここで重要な性質である次の性質を利用してみたいと思います。

Proposition 1:相対誤差に関する性質
2進p桁精度環境とする。このとき、
(丸め可能な)実数  x に対して ある  \epsilon_x  |\epsilon_x| \leq {2}^{-p} なるが存在し
 \mathrm{RN}(x) = (1 + \epsilon_x)x
と表せる.

(証明)

丸め誤差 RN の性質より ( 以前の記事Proposition2より)

 \displaystyle \left| \frac{  \mathrm{RN}(x) - x }{ x } \right| \leq {2}^{-p}

ここで  x > 0 とすると

 \displaystyle    -2^{-p} x \leq   \mathrm{RN}(x) - x \leq 2^{-p} x

なので

 \displaystyle    (1 -2^{-p}) x \leq   \mathrm{RN}(x) \leq (1 + 2^{-p}) x

となる。  x < 0 の場合も同様にすると

 \displaystyle    (1 +2^{-p}) x \leq   \mathrm{RN}(x) \leq (1 - 2^{-p}) x

となる。すなわち  \epsilon_x  |\epsilon_x| \leq {2}^{-p} とできる。


それでは評価に移ります。

再びわかりやすいように、変数を次のように定義します:

  •  e = (1 - a {x}_{1})
  •  e' = \mathrm{RN}(e)
  •  y = ({x}_{1} + {x}_{1} e')
  •  y' = \mathrm{RN}(y) (= x'_{2})

ここで上の 相対誤差に関する性質 を使い、

 \epsilon_y, \epsilon_e  |\epsilon_y| \leq {2}^{-p}, |\epsilon_e| \leq {2}^{-p} となるもので

 y' = (1 + \epsilon_y) y, e' = (1+\epsilon_e)e とできる。

  •  y' \displaystyle \frac{1}{a} の相対誤差

 = | 1 - a y' |

 = \left| 1 - a (1 + \epsilon_y) y \right|

 = \left| 1 - a (1 + \epsilon_y) ({x}_{1} + {x}_{1} e')\right|

 = \left| 1 - a (1 + \epsilon_y) \left({x}_{1} + {x}_{1} (1 + \epsilon_e) e \right) \right|

これをすべて展開して  e = (1 - a {x}_{1}) を使って  a {x}_{a} を削除すると

 = \left| {e}^{2} (1 + \epsilon_e + \epsilon_y + \epsilon_e \epsilon_y) + e (- \epsilon_e -\epsilon_e \epsilon_y ) - \epsilon_y \right|

ここで  |\epsilon_y| \leq {2}^{-p}, |\epsilon_e| \leq {2}^{-p} を利用して

 \leq  {|e|}^{2} (1 +  {2}^{-p} + {2}^{-p} + {2}^{-2p} ) +  |e|  ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

 = {|e|}^{2} (1 + {2}^{-p+1} + {2}^{-2p} ) +  |e|  ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

となる。

よって この方針で得られた

丸め誤差を含めた アルゴリズムその1 の評価式は

{|e|}^{2} (1 + {2}^{-p+1} + {2}^{-2p} ) +  |e|  ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

となりました。先ほどの評価より、式変形は簡単になりましたが、少しだけ結果が悪くなりました。

アルゴリズムその1 丸め誤差を含む評価まとめ

目的: 「精度向上アルゴリズムその1」を RN も含めて誤差評価を行う:
 e' = \mathrm{RN}(1 - a {x}_{1})
 x'_{2} = \mathrm{RN}({x}_{1} + {x}_{1} e')

ここで得られる  x'_{2} \frac{1}{a} に対する相対誤差は

  •   {|e|}^{2} (1 +  {2}^{-p} + {2}^{-2p}) +  |e|  (  {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

  • {|e|}^{2} (1 + {2}^{-p+1} + {2}^{-2p} ) +  |e|  ( {2}^{-p} + {2}^{-2p} ) + {2}^{-p}

の2種類が得られた。

ただし、  p は、この浮動小数点環境の精度であり、

 e は、もとの  {x}_{1}  \frac{1}{a} との相対誤差である。


2命令の誤差評価だけでもここまで大変 ということがわかった…。

前者の方針のほうが評価はよいのだが、

後者の方が式はあっさりしており、同じ方針で アルゴリズムその2・その3 についても評価を出すことができた。

しかしこれは長すぎるので、番外編で公開するか公開しないかもう少し時間がたってから考えることにする。

もし他のよい方法、よい近似式が得られた方は是非教えてください!

このブログはこういう、あまり誰もやりたがらない式評価に積極的に戦っていこうと思います。

次回考えたいこと

精度向上アルゴリズムを実施することで、

その丸め誤差の影響にも打ち勝ち、

その環境のFULL精度まで出すことができるのだろうか?

これについてはいくつか有名な事実があり、

上の精度向上アルゴリズム以外の事実を利用するのが一般的です。

これらについてまとまったら「逆数の近似命令と精度補正について (その2)」 を公開しようと思います。

たとえばこちらのブログの方は、精度向上アルゴリズムを3回かけても2回目と同じになりさらにずれが生じると述べており、 FULL 精度を出すのに失敗しているように見えます:

qiita.com

こうした内容にも助言できるよう目指したいと思います。

ちなみにですが

qiita.com

これら上の2つのQiitaブログは精度向上に

 {z}_{n+1} = {z}_{n} (2 - a {z}_{n})  という Newton 法の式を用いていますが、

このブログのアルゴリズムと本質的に一緒です。  {z}_{n} (2 - a {z}_{n}) = {z}_{n} + {z}_{n} (1 - a {z}_{n})

であり、  {z}_{n} (1 - a {z}_{n}) に このブログでは  e と名前をつけているだけです。

ただし、どちらのブログも3命令で実行していますが FMA があれば 2命令で実行でき、丸め誤差的にもそのほうがよいです。

この関係からもわかります通り、

今回の精度向上アルゴリズムは Newton法が元になっております。

この事実については以前次のブログでご紹介させていただきました(実数だけでなく整数版についても言及、数学的な証明付き。)

math-koshimizu.hatenablog.jp

興味がありましたらご覧ください。

それでは、この記事をここまで読まれた方で、本当に興味のある方は次回をお待ちください。

よろしくお願いします。

丸め誤差界の Hello World 的定理 -- Sterbenz の定理

小清水 (@curekoshimizu) です。

昨日は 数学カフェ に初めて参加させていただき、

機械学習について刺激を受けました!

connpass.com

その中で、

コイン投げは確率界の Hello World というフレーズを聞いて

なるほどなるほど!と感じました。

そして、

そういえば 丸め誤差界の Hello World はなんだろう?

と思い記事を書いてみました。

そもそも、 丸め誤差 は業界というには狭すぎて、

そもそも 丸め誤差に関する定理を聞く人も多いのでは?とも思います。

また、丸め誤差に関する本も少ないことなどから、

このセレクトがなるほど!となるかもわかりませんが、

私の主観で選択しました。

私的に思う、

丸め誤差界の Hello World はこちらになります!

「Sterbenzの定理」

この Sterbenz の定理、ゴールドバーグ先生の

what every computer scientist should know about floating-point arithmetic

の Theorem 11 にも載っていますので 皆さんもちろん知っていますよね!

とは言いつつも紹介します!

丸め誤差界の Hello World – Sterbenz の定理

Sterbenzの定理
 \beta浮動小数点環境とする。(2進でなくても構わない)
 x, y が有限値をとる浮動小数点数であり、
 \dfrac{y}{2} \leq x \leq 2y
を満たすならば、  x - y の結果を、丸め誤差なく正しく  \beta浮動小数点数として表現可能。

この定理は 有限値をとる浮動小数点数 という過程しかなく、

正規化数 という数以外にも、 非正規化数 というクラスにも成立することを言っています。

この点でかなり広いです。

この 正規化数・非正規化数 については、

後半の 「Prop. の証明 – 整数を用いた  \beta p 桁の浮動小数点数の表示」の章 に概略を記載しましたのでご参考ください。

どこに Hello World 性があるの?

浮動小数点の丸め誤差理論は、

誤差を抑えたり評価したりということをする理論だと思います。

その理論としては、

誤差なく計算できる! という 誤差のない究極性! があり、

なおかつ 証明が極めて簡単! というところに Hello World 性があると思います。

Sterbenz の定理の証明

(証明)

 \dfrac{y}{2} \leq x \leq 2y

であるから、  x, y は同符号であることがわかる。

対称性から  x , y \geq 0 で証明できればよい。

また、

 \dfrac{y}{2} \leq x \leq 2y

から

 \dfrac{x}{2} \leq y \leq 2x

と、  x, y の関係をひっくり返しても同じ不等式が成立する。

このため、  x \geq y で証明できれば  x \leq y も同様である。

以上をまとめると

 x \geq y \geq 0 で証明できれば充分であることがわかる。


この浮動小数点環境が

 \beta p 桁の浮動小数点環境 (指数部範囲:  e_{min}, e_{max}) とすると、

有限値をとる浮動小数点数  x, y \geq 0 はそれぞれ

 x = M_x \times \beta^{e_x - p + 1}
 y = M_y \times \beta^{e_y - p + 1}

と表せる。ただし  M_x, M_y は整数で  0 \leq M_x, M_y <  \beta^{p} を満たし、

 e_x, e_y は整数で  e_{min} \leq e_x, e_y \leq e_{max} である。 (*)(この事実について後で補足あり)


 M := M_x \beta^{e_x - e_y} - M_y とおくと

 x - y = M \times \beta^{e_y - p + 1}

となり、 M は整数なので  0 \leq M <  \beta^{p} が示されれば、

 x - y が この浮動小数点環境で正確に表示できていることがわかる。 (*)(この事実について後で補足あり)

そのため、このことを示す。


仮定より  x \geq y なので

 x - y = M \times \beta^{e_y - p + 1} \geq 0

であり  M \geq 0 がわかる。

最後に

仮定より  2y \geq x であることから

 x - y = M \times \beta^{e_y - p + 1} \leq y \quad (= M_y \times \beta^{e_y - p + 1} )

であるから

  M \leq M_y <  \beta^p

が示され  0 \leq M <  {\beta}^{p} が証明された。

故に証明完了。

上の証明への補足

上の証明で唯一難しい点として、下の事実を当然として利用しました。

この事実は頭の体操的なものであるものの、

おそらくこのブログの多くの箇所で参照する事実であることから、

解説を記載しておくこととします。

Prop. (整数を用いた  \beta p 桁の浮動小数点数の表示)
浮動小数点環境が  \beta p 桁の浮動小数点環境 (指数部範囲:  e_{min}, e_{max}) とすると、
有限値をとる浮動小数点数  x \geq 0
 x = M \times \beta^{e - p + 1}
と表せる。ただし  M は整数で  0 \leq M <  \beta^{p} を満たし、
 e は整数で  e\_{min} \leq e \leq e\_{max} である。
そして逆に、 その  x が上の表示であれば、 有限な浮動小数点数  x \geq 0 である。

Prop. の証明 – 整数を用いた  \beta p 桁の浮動小数点数の表示

 \beta p 桁の 有限値をとる浮動小数点数  x \geq 0

 x = (x_0 . x_1 x_2 \dotsb x_{p-1} )_{\beta} \times \beta^{e} = \left( x_0 + \dfrac{x_1}{\beta} + \dfrac{x_2}{\beta^2} + \dotsb + \dfrac{x_{p-1}}{\beta^{p-1}} \right) \times \beta^{e}

と表される。ただし  x_i は整数で  0 \leq x_i <  \beta を満たし、

 e は整数で  e_{min} \leq e \leq e_{max} である。

これは浮動小数点数の定義からである。

補足すると  x_0 \neq 0 とできる数が特に 正規化数 と呼ばれるものであった。

この数にとらわれない、例えば

 (0 . 1 1 \dotsb 1 )_{\beta} \times \beta^{e_{min}} = \left( \dfrac{1}{\beta} + \dfrac{1}{\beta^2} + \dotsb + \dfrac{1}{\beta^{p-1}} \right) \times \beta^{e_{min}}

という数は 非正規化数 と呼ばれる数である。この数は  x_0 \neq 0 の形で表示することができない。これは指数部が  e_{min} と既に最小であるからである。

上の  x は次のように変形できる:

 x = \left( x_0 \beta^{p-1} + x_1 \beta^{p-2} + x_2\beta^{p-3} + x_{p-1} \beta^{0} \right) \times \beta^{e-p+1}

ここで

 M := x_0 \beta^{p-1} + x_1 \beta^{p-2} + x_2\beta^{p-3} + x_{p-1} \beta^{0}

とおけばよく  M は 整数で  0 \leq M <  \beta^{p} を満たす。

逆に  M が整数で  0 \leq M <  \beta^{p} を満たすならば

 M = x_0 \beta^{p-1} + x_1 \beta^{p-2} + x_2\beta^{p-3} + x_{p-1} \beta^{0}

と表示できることから、逆も示される。

まとめ

丸め誤差界の Hello World 的定理である Sterbenz の定理 を紹介しました。

簡単な証明、簡単な条件なだけあり、

そこまで日の目を見ないかもしれませんが

そこは Hello World 、目を瞑りましょう。

浮動小数点や丸め誤差などに興味がわきましたら

math-koshimizu.hatenablog.jp

などを是非御覧ください。以上です。

エイプリルフールだし「1=2」を別の方法で証明してみた

小清水(@curekoshimizu) です。

久々のブログ投稿になります!

今日は ロマンティック数学ナイトボーイズ がありました!

f:id:curekoshimizu:20170401211158j:plain

そこで 2進数フレンズ という内容で発表させていただきまして、

その発表資料はこちらです:

今日はこの発表資料をつかって、

「1 = 2」を証明しようと思います!

「1 = 2」とは?

アンサイクロペディアでは有名な記事で、

おかしなことをやって 1 = 2 を証明するという内容になります。

http://ja.uncyclopedia.info/wiki/1%3D2

もっちょさんがこれらを分類などされています!

motcho.hateblo.jp

今回紹介する 「1 = 2」証明方法とは?

このもっちょさんの分類に属さない 新しい証明方法です!

発表スライドの 30ページ目 のこちらを御覧ください:

f:id:curekoshimizu:20170401212007p:plain

こちらを引用すると Excel の計算により

 {2}^{50} + 1 - {2}^{50} = 0

となったことがわかります。

また、Excel の計算により

 ({2}^{50} + 1 - {2}^{50}) = 1

ということもわかります。

これら 2式より

 (  0  ) = 1

となり、括弧を外して

 0 = 1

両辺に 1 を加えて

 1 = 2

が示された。

以上です!

エイプリルフールらしい記事になりました!

(しかし Excel の挙動はエイプリルフールではないことに注意。)

おまけ

今年の1月に 第8回日曜数学会 というところで

級数を大改造!! 劇的ビフォーアフター という内容で発表しましたので

その動画も掲載しておきます!

www.nicovideo.jp

こちらの発表の方が面白いかもしれません。

計算環境の精度を当てる ― 解説編 ―

小清水 です。 (@curekoshimizu)

前回の記事は

「計算環境の精度を当てる ― Intel CPU と Excel への応用 ―」というタイトルで、

 \beta p 浮動小数点環境

基数  \beta精度  p を当てる!

という話でした。 (一部 Excel の謎挙動を紹介)

math-koshimizu.hatenablog.jp

前回は証明を掲載していなかったため、

今回は証明に主眼をおいて紹介したいと思います。

最初に注意しますが、

 \beta は 2以上の整数 とします。

世の中には 「-2 進数」・「1+i進数」・「黄金進数」などというのがあったりもするのですが、

さすがに除外させてください。

(これらが採用された実用的なコンピューターを僕はまだ知りません。)

証明がわかりやすいよう、

記号  \circ (x) を次で定義しておきます:

 \circ (x) の定義
 \circ (x) を x について  \beta p桁環境 で計算した際の結果
と定義します。
ただし丸めの方向は環境に決っているものとし、
最もよく使われる「最近接偶数方向丸め」でなくとも構いません。

丸めの方向 、例えば、

最もよく使われる 最近接偶数方向丸め0方向丸め などについて知りたい方はココを御覧ください:

math-koshimizu.hatenablog.jp

3進数3桁の環境でアルゴリズムを振り返りながら+証明

実例があったほうがわかりやすいかと思いますので、 ありきたりな 2進数 ではなく 3進3桁 で考えましょう。

以前紹介していたように、

Step1・2・3 がありました。

Step のそれぞれのアルゴリズムを思い出し、

次にそれを 3進3桁で実例を動かし、

そのアルゴリズムで結局何を実行しようとしているのかとその証明を述べていこうと思います。

Step1. について

Step 1
 x_0 = \circ{(1)} = 1 から始め、
 x_{m+1} = \circ{(2\times x_m)} と計算していく。
つまり  {2}^{n}
 \beta p 桁で計算していくということになります。
そして、初めて  \circ{( \circ (x_m + 1) - x_m)} \neq 1 となる  x_m を求めます。

実例でみていきましょう。

Step1 を実例で考える

 {2}^{0} = 1 の 3進3桁表現は  {(1.00)}_{3} \times {3}^{0} です。

 {2}^{1} = 2 の 3進3桁表現は  {(2.00)}_{3} \times {3}^{0} です。

これを続けていくと、

f:id:curekoshimizu:20170129201119p:plain

となります。注意としては  {(1.21)}_{3} \times {3}^{2} の次 ( \times 2) は

本来は  {(1.012)}_{3} \times {3}^{3} なのですが、

これは 4桁精度 の表現なので、

3桁精度 に丸める必要があります。

このアルゴリズム丸めの方向に言及していない ことから、

近い2つの浮動小数点のどちらを選んでも問題ありません!

つまり候補としては 2つ あり、この例ですと

 \circ ({(1.012)}_{3} \times {3}^{3}) = {(1.01)}_{3} \times {3}^{3} または  {(1.02)}_{3} \times {3}^{3}

となります。

f:id:curekoshimizu:20170130000236p:plain


ここまで登場した、オレンジ以外の数 を  x としても

 (x + 1) - x の計算値は 1 になります。

例えば、ピンクの数である、

 x = {(1.21)}_{3} \times {3}^{2} で詳しくみてみましょう。

 x + 1 = {(1.21)} \times {3}^{2} + {(1.00)} \times {3}^{0}

ここで 桁あわせを行い

  = {(1.21)} \times {3}^{2} + {(0.01)} \times {3}^{2}

  = {(1.22)} \times {3}^{2}

と途中で丸めなのも必要なく、

 x+1 を正確に表すことができています。

しかしながら、オレンジの 指数3 の数はこうはなりません。

例えば  {(1.12)}_{3} \times {3}^{3} を選択した場合を例にあげると

f:id:curekoshimizu:20170130000335p:plain

ここでわかるのは 1 の効果が丸めによって(☆)、復元不能になっていることです。

終結果は丸めの方向によって、

  •  {(1.20)}_{3} \times {3}^{3}
  •  {(1.12)}_{3} \times {3}^{3}

の2つがでていますが、もとの  x = {(1.12)}_{3} \times {3}^{3} から見ると

差が 3 または 0 となっています。

すなわち、 1 という 指数0 の項の影響が、 指数1 の世界に持ち上がってしまったというイメージです。


 (x + 1) が正しく求められれば、  (x + 1) - x は当然 1になるのですが、

 (x + 1) が真の解からずれてしまうと、もはや  (x + 1) - x で 1 に復元することはできなくなります。

すなわち、 Step 1 からわかることは

3進3桁では  x_5 = {(1.12)}_{3} \times {3}^{3} で計算が合わなくなるということです。

ここからわかるアルゴリズムの目的と証明

Step1 のアルゴリズム

Step 1
 x_0 = \circ{(1)} = 1 から始め、
 x_{m+1} = \circ{(2\times x_m)} と計算していく。
つまり  {2}^{n}
 \beta p 桁で計算していくということになります。
そして、初めて  \circ{( \circ (x_m + 1) - x_m)} \neq 1 となる  x_m を求めます。

は何を求めるための Step になるのでしょうか。

 \beta p 桁表現の環境における、指数 p表現に属する数  x_m を求めるためのアルゴリズム

実はこれが主張になります。

先程の例でいうと、

f:id:curekoshimizu:20170130000335p:plain

  • 水色 (指数0型)
  • 紫色 (指数1型)
  • ピンク色 (指数2型)

を経て、

  • オレンジ色 (指数3型)

でちょうど Step1 が停止したのです。

この点を一般の場合でも証明しましょう。


Prop 1
整数  m m <  {\beta}^{p} を満たすならば、  m \beta p 桁で丸め誤差なく表現できる

(証明)

 m <  {\beta}^{p} なので  0 \leq x_m <  \beta なる自然数を使って

\begin{align*} m = x_{p-1} {\beta}^{p-1} + x_{p-2} {\beta}^{p-2} + \dotsb + x_{0} {\beta}^{0} \end{align*}

と表現できる。

ここで  x_i \neq 0 となる最大の添字を  k とすれば

\begin{align*} m = x_{k} {\beta}^{k} + x_{k- 1} {\beta}^{k- 1} + \dotsb + x_{0} {\beta}^{0} = (x_{k} . x_{k- 1} \dotsb x_0) \times {\beta}^{k} \end{align*}


この Prop 1 から  {\beta}^{p} を超えない限り  {2}^{n} \beta p桁 で丸め誤差なく表現できる。

 {2}^{n} <  {\beta}^{p} であるから  {2}^{n} + 1 \leq {\beta}^{p} である。

 {2}^{n} + 1 <  {\beta}^{p} の場合は Prop1 から  \beta p桁 で表現可能。

そうでない場合は  {2}^{n} + 1 = {\beta}^{p} であり

 {\beta}^{p}   = {(1.00\dotsb 0)} \times {\beta}^{p} と表されるので、 \beta p桁 で表現可能。

よって

 {2}^{n} <  {\beta}^{p} を満たす限り  (x + 1)  \beta p桁 で正確に表現可能であるから、

 \circ (x +1) = x + 1

となり、

ここから  \circ (\circ (x + 1) -x) = \circ( (x +1) -x ) = \circ (1) = 1

ということがわかった。つまり、

Prop 2
Step1 で求める  x_n x_n < {\beta}^{p} を満たす限り  x_n = {2}^{n} となり、さらに   \circ (\circ (x + 1) -x)  = 1 を満たす

ということがわかったことになる。


また、

Prop 3
 {\beta}^{p} \leq {2}^{n} <  {\beta}^{p+1} を満たす最小の自然数  N がとれる

(証明)

この式を変形すると  p\log_2{\beta} \leq {n} <  (p+1)\log_2{\beta} であり

この区間の端点間の差は

 (p+1)\log_2{\beta} - p\log_2{\beta} = \log_2{\beta} \geq 1

と 1以上のため、区間  \Big[p\log_2{\beta} , (p+1)\log_2{\beta} \Big) には必ず自然数が含まれる。

その自然数が存在し、その最小なものが証明する  N である。


この Prop3 より 自然数  m として

 {\beta}^{p} \leq {2}^{m} <  {\beta}^{p+1} を満たし  {2}^{m-1} <  {\beta}^{p} となるものがとれる。

この式から

 \circ ({\beta}^{p}) \leq \circ (2\times {2}^{m-1}) <  \circ ( {\beta}^{p+1}) となり、

 {2}^{m-1} <  {\beta}^{p} であることから Prop 2から  x_{m-1} = {2}^{m-1} となる。

すなわち  x_m =  \circ (2\times x_{m-1}) = \circ ( 2 \times {2}^{m-1} ) であり

 {\beta}^{p} \leq x_m <  {\beta}^{p+1} を得る。

 \circ (x_m + 1) = x_m + \beta , x_m である(両者は  {\beta} p 桁で表現可能な  x_N + 1 の両端点である。

以上から

Prop 4
 {\beta}^{p} \leq x_{m} <  {\beta}^{p+1} となり  x_{m-1} <  {\beta}^{p} となる 自然数  m が存在し、  \circ ( \circ (x_m + 1) - x_m = \beta , 0 を満たす。

ということが示された。

以上から Step1 のアルゴリズム

 \beta p 桁表現の環境における、指数 p表現に属する数  x_m を求めるためのアルゴリズム

となることは Prop2 と Prop4 から示されたことになります。

Step2. について

Step 2
Step1 で求めた  {x}_{m} \circ ( \circ ( x_m + 1) - x_m) \neq 1
であり、この 1 を動かして
 \circ ( \circ ( x_m + 2) - x_m) \neq 2
 \circ ( \circ ( x_m + 3) - x_m) \neq 3
と続き、 初めて等号となってしまう 整数  \beta を見つる
つまり、
 \circ ( \circ ( x_m + \beta) - x_m) = \beta
となる  \beta を見つける。

この Step2 から  \beta が見つかると

その環境が  \beta 進数計算環境であることがわかる

Step2 を実例で考える

先程の 3進3桁 では Step1 により  x_5 = {(1.12)}_{3} \times {3}^{3} が求まりました。

 \circ ( \circ ({(1.12)} \times {3}^{3} + 1) - {(1.12)} \times {3}^{3}) が 1 にならないことがわかりました。

今度は

 \circ ( \circ ({(1.12)} \times {3}^{3} + 2) - {(1.12)} \times {3}^{3})

 \circ ( \circ ({(1.12)} \times {3}^{3} + 3) - {(1.12)} \times {3}^{3})

を計算していこうという話になります。

計算してみるとわかるのですが下のようになります:

  •  \circ ( \circ ({(1.12)} \times {3}^{3} + 1) - {(1.12)} \times {3}^{3}) \neq 1 (Step1)

  •  \circ ( \circ ({(1.12)} \times {3}^{3} + 2) - {(1.12)} \times {3}^{3}) \neq 2

  •  \circ ( \circ ({(1.12)} \times {3}^{3} + 3) - {(1.12)} \times {3}^{3}) = 3

これは

3 が 3進3桁で  {(1.000)}_{3} \times {3}^{1}指数1 に初めてなる!

ことが関係します。それまでの数は 指数0 型の表現でした。

さきほど成立しなかった理由を思い出すと、

指数0 の数を足しても、桁数をあわせて足しても 3桁でおさまらないため、

丸める必要がありました。

今回は 指数1 なので、桁をあわせて足した数が 3桁でおさまります。

(x+1) がきちんと表現できてしまえばこちらのものです。

以上から、 3 が求められ、

3進3桁の 3進数 が求められます。

ここからわかるアルゴリズムの目的と証明

目的は上でも記載したとおり、基数  \beta を求めることです。

Step1. で  x_m であり

 {\beta}^{p} \leq x_m <  {\beta}^{p+1}

となるものが得られました。

この数と  1, 2, \dotsb , \beta - 1 の足し算は  \beta p 桁で正確に表現できないが、

 \beta  (1.000\dotsb ) \times {\beta}^{1} であるから  x_m + \beta は正確に表現可能。

よって  \circ ( \circ ( x_m + \beta) - x_m ) = \circ ( x_m + \beta - x_m) = \beta となり Step2 で  \beta が求まる。

Step3. について

Step3

Step 3
Step2 でわかった  \beta を使って  y_0 = \circ{(1)} = 1 から始め、
 y_{n+1} = \circ{(\beta \times x_n)} と計算していく。
そして、初めて  \circ ( \circ ( y_p + 1) - y_p) \neq 1 となる  p を求めます。

この  p精度  p を表すという主張です。

Step3 を実例で考える + 証明

 27 = {(1.000)}_{3} \times {3}^{3} のとき、

初めて  (27 + 1) - 27 の計算が 1 に一致しなくなります。

f:id:curekoshimizu:20170129212915p:plain

これは Step1 と同じ話で、

何が違うかというと、

Step1 は基数が不明だったので、  {2}^{p} という形でもって調べる必要がありましたが、

今度は  {\beta}^{p} 型なので、第  p 回で 指数  p 型になることがわかります。

これが証明であります。

あとは Step1 と同じです。

最後に: \beta素数じゃない場合でも?

基数  \beta素数じゃなくとも、きちんと求まります!

これは Step2 の証明を思い出すとわかります。

ということは

16進法 でもこの基数を当てられる!ということになります。

16 進法といえば 古典的 RISC アーキテクチャである (1970年頃)、

IBM System/370 などに採用されている IBM浮動小数点方式 では

「16進数」を採用していました。

「16進数」は一時期はコンピューターの歴史を彩っていた方式であり、

これを使う機会があれば、

「16進数」である証拠をこのアルゴリズムで確認できます!


以上1万字を超える、浮動小数点を含んだ証明の話でした。

前回の「計算環境の精度を当てる ― Intel CPU と Excel への応用 ―」という記事では、

math-koshimizu.hatenablog.jp

Excel の謎挙動 を紹介しました。これについては現在も調査中です。

大まかにわかってきたのですが、まとまりましたら紹介したいと思っております!

よろしくお願いします。

計算環境の精度を当てる ― Intel CPU と Excel への応用 ―

小清水 (@curekoshimizu) です。

計算ツールに対して

その 計算精度 がわかると便利ではないでしょうか?

例えば、コンピュータにのっている

その Excel はどのくらいの精度 で計算してくれる?

こうしたアプリケーションだけでなく、

そのコンピュータの中に入っている、

Intel の CPU はどのくらいの精度 で計算してくれる?

こういうことがわかると 非常に便利・ワクワク します!

特に CPU にはいろんな精度で計算できる命令が揃っており、

Intel x86_64 CPU ですと

  • 単精度
  • 倍精度
  • 拡張倍精度

の命令が載っています。

この記事を読むと、

実感をもって これらの精度がわかってしまいます。

f:id:curekoshimizu:20170106231338p:plain

しかも 数学的背景 もあるので納得もできます!

Abstract

1. CPU 単精度・倍精度 の 計算精度を当てる

2. CPU 拡張倍精度 の計算精度を当てる

3. Excel の計算精度を当てる!?

こういう話のとき Excel の例を最初に書きそうなものです。

なぜだか Excel の例が 最後 になっていますが

それは 困った話 が見つかってしまったからです。

1. CPU 単精度・倍精度 の 計算精度を当てる

プログラミングを勉強すると必ず習う計算環境 である、

単精度・倍精度浮動小数点環境の計算精度をまとめてみましょう。

単精度の浮動小数点 (float・binary32) 環境は 2進24桁精度

倍精度の浮動小数点 (double・binary64) 環境は 2進53桁精度

と言われています。これを実証してみましょう

精度を確かめるプログラムの方針

((x + 1) - x) と 1 が同じ

となるのが、

普通の実数の世界 であります。

当たり前です。


しかしながら、

普通じゃないのが 浮動小数点数 であり、

この2つの式が成立しない場合があります。

Step1

 {x}_{p} = {2}^{p} の形で

 p = 0 から初めて大きくしていき

初めて成立しなくなる  p を見つけます

Step2

その  p {x}_{p} = {2}^{p} に対して

 (( x_p + 1) - x_p) \neq 1

が成立しています。

この 1 を動かして

 (( x_p + 2) - x_p) \neq 2

 (( x_p + 3) - x_p) \neq 3

と続き、

初めて等号となってしまう 整数  \beta を見つけます

つまり、

 (( x_p + \beta) - x_p) = \beta

となる  \beta です。

この  \beta がわかると

 \beta 進数計算環境であることがわかる

Step3

Step2 でわかった  \beta を使って

もう一度 Step1 に似たアルゴリズムを実施します。

 {y}_{p} = {\beta}^{p} の形で

 p = 0 から初めて大きくしていき

初めて成立しなくなる  p を見つけます

つまり  (( y_p + 1) - y_p) \neq 1 となる

初めての  p精度  p を表す

という主張です。

補足的には、

2進数であることがわかっていればStep1とStep3 は等価になります。

このことについて

証明は次回の記事とさせていただき、計算例を見てみましょう

実際のプログラム

#include <stdio.h>

void check_float_precision(void)
{
    /* Step. 1 */
    float xp;
    for (xp = 1.0f; (xp + 1.0f) - xp == 1.0f; xp *= 2.0f) ;

    /* Step. 2 */
    float beta;
    for (beta = 1.0f ; (xp + beta) - xp != beta; beta += 1.0f) ;
    printf("base : %d\n", (int)beta);

    /* Step. 3 */
    int precision = 0;
    for (xp = 1.0f; (xp + 1.0f) - xp == 1.0f; xp *= beta) {
        precision++;
    }
    printf("precision : %d\n", precision);
}

void check_double_precision(void)
{
    /* Step. 1 */
    double xp;
    for (xp = 1.0; (xp + 1.0) - xp == 1.0f; xp *= 2.0) ;

    /* Step. 2 */
    double beta;
    for (beta = 1.0 ; (xp + beta) - xp != beta; beta += 1.0) ;
    printf("base : %d\n", (int)beta);

    /* Step. 3 */
    int precision = 0;
    for (xp = 1.0; (xp + 1.0) - xp == 1.0; xp *= beta) {
        precision++;
    }
    printf("precision : %d\n", precision);
}

int main(void)
{
    check_float_precision();
    check_double_precision();

    return 0;
}

実行結果

base : 2
precision : 24
base : 2
precision : 53

と先程紹介したとおり、

2進24桁・2進53桁精度 の 「2・24・2・53」が出力されております。

さらに次の例をみてみましょう。

2. CPU 拡張倍精度 の計算精度を当てる

CPU は過去互換性を大事にしており、

Intel x86_64 CPU には すごく昔によく用いられていた、

x87 命令といわれている、

拡張倍精度命令 というものをもっております。

拡張倍精度環境は 2進64桁精度 と言われています。

こうした 古めかしい命令 もわかってしまうのでしょうか?

gcc を使ったコンパイルで -mfpmath=387 オプションをつけると、

強制してこの x87命令 で計算することができます。

$ gcc -mfpmath=387 precision_check.c -o precision_check
$ ./precision_check

実行結果

base : 2
precision : 64
base : 2
precision : 64

float・double と宣言したプログラムであっても、

それが 無視されて 先程紹介したとおり、

2進64桁精度 の 「64」が出力されております。

3. Excel の計算精度を当てる!?

先程のアルゴリズムExcel でも計算してみましょう。

B列に  x_p = {2}^{p} を計算し、

その結果を使い  (x_p + n) - x_p を計算した表を作成します。

この表を使えば

f:id:curekoshimizu:20170122141737p:plain

となり、

f:id:curekoshimizu:20170122141743p:plain

このことから

Excelは 2進50桁 の計算環境 であることがわかります。

と思いますよね?

2進50桁ってなんだ?ということになるわけです。

なぜならあまり聞き慣れない精度のためです。

よく使われる精度は上でこれまででてきた

  • 2進24桁 (単精度)
  • 2進53桁 (倍精度)
  • 2進64桁 (拡張倍精度)

であります。

ためしに、上の表を次のように変えてみます。

 (x_p + n) - x_p を計算するのではなく

 ((x_p + n) - x_p) - n を計算してみます。

そして、正しく計算されたところは 0 になります。

f:id:curekoshimizu:20170122142854p:plain

もっと下まで計算しますと

f:id:curekoshimizu:20170122142859p:plain

となり、

同じく2進数であることは確かめられるのですが、

Excelは 2進53桁 の計算環境 になりました。

何がいいたい?

ある方法だと Excel は「2進50桁」

ある方法だと Excel は 「2進53桁」(倍精度) という結果になりました。

どうしてこうなったのかわからないため、

Excel詳しい方教えてください!

ちなみにですが CPU の倍精度環境で 同じ計算をすると、

 ((x_p + n) - x_p) - n の表と同じ結果を得ますので、

Excel はおそらく 倍精度精度で計算できるのだと思いますが…。

ということはアルゴリズムがおかしい?

このアルゴリズム浮動小数点環境の精度を求める方法として

よく知られたはずの方法なのですが

端的にいうと Excel の次のおかしな結果のせいで求められないことがわかります:

f:id:curekoshimizu:20170122150814p:plain

まとめ

CPU を例に 浮動小数点環境の精度をあてるプログラムを提示できました。

残念ながら Excel の場合はよくわからない結果になってしまった…。

(Excel が変な挙動を示す新しい例?)

記事が長くなってしまったため、

次回の記事はこのアルゴリズムの証明を紹介予定です。

浮動小数点数については下で詳しくまとめていますので

math-koshimizu.hatenablog.jp

こちらも是非御覧ください!


(追記)

本記事の証明記事書きました!

math-koshimizu.hatenablog.jp

expm1 や pow などの 指数・対数函数 への考察

小清水 (@curekoshimizu) です。

今日は x y 乗、

 {x}^{y} について考えてみたいと思います。

この定義ってどのように習いましたでしょうか?

 {e}^{x} \mathrm{log}(x) を定義した後に

 {x}^{y} := {e}^{y \mathrm{log}(x)} により定義していたかと思います。

この記事ではいったん複素数のことは忘れましょう。

多価関数が...とか言い出すと本質からずれますので。

この定義を考えれば

 {x}^{y} をコンピューターで計算した値

 {e}^{y \mathrm{log}(x)} をコンピューターで計算した値

同じ計算結果になるということなのでしょうか。

こうしたことを調べるために、

そもそも  {e}^{x} \mathrm{log}(x) などを

どのようにコンピューターで計算するのか

ということを考える必要があります。

ここではプログラミング言語

ライブラリとして標準的に提供している函数 を使い

計算していみることにしましょう。

そのため、まずはライブラリで提供されている函数を調べることにします。

Outline

1. C言語 の標準ライブラリで用意されている 指数・対数函数まとめ

2. 他のプログラミング言語 と比較

3.  {x}^{y} {e}^{y \mathrm{log}(x)} の計算結果比較

1. C言語 の標準ライブラリで用意されている 指数・対数函数まとめ

math.h で公開されている 指数函数対数函数をまとめてみましょう。

指数函数シリーズ

Cの標準ライブラリには 指数函数 を計算する標準函数シリーズとして

  •  {2}^{x} を計算する exp2

  •  {e}^{x} を計算する exp

  •  {e}^{x} - 1 を計算する expm1

  •  {x}^{y} を計算する pow

があります。

対数函数シリーズ

Cの標準ライブラリには 対数函数 を計算する標準函数シリーズとして

  •  \mathrm{log}_2(x) を計算する log2

  •  \mathrm{log}_{10}(x) を計算する log10

  •  \mathrm{log}_{e}(x) を計算する log

  •  \mathrm{log}_{e}(1+x) を計算する log1p

考察してみる

Taylor展開結果 をご存知の方は

 {e}^{x} だけでなく

 {e}^{x} - 1 が用意されていることや

 \mathrm{log}_{e} (x) だけでなく

 \mathrm{log}_{e}(1+x) が用意されている

理由がわかるかと思います。

これは使い分けなくてはいけないのですか?

 x が非常に小さいときは 当然 です!

#include <stdio.h>
#include <math.h>

int main(void)
{

    const double x = exp2(-100);

    printf("%e\n", exp(x) - 1.0);
    printf("%e\n", expm1(x));

    return 0;
}

の実行結果は

0.000000e+00
7.888609e-31

となり  {2}^{-100} の exp(x) - 1 と expm1(x) の計算結果は異なることがわかります。

また、興味深いとおもったことは

対数函数 には  \mathrm{log}_{10}(x) があるにも関わらず、

指数函数 には  {10}^{x} が提供されていない 点です。

これは私の調査ミスでしょうか?理由をご存じの方いらっしゃいましたら教えてください。

もちろん  {x}^{y} を使えば  {10}^{y} が計算できますが、

それではなぜ  {2}^{x} 系のみ存在するのでしょう。

もう一点気になるのは  {x}^{y} が提供されているにも関わらず

 \mathrm{log}_{x}(y) が提供されていないのか ということです。

2. 他のプログラミング言語と比較

上では C言語について調査してみたわけですが

他の言語ではどのようになっているのでしょうか。

例として pythonrubygolang について調べてみました。

注意として

例えば  {10}^{x} 函数が標準ライブラリで提供されていない場合、

それは  {x}^{y} で代用できるわけですが 代用可能(△) ということにしました。

特に rubyx ** y という記法で  {x}^{y} を表しますが、標準ライブラリには函数が存在しないので (○) としました。

p 2.72 ** 3.14

python もこの記法が使えますが pow 函数が提供されていましたので ○ としています。

また golang には Pow10 という函数があるのですが、

整数値限定であり実数用ではないため除外しています。

C python ruby golang
 {2}^{x}
 {10}^{x}
 {e}^{x}
 {e}^{x} - 1 ×
 {x}^{y} (○)
 \mathrm{log}_2(x)
 \mathrm{log}_{10}(x)
 \mathrm{log}_{e}(x)
 \mathrm{log}_{e}(1+x) ×
 \mathrm{log}_{x}(y) × ×

このように差があると思っていませんでしたので、

間違いだよというのがありましたらご指摘ください。

調査につかったページリンク

ここからわかる面白そうなこと

ruby {e}^{x} - 1 \mathrm{log}_{e}(1+x) がないので

この計算が必要とする場合、精度で問題になる例がつくれるかもしれない。

C言語Golang \mathrm{log}_{x}(y) がないため

この計算が必要とする場合、精度で問題になる例がつくれるかもしれない。

3.  {x}^{y} {e}^{y \mathrm{log}(x)} の計算結果比較

もともとこの 2つを計算して比べてみようという話でした。

ここでは

 {3}^{559} を調べてみます。

pow(3, 559)  \approx {3}^{559} の計算結果 (倍精度) は

513784961483922578512339180019698197951798048111926076589406942445158467436810000901599831502922238040737504312146764722557085503242923817731646777028029088263023249828259097129213987023170862839861791282352675466344683968308962273231524580225440726788249796697128960

pow(559 * log(3) )  \approx {e}^{559 \mathrm{log} (3)} の計算結果 (倍精度) は

513784961483977278818823467195598637829111367053443355111225788391489700605741415033850055549866587533685229352343834799219693792708265915735692787177526938262762638175334276719142227318286192430271216131815303492760406601731641300588913449010146062071623162601144320

すなわち

誤差  5.47... \times 10^{253} = 54700306484287175900439877313318941517278521818845946331233168931414132250224046944349492947725040197070076662608289465342098004046010149497849999739388347075179589928240295115329590409424849462628026415722633422679027357388868784705335283373365904015360

が生じています。

とんでもない大きさの誤差が生じてしまいました!

考察

 559 \;\mathrm{log}(3) の値が 浮動小数点演算により 近似されたことにより、

次の exp 函数へ誤差が伝播したことが原因と考えられます。

ということは

 \mathrm{log}_{x} (y) = \frac{ \mathrm{log}(y) }{ \mathrm{log}(x) }

という公式がありますが、C言語golang のように標準ライブラリがないものは

このような公式で計算してしまうと

多大な誤差を生む可能性がある ということを意味しているのかもしれません。

上で挙げた expm1 などもそうですが、

やはり標準ライブラリとして一つにまとまっている函数に対しては

丸め誤差への安心感が違います!

そんなお話でした。

素数大富豪だけじゃなく HEX もやろう! ― 紹介編―

小清水 です。 (@curekoshimizu)

twitter 上で 素数大富豪 が大変流行ってるのを目撃しております!

素数大富豪アドベントカレンダー

もつくられており楽しそうです!

HEX の提案!

そんな数学コミュニティに流行って欲しいゲームがあります。

それが HEX です。

この HEX は 数学的背景が大変おもしろいゲーム です。

そしてこれを考えたのは

ゲーム理論で有名な数学者 ジョン・ナッシュ であり、

ビューティフルマインド という映画でも有名な方です。

その ナッシュ学生時代 にこのゲームを考え、

そして 数学的面白さ(証明) を提起しました。

どうですか?ワクワクしてきませんでしょうか。

Abstract

1. HEX とはどんなゲームか?

2. HEX がもつ数学的面白さ

1. HEX とはどんなゲームか?

 n\times n マスで2人でやるゲームになります。

ここでは  5\times 5 マス版HEX で実践例を紹介します。

下のような六角形で敷き詰めた図を考えましょう。

f:id:curekoshimizu:20161222011212p:plain

 5\times 5 マスなので  5\times 5 HEX です。

図が対象的ですので、

先攻・後攻はどちらがどちらでも構いませんが

先攻を紫後攻をピンク

とします。

代わりばんこに塗っていき、

紫は左下から右上につなげると勝ち!

ピンクは右下から左上につなげると勝ち!

というゲームになります。

f:id:curekoshimizu:20161222010516p:plain

つながる?というイメージがわからないかと思うので、

ここで模擬プレイをしてみたいと思います。

HEX 模擬プレイ

先攻真ん中 がおすすめです!

○×ゲームでは真ん中を取ると思いますが、それくらいにオススメです!

奇数  \times 奇数 HEX に真ん中はありませんが、

その場合でも真ん中あたりが有力です。

f:id:curekoshimizu:20161222010541p:plain

この手がすごく強く、後手からすると鬱陶しいです。

続いての後手は適当にうってみましょう。

例えばこうです。

f:id:curekoshimizu:20161222010607p:plain

それを受けて先手はこう打つことにしてみましょう。

f:id:curekoshimizu:20161222010707p:plain

この手は激辛です。

なぜかというと

f:id:curekoshimizu:20161222010727p:plain

「A」・「B」のどちらかに ピンクが打ち込んできても

その反対側を塗ってしまえば

紫は右上までつながることが保証されるからです。

そのため、紫は 真ん中〜右上 までつながったようなものです。

これを実感してもらうと「A」にピンクが打ってきたらその反対側を塗る。

f:id:curekoshimizu:20161222100320p:plain

これを実感してもらうと「B」にピンクが打ってきたらその反対側を塗る。

f:id:curekoshimizu:20161222100327p:plain

真ん中から右上までつながりました!

ピンクの手として次にこの手を考えてみましょう。

f:id:curekoshimizu:20161222100333p:plain

これにはこの手が激痛です。「C」に対して先程と同じことが成立します。

f:id:curekoshimizu:20161222100339p:plain

紫勝利図

紫がつながりました!

f:id:curekoshimizu:20161222100349p:plain

途中経過は飛ばしますが、

例えばこういうつながり方でもOKです。

f:id:curekoshimizu:20161222013746p:plain

2. HEX がもつ数学的面白さ

1. (きちんとやれば)先手必勝ゲームである ことが証明できる

2. 引き分けがなく、盤面を埋めれば絶対に勝敗がきまる ことが証明できる

先手必勝手順の存在性

実は、先手必勝手順の 存在性 を証明できます。

だからといって面白くないと言っているわけではりません。

将棋もチェスも プロの勝率を考えると明らかに 先手優勢です。

具体的な手順が示せるわけではなく、

先手必勝手順があるらしい

という内容になります。

実数係数多項式複素数上で必ず存在する (代数学の基本定理) が、

その厳密解の求め方はわからないと似てますね。

とはいえ、 5\times 5 で見たように

 n\times n HEX で マス が小さいと

実感として 先手有利なことがわかります。

 11 \times 11 HEX までは 先手が 実感としてかなり有利 です。

個人的におすすめの大きさは

 15 \times 15 HEX になります。

f:id:curekoshimizu:20161222100312p:plain

印刷してプレイしてみてください!盛り上がりますよ。

先程のような先手一方戦にはならず、

後手からの反撃が始まります!

引き分けがない

引き分けがないというのがHEXの特徴です!

将棋もチェスも引き分けがあるゲームです。

一方で HEX は 数学的性質から 引き分けがないことが証明できます。

すなわち、

どうあがいても 最後まで塗ると必ずつながっていることになる ということになります。

このことをよくよく考えると、

HEX の逆 ができるということです。

すなわち、

つなげたら負け!というルールでも HEX は成立する

ということになります。

引き分けがないことの証明方法は?

  • Jordanの曲線定理
  • Schauderの不動点定理

を使うと証明ができます。

f:id:curekoshimizu:20161222095832p:plain:w400

準備が大変ですが、いつかこのブログで証明できればいいなと思っております。

今回は 紹介編 ということで証明はでてきません。

ブログ更新を楽しみにしてください!

Summary

数学の証明と紐付いたゲーム (Jordanの曲線定理・Schauderの不動点定理) である、

絶対に引き分けることができないゲーム、

HEX を紹介しました!

Android アプリ でも HEX はできます。

皆さん是非プレイしてみてください!

Newton法でつながるコンピューターと数学の隙間

小清水 (@curekoshimizu) です。

この記事は 日曜数学 Advent Calendar 2016 の 15日目の記事になります。

日曜数学 ということばを各所で聞いて、

僕も数学科時代の わくわく数学 を社会人になっても 感じたい伝えたい!

そう思うようになりました!

f:id:curekoshimizu:20161211183048p:plain:w500

現在はプログラマとして仕事をしており、

「コンピューターの世界と数学の世界」の間を埋めて

面白さを伝えたい! と思い始め、

本当につい最近、

重い腰を上げてこのブログを立ち上げました。

誰かの役に、または、面白く思って頂ければ幸いです。


Abstract

今回の 「コンピューターの世界と数学の世界」 をつなぐお話は Newton法 です。

非常にメジャーなテーマです。

Newton法 を 「除算」計算に使うとどうなるのか?

紹介したい話としては 実数空間だけじゃなく modulo (余り) の世界でも

活きてくる方法なんだ!っていうのが今回のアドベントカレンダーのオチです。

数値計算例もでてくる!

ので数式がわからなくてもフィーリングは伝わると思います!

  1. Newton法とは?
  2. 実数空間の除算とNewton法
  3.  Z/2^{n}Z (  2^{n} で割った余りで考えた空間 ) の除算と Newton法
  4. 今後ブログで伝えていきたいと考えている面白い内容!!!

1. Newton 法

Newton 法は 滑らかな函数  f のある解  \alpha (ただし  f'(\alpha)\neq 0) の近似解を求めるための技法です。

1669 年に Newton がその数年後に Raphson が一般化したため、 Newton-Raphson法 とも呼ばれます。

初期値  x_0 \alpha に近いと

 \displaystyle x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

という漸化式で  \displaystyle \lim_{n\to\infty} x_n = \alpha を得ることができます。

この漸化式は下の図を元に導かれています:

f:id:curekoshimizu:20161211183105p:plain:w300

収束証明

この節に書く証明は、明らかに modulo の世界である  Z/2^{n}Z には適用できない ことにご注意ください。

実数空間の証明は、先の図をみれば感覚的にはOKです!

(読み飛ばしても大丈夫です!)

きちんとした証明も

 f \alpha の近傍で  C^{2} 級であれば

 \alpha の閉近傍で  f'\neq 0 となるようとり、

 x_0 がこの近傍にいるようにすれば、

次式より帰納的に  x_n はこの閉近傍に収まっており、さらにこの不等式が証明できます。

 \displaystyle | x_{n+1} - \alpha |

 \displaystyle = \left| x_n - \alpha - \frac{f(x_n)}{f'(x_n)} \right|

Taylor の定理より

 \displaystyle = \left|x_n - \alpha -  \frac{ f(\alpha)+f'(x_n)(x_n - \alpha)+f''(\theta) \left( x_n - \alpha \right)^{2} }{f'(x_n)}   \right|

 \displaystyle = \left| \frac{ f''(\theta) \left( x_n - \alpha \right)^{2} }{f'(x_n)}   \right|

区間内で、連続函数の最大値が存在することから

 \displaystyle \leq \left\lVert \frac{ f'' }{f'}  \right\rVert_{\infty} \left( x_n - \alpha \right)^{2}

これは

精度が2乗ずつよくなっていく ことを示しています!

実際それを「除算」の例でみてみましょう。

2. 実数の除算 と Newton法

 \displaystyle f(x) = \frac{1}{x} - a に Newton法の公式を当てると

 \displaystyle x_{n+1} = x_{n} (2-ax_{n})

となります! 減算と乗算だけで除算の近似値が得られることを意味します!

生活の中で考えても、例えば、 小学生時代を思い出しましょう!

除算ってすごく計算が大変だったと思います。

とくに に何が立つのかという計算が実に厄介です。

次の計算は8がたつかな?と 予測 し、あー間違ってた…。

9じゃないかーという風に

筆算で除算を計算していると思います。

コンピューターにとっても除算は激しく大変です。

コンピューターでも上のような 予測 を入れながら計算する方法があるのですが、

Newton 法による計算は 予測なしに愚直に計算するだけ で得られます。

ためしに  a = 1.5 として 逆数の  2/3 の近似値を得てみましょう。

n  x_n 正しく得られた桁数
0 0.5 0桁
1 0.625 1桁
2 0.6640625 2桁
3 0.666656494140625 4桁
4 0.6666666665114462375640869140625 9桁
5 0.6666666666666666666 3052659425048318553308490 … 19桁
6 0.66666666666666666666666666666666666666470750 … 38桁

と先の証明の 精度が2乗ずつよくなっていく 様子がわかります。

続いて  Z/2^{n}Z へ全く同じ方法で計算するか計算してみましょう。

3.  Z/2^{n}Z の除算と Newton法

もちろん注意すべきは  Z/2^{n}Z ではないので一般に除算はできません。

しかしながら

 ax = 1 \quad\mathrm{mod}\; 2^{n}

という等式を  a が奇数に限り考えることができます。

偶数は明らかに解なしです。(左辺が偶数、右辺が奇数になってしまうためです。)

そのため、お話としては 逆元を Newton 法で求めてみようという話になります。

皆さん大好きな  a = 691

 Z/2^{32}Z の世界で Newton 法 にかけてみますと

次のようになります。

ここでは  x_0 = 1 として

 \displaystyle x_{n+1} = x_{n} (2-ax_{n})

を計算してきます。

n  x_n 2進法表現 正しく得られた桁数
0 1 00000000000000000000000000000001 1桁
1 4294966607 11111111111111111111110101001111 2桁
2 3966933707 11101100011100101001101011001011 4桁
3 3802448251 11100010101001001100000101111011 8桁
4 2476047483 10010011100101010111110001111011 16桁
5 2660269179 10011110100100000111110001111011 32桁
6 2660269179 10011110100100000111110001111011 32桁

最終的に収束して 不動点 が得られました。

精度が2乗ずつよくなっていく 様子も同じです。

実際

 691\times 2660269179 = 1 \quad\mathrm{mod}\; 2^{32}

となりますので 逆元が得られました。

たまたまなのでしょうか。

ちなみにですが、691 を持ち出した理由を知りたい方は

tsujimotterさんのブログ を御覧ください!

収束証明

 x_0 = 1 であり a は奇数であるから、

 ax_0 - 1 = 0 \quad\mathrm{mod}\; 2^{2^{0}}

ある n で

 ax_n - 1 = 0 \quad\mathrm{mod}\; 2^{2^{n}}

が満たされていたと仮定すると  a x_{n} - 1 = 2^{2^{n}} N となる 整数  N がとれる。

 a x_{n+1} - 1

 =ax_{n} (2 - ax_{n}) -1

 = (2^{2^{n}}N + 1) (1 - 2^{2^{n}}N) - 1

 = 2^{2^{n}} 2^{2^{n}} N^{2}
 = 2^{2^{n+1}} N^{2}

ゆえに

 ax_{n+1} - 1 = 0 \quad\mathrm{mod}\; 2^{2^{n+1}}

が満たされ、帰納法により示された。

なぜ  Z/2^{n}Z 系を扱った?

コンピューターでは 整数 を表現する際に

 Z/2^{n}Z として 扱われることが多いです。

特に  n = 8, 16, 32, 64 が頻出で、

プログラミング言語 を勉強しても登場します。

上の証明をよくよく見てみると、

これを発展させると p進数 にも応用できること簡単にわかりますよね!

Summary

Newton 法除算 につかうと

実は 「減算・乗算」 だけで計算できる式を導ける!

そんな、実数の 接線 と関係の深い Newton法

接線とは関係なさそうな

整数の世界 でも関係がある!

おまけにコンピュータとも関係あるなんて!

4. 今後ブログで伝えていきたいと考えている面白い内容!!!

このように 数学の知識とコンピューターの知識 の両方が登場する内容を書いて、

橋渡し的なブログになればいいなと思っております。

最近興味があるのは 浮動小数 という

コンピューターで実数を扱う際に

よく使われる計算体系を 数学的に見つめ直してみたい と思っており

その内容が多いかもしれません。

浮動小数点がどれくらい厄介?

浮動小数点体系では、例えば基本的にこうです…。

  •  (a + b) + c \neq a + (b + c)
  •  (a \times b) \times c \neq a \times (b \times c)

群・環・体を知っている人は

性質の悪さに 泣きそうになりますよね...。

非可換体なんて目じゃない難しさです。

この 難しい演算体系Excel などでは標準的に使われています。

下の内容は 三角函数の計算が定義通り計算できない!? 悲しい話になります:

math-koshimizu.hatenablog.jp

あなたも コンピュータと数学 に興味持ってみませんか?

以上、コンピューターと数学の間の話でした。


明日の 日曜数学 Advent Calendar 2016 の記事は

Toshiki Takahashi さんによる「リープグラフと複素確率」です。

確率って聞くだけでわくわくしちゃいます!

Excelでおかしな計算結果になった問題の正解値を求める

小清水 (@curekoshimizu) です。

以前紹介した記事に

math-koshimizu.hatenablog.jp

がありました。この内容はざっくりいいますと

 a=77617, b=33096 の場合の

 \displaystyle 333.75 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) +5.5 b^8 + \frac{a}{2b}

の値を Excel で計算すると 1.17260394...

真の解は -0.827396...

であると紹介しました。

正しい結果の計算方法を載せなかったために

正しい結果 を信じてくれない人がいましたので、

その計算方法を載せます。

どうやって計算する?

浮動小数点数で計算したのが間違いです。

この演算の登場したのはすべて 有理数のみ であることから

今回の場合、本当に正確に計算結果を求めたいのであれば

有理数ですべて計算すればいい のです。

当然といえば当然です。

計算してみよう

python有理数ライブラリを使って計算するのが簡単です。

#!/usr/bin/env python

from fractions import Fraction

a = Fraction(77617)
b = Fraction(33096)

ret = Fraction(333.75)* b*b*b*b*b*b + a*a*(Fraction(11) * a*a * b*b -b*b*b*b*b*b -Fraction(121)*b*b*b*b-Fraction(2)) + Fraction(5.5)*b*b*b*b*b*b*b*b + a/(2*b)

print ret
print float(ret)

この計算結果は

-54767/66192
-0.827396059947

なので有理数

 \displaystyle -\frac{54767}{66192} = -0.82739605994\cdots

が正解値なのがわかります。

明らかに 負の数 なのに Excel は正の数 を返してきたわけですが、

これで Excel の計算結果に不安 を感じてもらえたでしょうか?

短い記事でしたが以上です。

コンピュータでおかしなことになる計算例 (3) ―たった1度の型変換―

小清水 です。(@curekoshimizu)

今回は C言語のsin関数 で遊んでみましょう。

最後のオチまで気づかない人がいるかも!?

というわけで、

あーなるほどねと思っても読み進めて欲しいです!

今回紹介する例

  • sinf は 単精度用 (float) の sin関数
  • sin は倍精度用 (double) の sin関数
#include <stdio.h>
#include <math.h>

int main(int argc, char** argv)
{
    double x = 1.0e30;

    printf("%.10lf\n", sinf(x));
    printf("%.10lf\n", sin(x));

    return 0;
}

出力結果

-0.7911634445
0.0093314689

そう。全然結果が違う のです!

これは sin と sinf の違いによるものなのでしょうか?

すなわち、 倍精度のsin関数単精度のsin関数

実装の問題 なのでしょうか。

それは違います。

sinf は 単精度(float) への入力のみ受け付けることから、

暗黙の型変換 による 丸め誤差が1度 発生します。

検証プログラム

倍精度用 (double) 側にも 単精度(float) への型変換を加えて確認してみましょう。

#include <stdio.h>
#include <math.h>

int main(int argc, char** argv)
{
    double x = 1.0e30;

    printf("%.10lf\n", sinf(x));
    printf("%.10lf\n", sin((float)x)); // <--- sin(x) から sin((float)x)) にした

    return 0;
}

出力結果

-0.7911634445
-0.7911634385

ほぼほぼ近いです。先程出力された 0.0093314689 とは全然違います。

このことから sin と sinf に大きな違いがあったわけではありません。

違いがあったのは、

倍精度から単精度へのたった一度の丸め誤差が入ったかどうか です。

    printf("%.10lf\n", sinf(x)); // <--- x の float への暗黙の型変換

よくよく考えると当たり前だが少し考えさせられる

実はこのプログラム

 1.0 \times 10^{30} とういとてつもなく大きな値からスタートしています。

最終的に単精度(2進24桁)へ丸めが入りますので、

math-koshimizu.hatenablog.jp

ここで説明したルールの RN(x) が入ります。

これにより丸められた値は

 1.0000000150474662\times 10^{30} です。

この値で近似されるので

その誤差なんと

 1.50474662\times 10^{22} です。

そのため、ここまで大きな誤差が生じる例となっていた ということになります。

なので 当たり前といえば当たり前 の例かもしれません。

しかしながら、この短いプログラムでそんなことに気がつきましたでしょうか。

桁数が大きい場合の型変換は思っているより近くない

ため、プログラムによっては大ダメージになる ということに注意が必要です。

プログラミングをしているときに大きな桁になっているかも…ということに

気づくことはそうそうできないため、

少し考えさせられる例かもしれません。

最後のオチ

 1.0\times 10^{30} の sin計算を 単精度で計算したのが間違い!

 1.0\times 10^{30} の sin計算 は倍精度側で出力された値の

0.0093314689...

と捉えた人は間違いですよ。

 1.0\times 10^{30} の sin計算 は倍精度側で出力された値の

0.0093314689... ではありません!

倍精度大好きっ子はよくこの結論に達しがちなのですが、

    double x = 1.0e30; // <--- 実数から浮動小数点数への丸め

ここも暗黙の型変換があります!

さっきと同じ話で こんな大きな数字の倍精度誤差も、単精度同様 尋常 じゃありません。

 \mathrm{RN}_{倍精度}(1.0\times 10^{30}) = 1000000000000000019884624838656

の sin計算 が 0.0093314689... というのはあってます!

これがタイトルの「たった一度の型変換」の意味になります。

ちなみに、 1.0\times 10^{30} の sin計算結果は 負の数になりますので、全然違いますね。

勘違いしてしまった人は

math-koshimizu.hatenablog.jp

ここを読んで是非丸め誤差の勉強を一緒にしましょう!


そんな、当たり前のようで少し勘違いしそうな例を思いついた、 昼のひとときでした。

コンピュータでおかしなことになる計算例 (2) ― 三角函数を定義通りに ―

小清水 (@curekoshimizu) です。

以前の記事に Excelで計算が破綻する例 を紹介しました。

math-koshimizu.hatenablog.jp

紹介した数式は次でしたが、

 a=77617, b=33096
 \displaystyle 333.75 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) +5.5 b^8 + \frac{a}{2b}

実に 人工的 な感じがします。

Abstract

今回紹介する話は 数学上自然な定義 なのに

それが うまくいかない という話になります。

お題は 皆さん馴染みのある 三角函数 (sin) です。


三角函数の定義 といえばこの問題ですよね!

f:id:curekoshimizu:20161206083612p:plain:w450

‘99年東京大学入試問題。意表を突かれる入試問題です。

この問題への回答ではありませんが、

sin 函数の定義 ってなんでしょう。

大学以降、sin 函数は次で定義されることが多いです。

sin 函数の定義

 \displaystyle \mathrm{sin}(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^n \frac{x^{2n+1}}{(2n+1)!} + \cdots \quad (x \in \mathrm{R})

この式は 0周りの Talylor 展開ですが、

x : 実数全体 で (広義一様) 収束 する式になります。 (実は 複素数全体まで拡張できる)

非常に非常に自然な定義 です!

sin 函数を定義通りに コンピュータで計算してみよう

この式を使って計算をしてみると、

例えば次のようなプログラムになるかと思います。

(剰余項に関する評価は次の節で)

#include <stdio.h>
#include <stdint.h>
#include <math.h>

double sin_taylor(double x)
{
    const double minus_xx = -x*x;

    const double term_tho = 1.0e-16;

    /* Taylor展開第1項までは直接計算 (x) */
    double term = x;
    double ret = x;

    /* Taylor展開第3項以降をループで計算 */
    for (uint32_t k = 1; ; k++) {

        uint32_t n = 2*k + 1;

        /* Taylor 展開第x^n=x^{2k+1}項の計算式: (-1)^k x^n/n! */
        term *= minus_xx / ( (n - 1.0) *n );

        if (fabs(term) <= term_tho) {
            break; /* 充分収束したことを示す条件 (次の節を参照) */
        }

        ret += term;
    } 

    return ret;
}

int main(int argc, char** argv)
{

    for (uint32_t ix = 0; ix <= 50; ix += 5) {

        double ret = sin_taylor((double)ix);

        printf("%d, %lf\n", ix, ret);
    }

    return 0;
}

結果をまとめると

x sin 計算結果
0 0.000000
5 -0.958924
10 -0.544021
15 0.650288
20 -0.988004
25 -0.132351
30 -0.988004
35 -0.428035
40 1.069746
45 111.576759
50 -25547.371714

0 から離れると 非常におかしい結果に変わっていきます。

実数上では  | \mathrm{sin}(x) | \leq 1 満たさないといけません。

しかし  \mathrm{sin}(40.0), \mathrm{sin}(45.0), \mathrm{sin}(50.0) は明らかにおかしいです。

評価を間違えたのでしょうか?

評価方法を数学的に考える

 x > 0 : fixed としたとき、

Taylor の定理より

 \displaystyle \mathrm{sin}(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^{n-1} \frac{x^{2n-1}}{(2n-1)!} + (-1)^{n} \frac{\mathrm{sin}^{(2n+1)}(\theta)}{(2n+1)!} x^{2n+1}

となる  \theta \in (0, x) が存在することになります。

ここでは記法として  \mathrm{sin}^{(2n+1)} \mathrm{sin} (2n+1) 階導函数としました。

 \displaystyle \left| \mathrm{sin}(x) - \left( x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^{n-1} \right) \right| = \left| (-1)^{n} \frac{\mathrm{sin}^{(2n+1)}(\theta)}{(2n+1)!} x^{2n+1} \right| \leq \frac{x^{2n+1}}{(2n+1)!}

という評価を得ますので、

 \frac{x^{2n+1}}{(2n+1)!} \leq 10^{-16} となる  n まで計算すれば

 \displaystyle x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^{n-1} \frac{x^{2n-1}}{(2n-1)!}

は sin と少なくとも絶対誤差  10^{-16} となる解を得るはずです。

しかし計算結果は そうなっていません

Summary

何が言いたかったかというと

実数上で証明された式

コンピュータの中の 浮動小数点環境 において

意味のある数式になるとは限らない

ということです。

「0 から遠いのにそんな式を使うから悪い」

という人はいるかと思います。

確かに 0 まわりで近似された Taylor 展開の式なので、

その主張もわかります。

しかし今回の数式については、

どんな実数値 であっても 実数空間において 収束 が約束されていたはずです。

結局のところ、

どんな数式だったら収束するのか? と心配する必要がありそうです。

次の主張として

 \mathrm{sin}(x+2\pi) = \mathrm{sin}(x)

という数式を使って 0周りに近づけて計算すべきだとも思いますが、

この  2\pi がそもそも 浮動小数点数 では exact に表現できませんが大丈夫でしょうか。

この関係を使って0に近づけても やはりおかしなことは起こりうる ということを

を紹介する予定です。


この記事を読んで

浮動小数点・丸め が気になった方は是非こちらをどうぞ!

math-koshimizu.hatenablog.jp

(追記)

コンピュータでおかしなことになる計算例シリーズ第三弾できました!

math-koshimizu.hatenablog.jp

浮動小数点の丸めの方向と性質 (1)

小清水です。

前回の記事では 浮動小数点数による計算のため に、

計算が はちゃめちゃ・わやくちゃ になる例を紹介しました。

math-koshimizu.hatenablog.jp


Abstract

今回は そんな浮動小数点数

もう少し踏み込んだ話をするための準備記事になります。

実数を浮動小数点に近似する 丸め について詳しく踏み込んでみます。

4つの代表的な丸めを定義して考察を加えていきます。

四捨五入のような考え方がイケてない理由なども登場します。

浮動小数点数のざっくりした理解

2進p桁の浮動小数点数とは 次のように表現できる数:

 (-1)^s \times \left( 1 + \frac{ m_1 }{ 2^1} + \frac{ m_2 }{ 2^2} \dots + \frac{ m_{p-1}}{2^{p-1}} \right) \times 2^e \\
= (-1)^s \times (1.m_1 m_2 \dots m_{p-1})_2 \times 2^e \\
= (符号部) \times (仮数部) \times 2^{(指数部)}

になります。ただし、

 s, m_i \in \{0, 1\}, e \in \mathrm{Z}

例えば、2進3桁の浮動小数点の分布は次のようになります。

f:id:curekoshimizu:20161204201749p:plain

注意.

ここでざっくりと書いたのは、

ここでは e の区間を整数全体としているものの、

本来は特定の数の間で限定されていたり、

非正規化数の話が抜けているためです。

いつか正確な内容をブログにアップする予定です。

実数は基本的に浮動小数点数で表せない

上図で表した通り  \pi浮動小数点数で表せません。

そのため、浮動小数点数へと近似してあげる必要があります。

この、実数から浮動小数点への近似函数丸め といいます。

代表的な丸め函数

代表的な丸め方法としてこの4つがあります:

  • RN(x) : x の 最近接偶数方向丸め (round to nearest even)
  • RU(x) : x の 正の無限大方向丸め (round up)
  • RD(x) : x の 負の無限大方向丸め (round down)
  • RZ(x) : x の 0方向丸め (round toword zero)

f:id:curekoshimizu:20161204200848p:plain

先程の例の場合、

 \pi は RN・RD・RZ の場合に

 (1.00)_2 \times 2^1 = 3

へと丸められます。これは、

RN は最も近い浮動小数点数に丸められる から、

RD はその数以下の浮動小数点のうちで最も近い浮動小数点数に丸められる から、

RZ はその数が正であれば RD と同じ丸め方向に・負であれば RU と同じ方向に丸められる からです。

また、RU の場合には

 (1.01)_2 \times 2^1 = 3.5

へと丸められます。これは RU がその数以上の浮動小数点のうちで最も近い浮動小数点に丸められる からです。

RN の 最近接偶数方向丸めとは?

上の説明では RN の 最近接 たる所以、

つまり、最も近い浮動小数点に丸められるというところを説明しました。

問題は 偶数方向 の意味になります。

こちらは次のような 浮動小数点数のちょうど真ん中の数 に対する規則になります。

f:id:curekoshimizu:20161204213832p:plain

この図の説明通り、

最近接の2数の浮動小数点のうち、仮数部の最下位が 偶数(0) の方を選択する

という意味になります。

そのため、先の 3.25 だけでなく 2.75 についても 3 側に丸められるということになります。

f:id:curekoshimizu:20161204213843p:plain

偶数側ばかり選んで不平等に感じる?それは次を御覧ください。

RN は なぜ四捨五入のような規則ではない?

我々は 四捨五入 という規則を知っており、

それに近い規則を考えれば 零捨一入 という考えが普通に感じます。

このようにすれば 偶数値ばかり選択されず不平等さは生じません。

しかしながら、

そもそも 四捨五入 はある困った性質があり、

浮動小数点の丸めでは 四捨五入・零捨一入 というのは採用されませんでした。

最後に 2進法の場合 を示すとして、

ここでは我々にとって 身近な四捨五入直感に反してしまう例 を紹介して、

日常生活でも使えるウンチクにしたいと思います。

四捨五入を採用すると直感に反してしまう例

実数空間では

 a \quad + b - b \quad + b - b \quad + b - b\quad = a

というのは常識だと思います。

今回の例は、この等式が怪しくなるという話です。

10進3桁で考えることとして、

 a = (1.00)_{10} \times 10^0 = 1, b = (5.55)_{10} \times 10^{-1} = 0.555

a に b を足して引くということを繰り返すとどんどん大きくなるという例を示します。

四捨五入の丸め方法を R(x) と書くことにしますと

 
\mathrm{R}( \mathrm{R}(1.00 + 5.55\times 10^{-1}) - 5.55\times 10^{-1}) \\
= \mathrm{R}( \mathrm{R}(1.555) - 5.55\times 10^{-1}) \\
= \mathrm{R}( 1.56 \times 10^0 - 5.55\times 10^{-1}) \\
= \mathrm{R}( 1.005 \times 10^0 ) \\
= \mathrm{R}( 1.01 \times 10^0 )

a+b-b をやって 1.00 が 1.01 に大きくなりました。

さらに繰り返しましょう。

 
\mathrm{R}( \mathrm{R}(1.01 + 5.55\times 10^{-1}) - 5.55\times 10^{-1}) \\
= \mathrm{R}( \mathrm{R}(1.565) - 5.55\times 10^{-1}) \\
= \mathrm{R}( 1.57 \times 10^0 - 5.55\times 10^{-1}) \\
= \mathrm{R}( 1.015 \times 10^0 ) \\
= \mathrm{R}( 1.02 \times 10^0 )

先の結果に +b-b をすると 1.00 が 1.02 まで大きくなることが分かりました。

さらに繰り返しましょう。

 
\mathrm{R}( \mathrm{R}(1.02 + 5.55\times 10^{-1}) - 5.55\times 10^{-1}) \\
= \mathrm{R}( \mathrm{R}(1.575) - 5.55\times 10^{-1}) \\
= \mathrm{R}( 1.58 \times 10^0 - 5.55\times 10^{-1}) \\
= \mathrm{R}( 1.025 \times 10^0 ) \\
= \mathrm{R}( 1.03 \times 10^0 )

先の結果に +b-b をすると1.03 になりました。

実はまだまだずーっと続くのですが、このあたりで止めましょう。

これを drifting といいます。

注意. 2進法3桁で同じような例をつくるには

 a = (1.11)_2 \times 2^0 , b = (1.11)_2 \times 2^{-1}

で同じことをすればいいです。

Theorem 1. (RN は 連続した加減算で drifting が発生しない)

RN : 最近接偶数方向丸めに対して

 \mathrm{RN}(\mathrm{RN}(\mathrm{RN}(\mathrm{RN}(a+b) -b)+b)-b) = \mathrm{RN}(\mathrm{RN}(a+b)-b)

ということが知られています。

四捨五入・零捨一入という規則ではなく、

偶数方向丸めというものが採用されているのはこの性質のためです。

小清水は現時点でこの証明方法を知りません。

証明を知っている方は是非教えてください。

丸め誤差と誤差の関係をまとめる

再び浮動小数点の分布図を見ながら考えてみると、

f:id:curekoshimizu:20161204221456p:plain

RU・RD・RZ による丸め誤差隣合う浮動小数点の間隔 未満となります。

 \mathrm{RU}(\pi) はほぼ対岸側の浮動小数点に丸められていることを考えれば、

限りなくもう一方の浮動小数点数に近い実数を考えればこの事実がわかるかと思います。

一方、RN による丸め誤差は 隣り合う浮動小数点のうち 近い方 であるから、(隣り合う浮動小数点数の間隔)/2 が最大となります。

これらをまとめると

Proposition 1 (丸めと絶対誤差の関係)

2進p桁の浮動小数点への実数の丸めによる絶対誤差は

 {2}^{n} \leq |x| <  {2}^{n+1} \; (n \in \mathrm{Z}) ならば

  •  \displaystyle \left| x - \mathrm{RN}(x) \right| \leq {2}^{n-p} \left(= \frac{1}{{2}^{p-1}}\cdot {2}^{n} \cdot \frac{1}{2} \right)

  •  \displaystyle \left| x - \mathrm{RU}(x) \right| <  \displaystyle {2}^{n-p+1} \left(= \frac{1}{{2}^{p-1}}\cdot {2}^{n}  \right)

  •  \displaystyle \left| x - \mathrm{RD}(x) \right| <  \displaystyle {2}^{n-p+1} \left(= \frac{1}{{2}^{p-1}}\cdot {2}^{n}  \right)

  •  \displaystyle \left| x - \mathrm{RZ}(x) \right| <  \displaystyle {2}^{n-p+1} \left(= \frac{1}{{2}^{p-1}}\cdot {2}^n  \right)

という性質を満たす。

これがわかると簡単に次の内容もわかります!

Proposition 2 (丸めと相対誤差の関係)

2進p桁の浮動小数点への実数の丸めによる相対誤差は

  •  \displaystyle \left | \frac{ x - \mathrm{RN}(x) }{ x } \right| \leq {2}^{-p}

  •  \displaystyle \left | \frac{ x - \mathrm{RU}(x) }{ x } \right| <  {2}^{-p+1}

  •  \displaystyle \left | \frac{ x - \mathrm{RD}(x) }{ x } \right| <  {2}^{-p+1}

  •  \displaystyle \left | \frac{ x - \mathrm{RZ}(x) }{ x } \right| <  {2}^{-p+1}

という性質を満たす。

(証明)

 {2}^{n} \leq |x| <  {2}^{n+1} \; (n \in \mathrm{Z}) ならば  \displaystyle \left| x - \mathrm{RN}(x) \right| \leq {2}^{n-p} であるから

 \displaystyle \frac{ \left| x - \mathrm{RN}(x) \right| }{ |x| } \leq \frac{ {2}^{n-p} }{ {2}^{n} } = {2}^{-p}

となるため。他も同様。

最後に

ここまで代表的な丸めについて説明し

その誤差について説明しました。

ここまで説明できましたので

  • 誤差 についてさらに詳しく
  • RN・RU・RD・RZ の使い道

などを記載予定だったのですが、

長くなりすぎたので

次回のブログ記事にまわしたいと思います。

コンピュータでおかしなことになる計算例 (1)

小清水です。初ブログになります。

カタストロフィックに 計算結果がおかしくなる 例を紹介したいと思います。

Excel を使用した例と C言語の例です。


お膳立て:(飛ばしてOK)

コンピュータの中で 実数 は一般にどのように表されているのでしょうか?

浮動小数点数 と呼ばれる数で表されることが多いです。

例えば、 \pi はどうでしょうか?

この円周率は コンピュータ内部において、

 
\begin{align}
+(1.10010010000111111011011)_2 \times 2^1 
\end{align}

という 浮動小数点数 で近似されます。(単精度を使用した例)

2進法ではわかりづらいため、この数を 10進法 でかいてみると

3.141592 7410125732 という数になります。

本来の値は 3.141592 65358979 … なので 誤差 が生じたことになります。

この誤差を 丸め誤差 というわけです。

丸め誤差って何なのか、詳しく知りたい方はこちらをどうぞ!

math-koshimizu.hatenablog.jp


本題!おかしな例

今回の話は、 コンピュータで近似計算をしているため、 おかしな結果がでるよという話です。

とはいえ、色んな所で使われている Excel でもそんなことできるんでしょうか?

Excel による例

Excel 2011 (Mac) を使って次の式を計算したいと思います:

 a=77617, b=33096
 \displaystyle 333.75 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) +5.5 b^8 + \frac{a}{2b}

f:id:curekoshimizu:20161127225544p:plain

図の通り、 1.17260394 という計算ができました!

しかしながら実際の -0.827396... です。

もはや近いというレベルではありません。符号も違います。

(追記)

この -0.827396... の導出方法を解説しました。

math-koshimizu.hatenablog.jp

C言語による例

プログラマー向けに C言語によるプログラムも。よく使われる double で計算してみます。

#include <stdio.h>

double rump1(double a, double b)
{
    return 333.75 * b*b*b*b*b*b 
        + a*a * (11.0 * a* a* b* b - b*b*b*b*b*b - 121.0 * b*b*b*b - 2.0) 
        + 5.5 * b*b*b*b*b*b*b*b + a /(2.0*b);
}

double rump2(double a, double b)
{
    double b2 = b*b;
    double b4 = b2*b2;
    double b6 = b4*b2;

    return 333.75 * b*b*b*b*b*b
        + a*a * (11.0 * a*a* b*b - b6 - 121.0 * b*b*b*b - 2.0)
        + 5.5  *b*b*b*b*b*b*b*b + a/(2.0*b);
}

double rump3(double a, double b)
{
    double b2 = b*b;
    double b4 = b2*b2;
    double b6 = b4*b2;
    double b8 = b4*b4;
    double a2 = a*a;

    return 333.75 * b6
        + a2 * (11.0 * a2* b2 - b6 - 121.0 * b4 - 2.0)
        + 5.5  *b8 + a/(2.0*b);
}

int main(int argc, char** argv)
{
    const double a = 77617.0;
    const double b = 33096.0;

    printf("%f\n", rump1(a, b)); /* ---->  1.172604 */
    printf("%f\n", rump2(a, b)); /* ----> -2361183241434822606848.000000 */
    printf("%f\n", rump3(a, b)); /* ----> -1180591620717411303424.000000 */

    return 0;
}

と、計算の仕方によってはさっきの Excel のようになったり、 激しく大きな負数になったりします。

困りましたね。

端折りますが、単精度・倍精度・4倍精度でも計算は合いません。


まさにカタストロフィック!

この元ネタは 「Algorithms for verified inclusion」という1988年の論文です。 興味がある方は是非御覧ください。

この記事の魅力としては

Excel を使ってもC言語を使っても 1.17260394 という計算結果を得られたけど、本当の答えは-0.827396...という圧倒的に変な結果が得られた というわけで

世の中の人はExcelに対してとても不安になったに違いない!

浮動小数点数って難しい…って思ってもらうための話でした。

第二弾もあるよ。いつか書きます!

第二弾も書いたよ!

math-koshimizu.hatenablog.jp