数学ライブラリのあの関数，どんな計算してるの？

現代の大抵のプロセッサというのは四則演算を行う命令しか用意されていない．あってもそれらを組み合わせた積和演算ぐらい．下手すると乗算命令と除算命令が省かれていたりもする．だけど，$\sin$ とか $\log$ とか $\exp$ とか計算できてしまう．なぜか？

四則演算だけで計算する

とりあえず今は自然指数関数を作りたい．指数は実数としよう．
値を求めるための手法はいくつかあると思うが，ここではマクローリン展開を使うこととしよう．大体，こういうのはWikipediaを見れば分かるようになってる．

テイラー展開#マクローリン級数の一覧より

e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} \qquad \text{for all } x

無限級数の形になってはいるが，$n$ を大きくしていけば， $\frac{x^n}{n!}$ は多分 $0$ に近づいていくので，適当なところで計算をやめて大丈夫そう．$n$は自然数なので，$x^n$も$n!$も普通に乗算だけで計算できる．あとは除算して，加算するだけなので，四則演算だけで計算できている．これで完璧だね！

・・・・・・
・・・・
嘘つけ．

これ，$x=0$ だと $n=0$ で $0^0$ が出てきて計算できないじゃん．そこだけ場合分けが必要．

e^x = \begin{cases}
1  & (x = 0)\\
\sum_{n=0}^{\infty} \frac{x^n}{n!} & (\text{otherwise})
\end{cases}

四則演算だけで自然指数関数の計算ができそうなことが分かった．しかも，意外と簡単そうな計算だ．これならすぐ作れそう！

昔見た何かのライブラリも，アセンブリを追っていくと，やはりマクローリン展開を使って計算しているような雰囲気だった．

Wikipediaを見れば分かるように，$\sin$ も $\ln$ も，マクローリン展開が既に知られているので，同様に実装することができそうだ（今回はやらないけど）．

pythonで実装

今回はpythonで自作の自然指数関数my_expを実装してみる．もちろん，ライブラリは使わない．ついでなので，最近覚えたunittestも使ってみる．

まずは作ってみる

とりあえず作ってみた．dがMY_PREC以下になったら，計算結果が十分収束したと判断して，ループを終了する．importは省略してるわけではなく，importする必要が無かったので本当に何も書いていない．

mymath.py

MY_PREC = 0.00000000000001

def my_exp(x):
    if x == 0:
        return 1.0

    sum = 1.0
    d = 1.0
    n = 1
    while d > MY_PREC:
        d = d * x / n
        sum += d
        n += 1
    return sum

ね，簡単でしょ？

ユニットテストでmy_expの出来具合を確認する．比較する値は，mathライブラリを使って求める．求まった2つの値の差の絶対値が十分小さければ，同じ値とみなして，my_expは正しく計算できているものとする．さすがに，まったく同じ計算精度で計算ができるとは期待していないため，このようにした．x=0の場合は厳密に1であるべきなので，差の絶対値ではなく，assertEqualを使って，1.0と比較した．

test_mymath.py

import unittest
from mymath import my_exp
import math

TEST_PREC = .0000000000001

class TestMyExp(unittest.TestCase):
    def test_my_exp(self):
        self.assertEqual(1.0, my_exp(0))
        self.assertLess(self.__cmp(1/2), TEST_PREC)
        self.assertLess(self.__cmp(1), TEST_PREC)
        self.assertLess(self.__cmp(-1), TEST_PREC)

    def __cmp(self, x):
        t = math.exp(x)
        myexp = my_exp(x)
        
        d = t - myexp
        return d if d > 0 else -d

if __name__ == "__main__":
    unittest.main()

このユニットテストを走らせると，次のような結果を得る．

> python .\test_mymath.py
F
======================================================================
FAIL: test_my_exp (__main__.TestMyExp)
----------------------------------------------------------------------
Traceback (most recent call last):
  File ".\test_mymath.py", line 11, in test_my_exp
    self.assertLess(self.__cmp(-1), TEST_PREC)
