読者です 読者をやめる 読者になる 読者になる

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

コンピューター・数学 に関することを書きます

最近接偶数方向丸めの統計的性質を調べてみた

小清水 (@curekoshimizu) です。

以前のこちらの記事では、

math-koshimizu.hatenablog.jp

浮動小数点の丸め について説明しました!

今回は 「絶対誤差」・「相対誤差」の性質 について紹介したいと思います。

簡単に復習

この記事で登場した

RN(x) : 最近接偶数方向丸め函数

次のような 最も近い浮動小数点数で近似する函数 でした。

例えば、次のような 2進数3桁精度 のとき、

 \mathrm{RN}(\pi) は最も近い浮動小数点数である、

 (1.10)_{2} \times {2}^{1} = 3 を返します。

f:id:curekoshimizu:20170318224031p:plain

ちょうど浮動小数点数の 2数の間を取った場合の RN の定義や、

その他の丸め誤差函数については上のブログ記事を御覧ください。

今回のお話は、この函数  \mathrm{RN} の詳細な挙動についてです!

RN と 誤差

 \mathrm{RN}(x) は 真値  x とどの程度誤差があるのでしょうか?

浮動小数点の誤差は

  • 絶対誤差 (Absolute error)
  • 相対誤差 (Relative error)

の 2つで語られることが多いです。

真値 ( {x}_{Exact} ) と 近似値 ( x_{Approx}) があった際に

絶対誤差は

\big| x_{Exact} - x_{Approx} \big|

絶対誤差は

 \displaystyle\left| \frac{ x_{Exact} - x_{Approx} }{ x_{Exact}} \right|

で定義されます。

RN と 絶対誤差

最近接遇数方向丸め函数 RN の 絶対誤差 は次のような挙動となります (この例は2進数3桁精度の例):

f:id:curekoshimizu:20170318224511p:plain

指数部の大きさによってその誤差が異なります。

 \mathrm{AbsErr}_{e,p}(x) = \mathrm{T}_D (x; {\beta}^{e}) + \mathrm{T}_D (x; {\beta}^{e}+D) + \mathrm{T}_D (x; {\beta}^{e}+2D) + \dotsb +  \mathrm{T}_D (x; {\beta}^{e+1}-D)

参考:用いた計算結果

 \displaystyle \int_{-\infty}^{\infty} T_D (x; S)\; dx = \frac{{D}^{2}}{4}

 \displaystyle \int_{-\infty}^{\infty} x T_D (x; S)\; dx = \frac{{D}^{2}}{8} (D+2S)

 \displaystyle \int_{-\infty}^{\infty} {x}^{2} T_D (x; S)\; dx = \frac{{D}^{2}}{96} (7 {D}^{2}+24DS+24{S}^{2})

 \displaystyle \int_{-\infty}^{\infty} \frac{1}{x} T_D (x; S)\; dx = -D \log{ \left(\frac{D}{2}+S \right)} + D\log{(D+S)}  +S\log{S} - 2S\log{ \left(\frac{D}{2}+S \right)} + S\log{(D+S)}

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

小清水 です。 (@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 の計算結果に不安 を感じてもらえたでしょうか?

短い記事でしたが以上です。