2009年02月06日 11:30 [Edit]

perl - Math - Mersenne Twister を Pure Perlで

camel

Refactor してみた。

PurePerlでメルセンヌ・ツイスタな話。Ver.0.1 - 永字八法
これに対する解決として、use bigintプラグマの導入と言う激烈馬鹿な手段を選択してしまい、MTの利点である「高速性」を大きく損なうことに成功。ダメじゃん!

まずは結果を。

[Run via CodePad]
#!/usr/local/bin/perl;
package MTpp;
use strict;
use warnings;
use Math::BigInt;

use constant N          => 624;
use constant M          => 397;
use constant MATRIX_A   => 0x9908b0df;
use constant UPPER_MASK => 0x80000000;
use constant LOWER_MASK => 0x7fffffff;

my @mt  = ();
my $mti = N + 1;

sub mult32 {
    use Math::BigInt;
    my $bi = Math::BigInt->new( $_[0] );
    $bi->bmul( $_[1] );
    $bi->band(0xffffffff);
    return $bi->numify();
}

sub init_genrand {
    my $s = shift;
    $mt[0] = $s & 0xffffffff;
    foreach my $mti ( 1 .. N - 1 ) {
        my $n = $mt[ $mti - 1 ];
        $n ^= ( $mt[ $mti - 1 ] >> 30 );
        $mt[$mti] = mult32( 1812433253, $n ) + $mti;
    }
}

sub genrand_int32 {
    my $y = 0;
    my @mag01 = ( 0, MATRIX_A );
    if ( $mti >= N ) {
        init_genrand(5489) if ( $mti == N + 1 );
        foreach my $kk ( 0 .. ( N - M - 1 ) ) {
            $y = ( $mt[$kk] & UPPER_MASK ) | ( $mt[ $kk + 1 ] & LOWER_MASK );
            $mt[$kk] = $mt[ $kk + M ] ^ ( $y >> 1 ) ^ $mag01[ $y & 0x1 ];
        }
        foreach my $kk ( ( N - M ) .. ( N - 1 - 1 ) ) {
            $y = ( $mt[$kk] & UPPER_MASK ) | ( $mt[ $kk + 1 ] & LOWER_MASK );
            $mt[$kk] =
              $mt[ $kk + ( M - N ) ] ^ ( $y >> 1 ) ^ $mag01[ $y & 0x1 ];
        }
        $y = ( $mt[ N - 1 ] & UPPER_MASK ) | ( $mt[0] & LOWER_MASK );
        $mt[ N - 1 ] = $mt[ M - 1 ] ^ ( $y >> 1 ) ^ $mag01[ $y & 0x1 ];
        $mti = 0;
    }
    $y = $mt[ $mti++ ];
    $y ^= ( $y >> 11 );
    $y ^= ( $y << 7 ) & 0x9d2c5680;
    $y ^= ( $y << 15 ) & 0xefc60000;
    $y ^= ( $y >> 18 );
    return $y;
}

if ( $0 eq __FILE__ ) {
    print genrand_int32(), "\n" for ( 1 .. $ARGV[0] || 1000 );
}
1;

計算でbigintが必要だったのは、結局かけ算だけだったので、そこのみMath::BigIntを明示的に使って計算した。一端@mtが埋まってしまえば、あとはbigintは不要なのでまっとうな速度になる。

あと、genrand_int32()にオリジナルのmt19937ar.cと不一致な点がいくつかあって、結果が一致しなかったのでそこも書き直し。結果が一致することも確認済み。

Enjoy!

Dan the Mersenne Twisted

追記:

OS X v10.5.6 / Perl 5.10.0 で Math::Random::MT のコンパイルがこけるなあ....

追^2記:

報告済みの模様


この記事へのトラックバックURL