AssertionError: 0.36787944117144233 not less than 1e-13

----------------------------------------------------------------------
Ran 1 test in 0.001s

FAILED (failures=1)

x=0のときに問題ないのはソースを見れば明らかなので，わざわざテストしなくて良いような気がしてしまうが，そうではない．バグというのはテストしなかった部分に存在するということを忘れてはならない．全ての値でテストするのは現実的でない以上，何をテストデータとして選択するかは，よく考えて選ばなければならない．¹

x=1はネイピア数そのものが求まるわけだが，ちゃんと計算できているようだ．

しかし，x=-1で値が大きくずれているようだ．

xが負の場合

よく考えると，xが負の場合，dも負の値を取り得る．しかしながら，whileの条件はdが負の場合を想定していない．n=1でいきなりdは負の値になるので，ループがすぐに終了してしまうことが分かる．

$e^x=\frac{1}{e^{-x}}$だから，$x$が負の場合は$e^{-x}$を計算して，後で逆数にしてやることにする．

xが負の場合に対応したmy_expは次のようになる．

mymath.py

def my_exp(x):
    if x == 0:
        return 1.0

    sig = 0
    if x < 0:
        sig = 1
        x = -x

    d = 1.0
    sum = d
    n = 1
    while d > MY_PREC:
        d = d * x / n
        sum += d
        n += 1

    return sum if sig == 0 else 1 / sum

ついでなので，ユニットテストのデータも少し追加してみる．

test_mymath.py

        self.assertLess(self.__cmp(2), TEST_PREC)
        self.assertLess(self.__cmp(-2), TEST_PREC)
        self.assertLess(self.__cmp(100), TEST_PREC)

実行結果は次のとおり．

F
======================================================================
FAIL: test_my_exp (__main__.TestMyExp)
----------------------------------------------------------------------
Traceback (most recent call last):
  File ".\test_mymath.py", line 14, in test_my_exp
    self.assertLess(self.__cmp(100), TEST_PREC)
AssertionError: 1.9807040628566084e+28 not less than 1e-13

----------------------------------------------------------------------
Ran 1 test in 0.003s

FAILED (failures=1)

x=100で非常に大きく値が違っていることが分かる．指数関数の計算を間違えたのだろうか？

ユニットテストの改良

計算結果を見てみると，次のようになっていた．

>>> mymath.my_exp(100)
2.6881171418161336e+43
>>> math.exp(100)
2.6881171418161356e+43

似たような結果が得られている．しかし，指数部が43なので，単純に差の絶対値を計算すれば，それはもう非常に大きな値になってしまう．

ここではユニットテストの比較のやり方を改良する．mathライブラリで得られた結果に対して，差がどれくらいの比になるかを計算する．この比が十分小さければ，my_exp関数は正しく実装できているものとする．これであれば，計算結果の大小にかかわらず，定量的な評価ができそうだ．

test_mymath.py

    def __cmp(self, x):
        t = math.exp(x)
        myexp = my_exp(x)
        d = t - myexp
        d = d if d > 0 else -d
        res = d / t
        return res

実行してみる．

.
----------------------------------------------------------------------
Ran 1 test in 0.000s

OK

これで完成だ！
　
　
　

・・・・・・本当に？

もっとユニットテスト

ユニットテストにx=708の場合を追加してみる．

test_mymath.py

        self.assertLess(self.__cmp(708), TEST_PREC)

実行すると，次の通り・・・

やばい，結果が返ってこないだと！？pythonの対話モードで確認してみる．

>>> math.exp(708)
3.023383144276055e+307
>>> mymath.my_exp(708)

mathライブラリはちゃんと計算できているが，my_expは結果が返ってこない．これの原因は，dの値が大きくなりすぎるから．$n$が小さい間は$n!$に比べて，$x^n$の増加が非常に大きく，dが収束する前に無限大になってしまう．そうなると，dは収束しなくなるのでループを抜けられなくなってしまう．実際に確認すると，n=681でdは無限大になっていた．

