標準添付の Math:: モジュールのご紹介

基礎寄りの話になりますが、Perl5 でもコアモジュール（最初から入っているモジュール）で大きな数の計算などを簡単に行うことができます。

Perlより後発の言語がたくさん出てきて、Perlを主に書いていると「数値計算なら○○です」といったことを言われることが多いです。でも、ちょっとしたことであればPerlでも書けますし、CPANモジュールで結構まかなえたりします。

Perlの数学や数値計算のお話は、意外とメジャーじゃないのかなと思っ て、今回このテーマで書くことにしました。

より高度なモジュールについては割愛して、Perl を入れてすぐ使えるコアモジュールに焦点をあててみることにしました。

Math:: モジュール

Math::なんとか といった名前空間のモジュールをここでは「Math:: モジュール」と書いています。

perldoc perlmodlib を見ると、数学っぽいコアモジュールをいくつも見つけることが出来ます。しかもだいたいは10年くらい前の Perl5.8 から使うことが出来るものです。

• Math::BigFloat
• Math::BigInt
• Math::BigInt::Calc
• Math::BigInt::CalcEmu
• Math::BigInt::FastCalc
• Math::BigRat
• Math::Complex
• Math::Trig

Math::Big で始まるものは、主に精度を上げてくれるものです。Int は整数、Float は浮動小数点、Rat は有理数です。

これには対応するプラグマ

• bigint
• bigfloat
• bigrat

があって、対応するものの精度を透過的に広げます。

簡単に試してみるのであれば

$ perl -E 'say 2**100'
1.26765060022823e+30
$ perl -Mbigint -E 'say 2**100'
1267650600228229401496703205376

といった感じです。だいたい有効数字14桁だったのが、巨大な数でも（メモリなどの資源が潤沢であれば）どんとこいとなります。

Math::BigRat は、小数点ではなく分数を保持するモジュールです。巨大な有理数の計算だけでなく、ちょっとした分数を計算したいというときにも使えます。

$ perl -MMath::BigRat -E 'my ($x, $y) = map { Math::BigRat->new($_) } qw(1/3 1/2); say $x+$y;'
5/6

Math::Complex は複素数を扱います。虚数単位として定数 i がインポートされます。

$ perl -MMath::Complex -E 'my $z = i; say $z**2;'
-1+1.22464679914735e-16i

本当は虚数部分が 0 なんですが、ちょっと惜しいですね。

$ perl -MMath::Complex -E 'my $z = 2 + 2*i; say $z**2;'
4.89858719658941e-16+8i

こちらも実数部分が 0 になってほしいのに惜しい。

新しめの Ruby だと、2i といった表現で純虚数を書いたりできるようで、Perlにも取り入れられないかな（他の書き方とバッティングしないし）と思います。

Math::Trig は三角関数ライブラリです。標準でも sin cos atan2 が用意されていますが、その他の三角関数を提供してくれるものです。

最初に言ったことと相反するようですが、コアモジュールだと必要最小限の数学ができるといった感じです。

この後は各々が求める数学的題材は結構違ってくると思うので、Perlがコアで何を導入するかといったことはなかなか難しいのかなと感じます。ある人は行列の固有値を計算したいのかもしれないし、また別の人は微分方程式を解きたいのかもしれない。何でも取り入れてもらえると嬉しい場面もあることは確かですが、数学を一切使わない多くの人にとって perl のコアが無意味に重たくなることを考えたら、次は CPAN を頼ってねというのは Perl5 の良い戦略なのかもしれないですね。

素数判定プログラムを書いてみる

ありきたりな題材ですが、以前必要があって素数判定プログラムを書いてみました。素数判定といっても偉そうなものでなく、いわゆる エラトステネスのふるい と呼ばれる、地道に割っていくというものです。

#!/usr/bin/perl

use strict;
use warnings;

# 超巨大な整数を扱う
use bigint;
use Math::BigInt;

# memoize すれば CPUコストは減らせるけれど、強烈な大きさの数の検査をするとメモリを食い潰しそうなのでなんとも
use Memoize;
memoize('is_prime');

# エラトステネスの篩
sub is_prime {
    my $n = shift;
    my $cb = shift; # 状況確認フック
    my $max_try = int(sqrt($n));
    # $max_try があまりに大きすぎると範囲演算子で膨大なリストを作れないので
    for ( my $try_divider = 2;  $try_divider <= $max_try; $try_divider++ ) {
        next if !is_prime($try_divider); # 合成数で試す必要はない
        if ( $cb && ref $cb eq 'CODE' ) {
            $cb->($try_divider);
        }
        if ( $n % $try_divider == 0 ) {
            return; # 合成数だ
        }
    }
    return 1; # 素数だ
}

