view release on metacpan or  search on metacpan

scripts/boxmuller  view on Meta::CPAN

#!/usr/bin/perl
use 5.014 ; use strict ; use warnings ;  # the functions requires 5.10 for "state", 5.14 for srand. 
use Getopt::Std ; getopts ':.:@:1d:g:Lm:s:v:', \my%o ;  
use Math::Trig qw/pi/ ; # 5.4から
use Scalar::Util qw/looks_like_number/ ; # 5.7.3から
use Term::ANSIColor qw/:constants color/ ;  $Term::ANSIColor::AUTORESET = 1 ;
use Time::HiRes qw/sleep usleep gettimeofday tv_interval/ ; # 5.7.3から

$SIG{INT} = sub { & SecondInfo ; exit 130 } ;

my $time0 = [ gettimeofday ] ;
my ( $mu , $sd ) ;  #mu : 平均 , sd : 標準偏差
my ( $s1 , $s2 )  = (0,0) ; # 1乗和 と 2乗和
my $count = 0 ; # 出力した個数
my $upto = $o{g} // 6 ;  # 出力要素数
& init ; 
& main ; 
& SecondInfo ; 
exit 0 ;

sub init ( ) {   #オプションを使った設定
   $o{s} = defined $o{s} ? srand $o{s} : srand ; # 乱数シードの保管/設定

   sub LLN ( $ ) ; * LLN = * looks_like_number ; # 関数名が長すぎるので、短くした。
   sub printErr( $ ){ print STDERR BRIGHT_RED "Option -$_[0] should have a numeric specification.\n" ; exit 1 }
   $mu = $o{m} ? LLN $o{m} ? $o{m} : printErr "m" : 0 ;  #m : 平均
   $sd = $o{d} ? LLN $o{d} ? $o{d} : printErr "d" : $o{v} ? LLN $o{v} ? sqrt  $o{v} : printErr "v" : 1 ;  #sd：標準偏差 
}

sub main ( ) {  #  乱数の出力
   sub getrand  ;
   sub boxmuller ( $$ ) ;  #ボックスミュラー法によるガウス乱数の作成
   * getrand = * boxmuller ; 
   * getrand = * lognormal if $o{L} ;  # 対数正規分布の指定があった場合。

   while ( $count < $upto ) { 
   	   sleep $o{'@'} if defined $o{'@'} ;
       my $x = getrand $mu, $sd ; 
       $x = sprintf "% .$o{'.'}f" , $x if defined $o{'.'} ; # <-- May be efficientized. 
       $s1 += $x ; 
       $s2 += $x ** 2 ; 
       $count ++ ; # 出力個数を計数するので個々を選んだ。
       print "$count\t" if $o{':'} ;  # <-- Maybe effiecientized by other code structure.
       print "$x\n" ; 
   }

  sub boxmuller ( $$ ) {  #  ガウス乱数の生成
      state $piW = atan2(1,1)* 8  ;  #6.28=3.14*2
      state $z = undef ; 
      do { my $t = $z ; undef $z ; return $t * $_[1] + $_[0] }  if defined $z ; 
      
      my $r1 = 1 - rand ; # ; $r1 = rand until $r1 ; 値が0になることを阻止。
      my $r2 = $piW * rand ;
      my $r1R = sqrt ( -2 * log $r1 ); 
      my $r2S = sin $r2 ; 
      my $r2C = cos $r2 ; 
      ( my $t , $z ) = ( $r1R * $r2S , $r1R * $r2C ) ; 
      return $t * $_[1] + $_[0] ; 
  }

  sub lognormal ( $$ ) { 
  	return exp boxmuller $_[0], $_[1] ;
  }
}

sub SecondInfo( ) {   #  処理したことについての二次情報を出力
    return if $o{1} ;
    use FindBin qw [ $Script ] ; 
    my $cmd = "$Script -m $mu -d $sd" ; 
    $cmd .= ' -L' if $o{L} ;
    print STDERR 
       CYAN "printed lines: ", BRIGHT_CYAN $count ,
       CYAN " , used random seed: " , BRIGHT_CYAN  $o{s} ,
       CYAN " , elapsed seconds: " , BRIGHT_CYAN  tv_interval ($time0) ,
       RESET "\n" , 
       CYAN "sum = " , BRIGHT_CYAN  sprintf("%g", $s1 ) ,
       CYAN " , squared sum = " , BRIGHT_CYAN  sprintf( "%g" , $s2 ) , 
       CYAN " ($cmd) " , RESET "\n" ;
 }

## ヘルプとバージョン情報
BEGIN {
  our $VERSION = 0.24 ;
  $Getopt::Std::STANDARD_HELP_VERSION = 1 ;