ここでは改良のために，$e^{a+b}={e^a}{e^b}$という性質を使う．指数を整数+小数に分けて，整数部は普通に乗算のループで，小数部をマクローリン級数で求めて，最後にこの2つの値を乗算する．あと，ネイピア数は数値計算だけで求められるといっても，さすがにこれは定数として定義しておくことにする．

mymath.py

MY_PREC = 0.00000000000001
my_e = 2.718281828459045

def my_exp(x):
    if x == 0:
        return 1.0

    sig = 0
    if x < 0:
        sig = 1
        x = -x

    x1 = int(x)
    x2 = x - x1

    res = 1.0
    for i in range(x1):
        res *= my_e

    if x2 > 0:
        d = 1.0
        sum = d
        n = 1
        while d > MY_PREC:
            d = d * x2 / n
            sum += d
            n += 1
        res *= sum
    
    return res if sig == 0 else 1 / res

この場合，xを整数にすると，マクローリン級数の計算が全くなされないことになるので，ユニットテストでは必ず小数点以下が0でない値で確認しなければならない．テストをさらに追加する．

test_mymath.py

        self.assertLess(self.__cmp(708.1), TEST_PREC)
        self.assertLess(self.__cmp(708.5), TEST_PREC)
        self.assertLess(self.__cmp(708.7), TEST_PREC)

このユニットテストは成功した．
今度こそ完成だね！

さらにユニットテスト

次のテストを追加してみる．

test_mymath.py

        self.assertLess(self.__cmp(710), TEST_PREC)

実行結果は次の通り．

E
======================================================================
ERROR: test_my_exp (__main__.TestMyExp)
----------------------------------------------------------------------
Traceback (most recent call last):
  File ".\test_mymath.py", line 19, in test_my_exp
    self.assertLess(self.__cmp(710), TEST_PREC)
  File ".\test_mymath.py", line 22, in __cmp
    t = math.exp(x)
OverflowError: math range error

----------------------------------------------------------------------
Ran 1 test in 0.001s

FAILED (errors=1)

今度はmathライブラリの方でオーバフロー例外が発生している．ということは，おそらく，my_expの方もオーバフローが発生する．my_expにも例外出力が必要だし，ユニットテストもオーバフローの例外処理が必要そうだ．

my_exp関数は，整数部の計算の後と，最終的な結果を求めた後の2カ所で，無限大の判定をし，そうであるなら例外が発生するようにした．まあ，最後の1カ所だけでも良いのだが，ムダにループ回す必要も無いかと思って2カ所にした．

mymath.py

def my_exp(x):
    if x == 0:
        return 1.0

    sig = 0
    if x < 0:
        sig = 1
        x = -x

    x1 = int(x)
    x2 = x - x1

    res = 1.0
    for i in range(x1):
        res *= my_e

    if res == float('inf'):
        raise OverflowError

    if x2 > 0:
        d = 1.0
        sum = d
        n = 1
        while d > MY_PREC:
            d = d * x2 / n
            sum += d
            n += 1

        res *= sum

    if res == float('inf'):
        raise OverflowError
    
    return res if sig == 0 else 1 / res

ユニットテストは例外処理を追加した．mathライブラリもmy_expもオーバフローが発生したなら，値は一致したものとする．どちらか一方がオーバフローならば，一応，不一致としておく．両方ともオーバフローが発生しなければ，今まで通り．$e^x$は全域で正の値を取ることが分かっているので，オーバフロー時は結果を-1とすることで，容易に判別できる．

