exp関数の (実装の) 限界を体感

c 数値計算

@shuyo さんがブログを書いた数日後に奇しくもロジスティック回帰を実装してexpのオーバーフローではまっていたので，expの限界について，自分自身で体感してみることにした．

@shuyo さん曰く，だいたい710くらいでinfになりまっせ，とのこと．ふむふむこれは勉強になる．

"exp オーバーフロー" とか検索すると，MSDNでは709.ほげほげでオーバーフローしまっせとかいう記述を見つけることができる．


こんなコードを書いて限界を調べてみた．

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int
main (void)
{

  double x;
  for (x = 0.0; x < 1000.0; x += 1.0) {
    printf("x=%f: exp(x)=%e\n", x, exp( x ));
  }

  return 0;
}

実行

% gcc -lm exp.c && ./a.out

...
x=706.000000: exp(x)=4.091704e+306
x=707.000000: exp(x)=1.112241e+307
x=708.000000: exp(x)=3.023383e+307
x=709.000000: exp(x)=8.218407e+307
x=710.000000: exp(x)=inf
x=711.000000: exp(x)=inf
x=712.000000: exp(x)=inf
x=713.000000: exp(x)=inf
...

おぉ，やはり709と710の間でオーバーフローしてしまった!
オーバーフローするとinfになるらしい．


今度は，アンダーフローのチェック

  double x2;
  for (x2 = 0.0; x2 > - 1000.0; x2 -= 1.0) {
    printf("x=%f: exp(x)=%e\n", x2, exp( x2 ));
  }


実行結果

...
x=-744.000000: exp(x)=9.881313e-324
x=-745.000000: exp(x)=4.940656e-324
x=-746.000000: exp(x)=0.000000e+00
x=-747.000000: exp(x)=0.000000e+00
...

今度は，745と746の間でアンダーフローを起こして0になっている．仕様によって違いはあるかもしれないが，0を返すので一見問題がなさそう．デバッグしづらそう．．．



値が小さいときは，どこらへんでサチるかチェック

  double y;
  for (y = 1.0; y > 0; y /= 10.0) {
    printf("y=%e: exp(y)=%e\n", y, exp( y ));
  }

結果

y=1.000000e+00: exp(y)=2.718282e+00
y=1.000000e-01: exp(y)=1.105171e+00
y=1.000000e-02: exp(y)=1.010050e+00
y=1.000000e-03: exp(y)=1.001001e+00
y=1.000000e-04: exp(y)=1.000100e+00
y=1.000000e-05: exp(y)=1.000010e+00
y=1.000000e-06: exp(y)=1.000001e+00
y=1.000000e-07: exp(y)=1.000000e+00
y=1.000000e-08: exp(y)=1.000000e+00
y=1.000000e-09: exp(y)=1.000000e+00
...

大体，10^(-6)と10^(-7)で1に落ち着いていることがわかる．


数式では簡単に肩の上に載せたり降ろしたりということをしてしまうけれど，計算機はこんなに大変な思いをしているんだよ，ということを少しだけ体感することができた．

さて，expがオーバーフローするよウワァァンというときにはlogsumexpの出番なのだけれど，これについては既に丁寧な解説があるので割愛．というか logsumexp も言葉は知っているけれど一度も使ったことがないので解説できないのが本音．


さて，実際にexpをどう計算するかということまで踏み込んだ良記事．奥が深い．．．