TSperl 発言

極悪です。

HIDAKA Takahiro さんの [TSperl:60] Re: Let's野ぷ野ぷ（その３） から

>　いや、それだけじゃなくて、例えば RAND_MAX が 2^15 の
>環境で、rand() / (RAND_MAX + 1.0) で 0 〜 1 の乱数を求めているの
>だとすると、どうがんばっても rand() として 0 〜 32767 しか
>出てこなくて、結果として 32768 通りの数値しか取得できない場合が
>ありますよね。
>
>　このときでも、周期性としては 2^100 の乱数だって作れるはず
>ですが、モンテカルロ法としては、 32768 回以上の繰り返しは
>無意味になるはずじゃないかと思ったのです。

なるほど、すばらしいです。実は乱数ってスカスカだったんですね。

昔書いたスクリプトが @nifty の会議室に残ってました。その後ハ
ードディスクがとんで いろんなものがなくなっちゃったんですが、
こういうときはネットがちょっとしたバックアップになって便利だ
ぁ。



--^ PerlPi.pl
#!/usr/local/bin/perl
use strict;
use Math::BigInt;
my $VERBOSE = 0;



#
# 展開公式の係数と、atan の中の値の分母とを順に指定する。
# π/4 =   8 atan 1/10 -   atan 1/239 -  4 atan 1/515 より
#   π =  32 atan 1/10 - 4 atan 1/239 - 16 atan 1/515
#        -A-       -B-  -C-       -D-  -E-        -F-
#
my @param   = (
     32,    # A
     10,    # B
     -4,    # C
    239,    # D
    -16,    # E
    515,    # F
);;;



#
# y = a * atan(1/x) を計算する。
#   = a[(1/x^1)/1 -(1/x^3)/3 +(1/x^5)/5 +...]
#
sub atan{
    my $a  = Math::BigInt->new(shift);
    my $x  = Math::BigInt->new(shift);
    my $y  = Math::BigInt->new(0);
    my $X  = Math::BigInt->new($a/$x);
    my $x2 = Math::BigInt->new($x*$x);
    my $d  = Math::BigInt->new(1);
    
    my $temp = Math::BigInt->new(0);
    while($temp = $X/$d){
        $y += $temp;    $d += 2;    $X /= $x2;
        $y -= $X/$d;    $d += 2;    $X /= $x2;
    }
    return $y;
}



#
# y = a * atan(1/x) ただし x=10 を計算する。
#   = a[(1/x^1)/1 -(1/x^3)/3 +(1/x^5)/5 +...]
# 1/x^n を文字列に変換して２回 chop すると 1/x^(n+2) が求まる。
#
sub atan10{
    my $a  = Math::BigInt->new(shift);
    my $x  = Math::BigInt->new(10);
    my $y  = Math::BigInt->new(0);
    my $X  = Math::BigInt->new($a/$x);
    my $x2 = Math::BigInt->new($x*$x);
    my $d  = Math::BigInt->new(1);
    
    my $temp = Math::BigInt->new(0);
    while($temp = $X/$d){
        $y += $temp;    $d += 2;    $X =~ s/..$//;
        $y -= $X/$d;    $d += 2;    $X =~ s/..$//;
    }
    return $y;
}



#
# atan の各項を計算して最後に項の合計を表示する。
#
main:{
    my $scale = shift @ARGV || 100;

    # +/-xxx    => +/-0.00xxx
    # +/-yxxxxx => +/-y.xxxxx
    my $unscale = sub{
        local $_ = shift;
        my $sgn = $&    if s/^[-+]//;
        $_ = ('0' x ($scale-length($_) +1)).$_;
        s/^\d/$&./;
        s/\d\d\d\d\d/$& /g;
        return $sgn.$_;
    };
    
    my @term;
    while(1){
        my $a = shift @param;   last unless($a);
        my $x = shift @param;   last unless($x);

        $VERBOSE && printf("%+4d x atan( 1/+%4d )\t= ",$a,$x);
        my $y = ($x==10) ? atan10($a.('0' x $scale)):
                             atan($a.('0' x $scale),$x);
        push(@term,$y);
        $VERBOSE && print &$unscale($y),"\n";
    }

    my $pi   = Math::BigInt->new(0);
    my $term = Math::BigInt->new(0);
    $pi += $term while($term = pop @term);
    $VERBOSE && print "\t\tpi\t= ",&$unscale($pi),"\n";

    $pi = &$unscale($pi);
    $pi =~ s/(\d{5}\D){10}/\n\t$&/g;
    print "pi = ",$pi,"\n";

}

__END__
--$



使い方は引数として桁数を指定します。省略すると 100 桁。

%perl PerlPi.pl
pi = +3.
        14159 26535 89793 23846 26433 83279 50288 41971 69399 37510
        58209 74944 59230 78164 06286 20899 86280 34825 34211 70683

%

あくまでも 100 桁は内部で計算に使っている桁数なので、その結
果の最後 2 桁くらいはいつも間違ってます。時間を気にしなけれ
ば 32766 桁あたりまでイケたはず。
-- 
FZH01112@..., http://www1.u-netsurf.ne.jp/~dune/