view release on metacpan or  search on metacpan

t/test_pi.pl  view on Meta::CPAN

use Acme::Tools;
print "$Acme::Tools::PI\n";
my $pi=bigf($PI);
my $pi_big=bigf(repl(<<'',qr/\s/));
      3.1415926535897932384626433832795028841971693993751058209749445923078
        164062862089986280348253421170679821480865132823066470938446095505822317
        253594081284811174502841027019385211055596446229489549303819644288109756
        659334461284756482337867831652712019091456485669234603486104543266482133
        936072602491412737245870066063155881748815209209628292540917153643678925
        903600113305305488204665213841469519415116094330572703657595919530921861
        173819326117931051185480744623799627495673518857527248912279381830119491
        298336733624406566430860213949463952247371907021798609437027705392171762
        931767523846748184676694051320005681271452635608277857713427577896091736
        371787214684409012249534301465495853710507922796892589235420199561121290
        219608640344181598136297747713099605187072113499999983729780499510597317
        328160963185950244594553469083026425223082533446850352619311881710100031
        378387528865875332083814206171776691473035982534904287554687311595628638
        823537875937519577818577805321712268066130019278766111959092164201989380
        952572010654858632788659361533818279682303019520353018529689957736225994
        138912497217752834791315155748572424541506959508295331168617278558890750
        983817546374649393192550604009277016711390098488240128583616035637076601
        047101819429555961989467678374494482553797747268471040475346462080466842
        590694912933136770289891521047521620569660240580381501935112533824300355
        876402474964732639141992726042699227967823547816360093417216412199245863
        150302861829745557067498385054945885869269956909272107975093029553211653
        449872027559602364806654991198818347977535663698074265425278625518184175
        746728909777727938000816470600161452491921732172147723501414419735685481


pi_bin();

sub pi_1 { #  pi = 4 sigma(0..inf) -1^k/(2k+1)
  for my $n (map 10**$_,1..18){
    my($start,$sum,$one,$e)=(time_fp(),0,bigf(-1),0);
    $sum+=($one*=-1)/(2*$_+1) for 0..$n;
    my $mypi=4*$sum;
    printf "%7d: ".("%30.25f" x 5)."  %5.2fs\n",
      $n,
      $mypi,
      $pi-$mypi,
      $pi-($mypi - 1/$n**1),
      $pi-($mypi - 1/$n**1 + 1/$n**2),
      $pi-($mypi - 1/$n**1 + 1/$n**2 - 0.75/$n**3),
      time_fp()-$start;
  }
}

sub pi_2 { # pi^2/6 = 1/1**2 + 1/2**2 + 1/3**2 + 1/4**2 ...
  for my $n (map 10**$_,0..8){
    my($start,$sum)=(time_fp(),0);
    $sum+=1/$_**2 for 1..$n;
    my $mypi=sqrt(6*$sum);
    printf "%9d: ".("%30.25f" x 2)."  %5.2fs\n",  $n, $mypi, $pi-$mypi, time_fp()-$start;
  }
}

sub pi_3 { # dart and pythagoras
    for my $n (map 10**$_,0..8){
	my($start,$s)=(time_fp(),0);
	for(1..$n){
	    my($x,$y)=(rand(),rand()); #throw dart
	    ++$s if sqrt($x*$x + $y*$y) < 1;
	}
	my $mypi=4*$s/$n;
	printf "%9d: %30.25f  %30.25f  %5.2fs\n", $n, $mypi, $pi-$mypi, time_fp()-$start;
    }
}

#use Math::BigFloat lib=>"GMP";# if !$INC{'Math/BigFloat.pm'};
sub pi_4 { # ramaputramama...
    #use Math::BigFloat ':constant';
    my @fak; $fak[$_]=$_?$fak[$_-1]*$_:bigf(1) for 0..1000; #die join("\n",@fak)."\n";
    bigscale(1000); #hm
    my $pi_bigger=Math::BigFloat->bpi(1000);
    for my $n (30..50){
	my($start,$s)=(time_fp(),bigf(0));
	for my $k (0..$n) {
	    my $kf=bigf($k);
	    $s+=  $fak[$k*4] / $fak[$k]**4
		* (1103 + 26390*$kf) / 396**($kf*4)
	}
	$s*=2*sqrt(bigf(2))/9801;