test_mymath.py

    def __cmp(self, x):
        try:
            t = math.exp(x)
        except OverflowError:
            t = -1
        
        try:
            myexp = my_exp(x)
        except OverflowError:
            myexp = -1
        
        if t < 0:
            if myexp < 0:
                res = 0
                print('inf x={0}:math.exp my_exp'.format(x))
            else:
                res = 1
                print('inf x={0}:math.exp'.format(x))
        else:
            if myexp < 0:
                res = 1
                print('inf x={0}:my_exp'.format(x))
            else:
                d = t - myexp
                d = d if d > 0 else -d
                res = d / t
        return res

最終的に，次のようなテストを行った．

test_mymath.py

    def test_my_exp(self):
        self.assertEqual(1.0, my_exp(0))
        self.assertLess(self.__cmp(1/16), TEST_PREC)
        self.assertLess(self.__cmp(1/8), TEST_PREC)
        self.assertLess(self.__cmp(1/4), TEST_PREC)
        self.assertLess(self.__cmp(1/2), TEST_PREC)
        self.assertLess(self.__cmp(1), TEST_PREC)
        self.assertLess(self.__cmp(-1), TEST_PREC)
        self.assertLess(self.__cmp(2), TEST_PREC)
        self.assertLess(self.__cmp(-2), TEST_PREC)
        self.assertLess(self.__cmp(3.1415), TEST_PREC)
        self.assertLess(self.__cmp(10), TEST_PREC)
        self.assertLess(self.__cmp(100), TEST_PREC)
        self.assertLess(self.__cmp(707), TEST_PREC)
        self.assertLess(self.__cmp(708), TEST_PREC)
        self.assertLess(self.__cmp(709), TEST_PREC)
        self.assertLess(self.__cmp(-709), TEST_PREC)
        self.assertLess(self.__cmp(709.5), TEST_PREC)
        self.assertLess(self.__cmp(709.8), TEST_PREC)
        self.assertLess(self.__cmp(709.9), TEST_PREC)
        self.assertLess(self.__cmp(710), TEST_PREC)

結果は次の通り．

inf x=709.8:math.exp my_exp
inf x=709.9:math.exp my_exp
inf x=710:math.exp my_exp
.
----------------------------------------------------------------------
Ran 1 test in 0.001s

OK

x=709.8, 709.9, 710で両方ともオーバフローが発生していることが分かる．片方だけオーバフローというのは，少なくともこのテストデータの中には無かったようだ．本当にそういう場合があるのかどうかまでは検証しない．単にめんどすぎるので．

まとめ

Pythonを使って，mathライブラリを使わず，四則演算だけで自然指数関数my_expを実装した．計算結果はmathライブラリのexpと比較して，計算精度の細かな違いはあるものの，ほぼ同じものができたと考えられる．

ついでなので，最近覚えた，unittestを使って検証を行った．ユニットテストのおかげで，かなり楽に検証が行えるようになったし，デバッグにも十分役に立ってくれた．元の数式は単純に見えても，それを実際にプログラムを組んでテストを重ねていくと，元の数式だけでは見えてこなかった問題点をいくつも発見することができた．

my_exp関数の完成し具合としてはまあこんなもんだろう．細かい改良点はあるだろうが，別にこれで何をするってわけでもなし，これで実装完了としたい．

このプログラムの場合，最初のif文を削除するとバグに・・・と思ったが，実は正しい結果が得られる．$x\ne0$ なら，$n=0$ のとき，必ず $\frac{x^n}{n!}=1$ になるため，最初からsum=1.0(=d)となるようにハードコーディングしてある．したがって，この後のループでは $n\geq1$ のときのみを計算していることになる．ところが，$x=0$ の場合，ハードコーティング部は特に条件が無いので　$x\ne0$ のときとまったく同じ処理が進み，ループへ到達する．このループでdはすぐに0となり，sum(=1.0)に加算しても結果は変わらず，明らかにd < MY_PRECなので，すぐにループを抜けてしまう．当然，$n\geq2$ で $\frac{x^n}{n!}=0$ なのは明らかなので，計算しなくても特に支障は無い．したがって，sum=1.0のまま結果を返すことになるが，これはx=0のときに欲しい値と厳密に一致することになる． ↩