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

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

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

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

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

そんなお話でした。