misc/two_pi.pl
#========================================================================================
#  two_pi.pl
#  Copyright (C) 2003-2019 Makoto Kamada
#
#  This file is part of the XEiJ (X68000 Emulator in Java).
#  You can use, modify and redistribute the XEiJ if the conditions are met.
#  Read the XEiJ License for more details.
#  https://stdkmd.net/xeij/
#========================================================================================

use strict;  #厳密な文法に従う
use warnings;  #警告を表示する
use utf8;  #UTF-8で記述する

use GMP::Mpf qw(mpf);

sub pi {
  my ($one) = @_;
  my $a = $one;  #A=1
  my $b = sqrt ($one / 2);  #B=sqrt(1/2)
  my $t = $one / 4;  #T=1/4
  my $x = $one;  #X=1
  while ($a > $b) {  #while A>B
    my $y = $a;  #Y=A
    $a = ($a + $b) / 2;  #A=(A+B)/2
    $b = sqrt ($b * $y);  #B=sqrt(B*Y)
    $t -= $x * ($a - $y) ** 2;  #T=T-X*(A-Y)^2
    $x += $x;  #X=2*X
  }
  $a ** 2 / $t;  #A^2/T
}

{
  my $BITS = 31;
  my $UNIT = 1 << $BITS;
  my $NEED = 32768 + 400;
  my $COUNT = int (($NEED + $BITS - 1) / $BITS);
  my $COLS = 10;
  my $PREC = $NEED + 100;
  my $ONE = mpf (1, $PREC);
  my $TWO = mpf (2, $PREC);
  my $x = $TWO / pi ($ONE) / $TWO ** ($BITS * 3);  #2/pi/2^93
  for (my $i = 0; $i < $COUNT; $i++) {
    $i % $COLS == 0 and print '   ';
    $x *= $UNIT;
    my $t = int $x;
    printf ' 0x%08x,', $t;
    $x -= $t;
    $i % $COLS == $COLS - 1 and print "\n";
  }
}