文字列で 10 進任意精度整数の四則演算 その 4 ベルヌーイ数

ネタとはいえ、巨大整数を計算できるようになったら、ベルヌーイ数を計算してみたくなるものです。

プログラミング言語年表 - Wikipedia

チャールズ・バベッジが提案した最新鋭の機械である解析機関についての記事をイタリアの数学者ルイジ・メナブレアが執筆し、1842年から1843年の9ヶ月間にエイダ・ラブレスがそれを翻訳した。この記事の中で彼女はこの機械でベルヌーイ数を計算する完全なプログラムを掲載した。これは世界初のコンピュータプログラムであると言われている。

ベルヌーイ数は浮動小数点演算で求めるには向いておらず、都合の良いことに有理数なので、分数計算によって正確な値が求まります。分数をその都度、約分しながら計算すると遅いので、ここでは、分母・分子をそれぞれ求めていって、最後に約分することにします。

sub fraction_reduction {
    my($xnum, $xden) = split m{/}, $_[0], 2;
    my $snum = $xnum =~ s/\A-//msx ? '-' : '';
    my $sden = $xden =~ s/\A-//msx ? '-' : '';
    my $sign = $snum ^ $sden;
    my $d = _gcd($xnum, $xden);
    return join '/', $sign . _div($xnum, $d), _div($xden, $d);    
}

sub _gcd {
    my($m, $n) = @_;
    $m =~ s/\A-//msx; $n =~ s/\A-//msx; 
    return $n if $m eq '0';
    return $m if $n eq '0';
    if (_cmp($m, $n) < 0) {
        my $t = $m; $m = $n; $n = $t;
    }
    while ($n ne '0') {
        my($q, $r) = _divmod($m, $n);
        $m = $n;
        $n = $r;
    }
    return $m;
}

ベルヌーイ数の計算方法は様々なやりかたがありますが、ここでは書きやすいけれど最も遅い式で求めてみます。

sub bernoulli_number {
    my($n) = @_;
    my($bn, $bd) = ('0', '1'); # 求めたいベルヌーイ数の分子と分母
    my $f = '1';
    for my $j (1 .. $n) {
        my($sn, $sd) = ('0', '1'); # 入れ子の級数を求める
        for my $m (1 .. $n - $j + 1) {
            my $r = '1';
            for my $i ($m .. $m + $j - 1) {
                $r = _mul($r, $i);
            }
            $sn = _add(_mul($sn, $m + $j), _mul($sd, $r)); $sd = _mul($sd, $m + $j);
        }
        my $p = '1';        # j**n
        for my $i (1 .. $n) {
            $p = _mul($p, $j);
        }
        $f = _mul($f, -$j); # (-1)**j * j!
        $sn = _mul($p, $sn); $sd = _mul($f, $sd);
        $bn = _add(_mul($bn, $sd), _mul($bd, $sn)); $bd = _mul($bd, $sd);
    }
    return fraction_reduction(join '/', $bn, $bd);
}

say bernoulli_number(22); # B[22]=>854513/138

22 番目を求めるのに数秒かかるかと予想していたら、一瞬間があった程度で、意外とすぐに計算結果が出てきました。