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

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

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

コンピュータでおかしなことになる計算例 (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