マクローリン展開で sin(x) を求める

sinｘの数値計算 -OKWave にインスパイヤされてマクローリン展開で sin(x) を求めてみます。 テイラー展開 - Wikipedia によるとマクローリン展開とはテイラー展開の一部であり、sin(x) のマクローリン展開も同ページに掲載されています。

#! /usr/local/bin/nawk -f
# maclaurin_series.awk
# sin(x) のマクローリン展開
# usage: nawk -f maclaurin_series.awk

BEGIN {
    x = ARGV[1];

    print "         n:     Maclaurin:       Builtin";

    # 無間ループ
    for (n = 0;; n++) {
        sin_x = sin_x + ((-1) ^ n) / fact(2 * n + 1) * (x ^ (2 * n + 1));

        printf("%10d: %13.10f: %13.10f\n", n, sin_x, sin(x));

        if (abs(sin_x - tmp_sin_x) < 1e-8) {
            exit;
        }
        tmp_sin_x = sin_x;
    }
}

# fact():   階乗を計算する
#   in:     数値 (n)
#   out:    入力された数値の階乗を返す
function fact(n) {
    if (n >= 2) {
        return n * fact(n - 1);
    } else {
        return 1;
    }
}

# abs():    絶対値を返す
#   in:     数値 (n)
#   out:    数値 (n) の絶対値
function abs(n) {
    return (n >= 0) ? n : -n;
}

特に難しいところはないと思いますが、ひとつ前の計算結果を何らかの形で保持しなければならないくらいでしょうか。

$ nawk -f maclaurin_series.awk 0
         n:     Maclaurin:       Builtin
         0:  0.0000000000:  0.0000000000

$ nawk -f maclaurin_series.awk 1
         n:     Maclaurin:       Builtin
         0:  1.0000000000:  0.8414709848
         1:  0.8333333333:  0.8414709848
         2:  0.8416666667:  0.8414709848
         3:  0.8414682540:  0.8414709848
         4:  0.8414710097:  0.8414709848
         5:  0.8414709846:  0.8414709848
         6:  0.8414709848:  0.8414709848

$ nawk -f maclaurin_series.awk 3.14159265
         n:     Maclaurin:       Builtin
         0:  3.1415926500:  0.0000000036
         1: -2.0261201123:  0.0000000036
         2:  0.5240439130:  0.0000000036
         3: -0.0752206116:  0.0000000036
         4:  0.0069252742:  0.0000000036
         5: -0.0004451566:  0.0000000036
         6:  0.0000211462:  0.0000000036
         7: -0.0000007692:  0.0000000036
         8:  0.0000000260:  0.0000000036
         9:  0.0000000031:  0.0000000036
        10:  0.0000000036:  0.0000000036

今回は隣に組込関数の sin() を併記してありますので、ターゲットとなる値もわかりやすいとでしょう。