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

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

素数大富豪だけじゃなく 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