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

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

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