巨大な数の検査にも使えるよう、数値リテラル（ソースコードに直接書く数値）なども巨大な数に耐えうるよう、bigint プラグマと Math::BigInt モジュールを使っています。is_prime に与えられる $n も Math::BigInt オブジェクトであることを想定しています。

割り算がいつまでたっても終わらないことを観察するために、状況確認フックを渡せるようにもなっています。

これを書いたのは、2 から 89 までの素数で $p^p$ という形を作った和が素数であるという話題を計算しようとしている方がいたからです。

詳細な形で書くと、$n$ 番目の素数を $p_n$ と書くと、$2 = p_1$ から $89 = p_{24}$ までについての以下の和

\sum_{n=1}^{24} p_n^{p_n}

が素数であるかということです。

私も素数についての知識が少なかったのですが、実際に bigint プラグまで精度気にせず割り算したらどうなるのか書いてみたときに書いたのが上記のサブルーチンでした。

その時書いたもの全部は以下のソースコードになります。

#!/usr/bin/perl

use strict;
use warnings;

# 超巨大な整数を扱う
use bigint;
use Math::BigInt;

# memoize すれば CPUコストは減らせるけれど、強烈な大きさの数の検査をするとメモリを食い潰しそうなのでなんとも
use Memoize;
memoize('is_prime');

# エラトステネスの篩
sub is_prime {
    my $n = shift;
    my $cb = shift; # 状況確認フック
    my $max_try = int(sqrt($n));
    # $max_try があまりに大きすぎると範囲演算子で膨大なリストを作れないので
    for ( my $try_divider = 2;  $try_divider <= $max_try; $try_divider++ ) {
        next if !is_prime($try_divider); # 合成数で試す必要はない
        if ( $cb && ref $cb eq 'CODE' ) {
            $cb->($try_divider);
        }
        if ( $n % $try_divider == 0 ) {
            return; # 合成数だ
        }
    }
    return 1; # 素数だ
}

sub keta {
    my $x = shift;
    #return int(log($x)/log(10)) + 1;
    return length($x);
}

my $current_check_keta = 1;
sub progress {
    my $x = shift;
    my $keta = keta($x);
    if ( $current_check_keta < $keta ) {
        print "$keta 桁に突入しました\n";
        $current_check_keta = $keta;
    }
}

#my $sum = 0;
my $sum = Math::BigInt->new(0);
for my $p ( grep { is_prime($_) } (2..89) ) {
    my $pp = Math::BigInt->new($p);
    # $p_pow_p = $p ** $p;
    my $p_pow_p = $pp->bpow($pp);
    #print "$p^$p = $p_pow_p\n";
    #$sum += $p_pow_p;
    $sum->badd($p_pow_p);
}

print "$sum は素数らしい\n";

my $keta = keta($sum);
print "$keta 桁ある数、検査してみます。\n";

$SIG{ALRM} = sub {
    print "終わらないと思うし、中止しますか？ [Y/n]\n";
    chomp( my $input = <STDIN> );
    exit if length $input == 0 || $input =~ /^y(?:es)?/i;
};
alarm 3600;

if ( is_prime($sum, \&progress) ) {
    print "うんやっぱり素数だ\n";
} else {
    print "違うんじゃない？\n";
}

この和、200桁近くあることは対数を使って分かるので、こりゃ手元の MacBook Air で素因数分解するレベルじゃないなぁと思っていたのですが、当然ながら全く終わる気配がありませんでした（状態確認フックを挟むとよくわかります）。

学生時代は数学をやっていたのですが、最近までコンピュータを使った数値計算というものに関心を持ってこなかったので、このあたりの知識はほとんど素人です。

最近知って感心した素数判定アルゴリズムは、ピエール・デザルトによる1999年の成果です。

p_n > n(\log n + \log\log n - 1), \, n>6

ある場面においてこの成果を効果的に使うことにより、かなり素数判定が速くなるそうです。

その他の Perl の数値計算モジュール

CPAN を探せば Math:: から始まるモジュールがいっぱいあって、簡単なことであればだいたい間に合うでしょう。また、GMP や PDL、また Pari/GP などの高度な数値計算ライブラリの Perl バインディングなども用意されています。

ちょっとした巨大な数の計算であれば、bigint などのコアモジュールを使って計算できることを知っておくと、古い Perl くらいしか入っていないサーバでも巨大な数の計算を手軽にさせることができて便利なこともあると思います。

よく勉強会の懇親会で「Perlで数値計算をやりたいんですが、○○といったモジュールが用意されていなくて…」という人の話を聞いたりします。でも、そういうことを考えている人が何人もいることが分かれば、誰か手の空いているPerlプログラマーが書いてくれるというのは夢物語でもない気がするので、とりあえず数値計算を大好きなPerlでやりたいけれど…という方は、Twitterなどで想いを言葉にしてみるといいかもしれません。