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

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

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 などもそうですが、

やはり標準ライブラリとして一つにまとまっている函数に対しては

丸め誤差への安心感が違います!

そんなお話でした。