awk で多倍長演算

404 Blog Not Found:javascript - Math.BigInt で多倍長整数演算 にインスパイヤされて、xgawk の MPFR 機能を用いて計算を行ってみます。 MPFR は整数演算ではなく任意精度浮動小数点演算 (arbitrary-precision binary floating-point computation) なので、意味合いは元と異なりますが、awk でこうした演算ができるのは嬉しい限りです。

元ネタのものと全て揃えるのは疲れるので、一部省略していますが、500 までのフィボナッチ数も求めています。

#! /usr/local/bin/xgawk -f
# bigint.awk
# 多倍長演算
# usage: xgawk -f bigint.awk

@load mpfr

BEGIN {
    MPFR_BASE = 10;
    MPFR_PRECISION = 400;

    # そもそも awk で数値として扱えないため、ダブルクォートで括って
    # 文字列にしておく必要がある。
    x = "1234567890123456789012345678901234567890";
    y = "12345678901234567890";

    print mpfr_add(x, y);
    print mpfr_sub(x, y);
    print mpfr_mul(x, y);
    print mpfr_div(x, y);
    print mpfr_cmp(x, y);
    print mpfr_neg(x);
    print mpfr_div(fact(101), fact(100));

    for (i = 1; i <= 500; i++) {
        print fibonacci_memo(i);
    }
}

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

# fibonacci_mem()
#   in:    number
#   out:   fibonacci number
function fibonacci_memo(n) {

    if (n <= 2) {

        return 1;

    } else if (memo[n] != 0) {

        return memo[n];

    } else {

        memo[n] = mpfr_add(fibonacci_memo(n - 2), fibonacci_memo(n - 1));

        return memo[n];
    }
}

ここでフィボナッチ数は AWK Users JP :: フィボナッチ数列 で使ったメモ化を使っていますので、非常に高速に処理することができます。

$ time xgawk -f bigint.awk
1.2345678901234567890246913578024691357800000000000000000000000000000000000000000000000000000000000000000000000000000000000E39
1.2345678901234567890000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E39
1.5241578753238836750342935775034293577501905199875019052100000000000000000000000000000000000000000000000000000000000000000E58
1.0000000000000000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E20
1
-1.2345678901234567890123456789012345678900000000000000000000000000000000000000000000000000000000000000000000000000000000000E39
1.0100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2
1
1
2.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
3.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
5.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
8.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
1.3000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1
2.1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1
3.4000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1
5.5000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1
8.9000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E1
1.4400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2
2.3300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2
<snip>
8.6168291600238450732788312165664788095941068326060883324529903470149056115823592713458328176574447204501000000000000000000E103
1.3942322456169788013972438287040728395007025658769730726410896294832557162286329069155765887622252129412500000000000000000E104
xgawk -f bigint.awk  0.04s user 0.01s system 32% cpu 0.146 total

MPFR_PRECISION を上げているため、ちょっと見づらいのが難点です。

ちなみに gawk のロードマップが示されていますが、この中に MPFR も予定されているようです。

tag_xgawk.png