#!/usr/bin/perl
#
# Compute PI to arbitrary precision with no external lib
# dependency (e.g Math::BigFloat). Just
#
# francesco.prelz@mi.infn.it, November 27th, 2020
#

our @p;
our @t;
our $q;

sub add()
{
  my $j;

  for ($j = $q; $j >= 0; $j--)
   {
    if ($t[$j] + $p[$j] > 9)
    {
      $p[$j] += $t[$j] - 10;
      $p[$j-1] += 1;
    }
    else
     {
      $p[$j] += $t[$j];
     }
   }
}

sub subt()
{
  my $j;

  for ($j = $q; $j >= 0; $j--)
   {
    if ($p[$j] < $t[$j])
     {
      $p[$j] -= $t[$j] -10;
      $p[$j-1] -= 1;
     }
    else
     {
      $p[$j] -= $t[$j];
     }
   }
}

sub mul($)
{
  my $multiplier;
  my $b;
  my $i;
  my $carry;
  my $digit;

  $multiplier = shift;
  $carry = 0;
  for ($i = $q; $i >= 0; $i--)
  {
    $b = ($t[$i] * $multiplier + $carry);
    $digit = $b % 10;
    $carry = int($b / 10);
    $t[$i] = $digit;
  }
}

sub divide($)
{
  my $divisor;
  my $i;
  my $b;
  my $quotient;
  my $remainder;

  $divisor = shift;

  $remainder = 0;
  for ($i = 0; $i <= $q; $i++)
  {
    $b = (10 * $remainder + $t[$i]);
    $quotient = int($b / $divisor);
    $remainder = $b % $divisor;
    $t[$i] = $quotient;
  }
}

sub div4()
{
  my $i;
  my $c;
  my $d;

  $d = 0;
  for ($i = 0; $i <= $q; $i++)
  {
    $c = (10 * $d + $p[$i]) / 4;
    $d = (10 * $d + $p[$i]) % 4;
    $p[$i] = $c;
  }
}

sub mul4()
{
  my $i;
  my $c;
  my $d;

  $c = 0;
  for ($i = $q; $i >= 0; $i--)
  {
    $d = ($p[$i] * 4 + $c) % 10;
    $c = int(($p[$i] * 4 + $c) / 10);
    $p[$i] = $d;
  }
}

sub tiszero()
{
  my $k;

  for ($k = 0; $k <= $q; $k++)
   {
    if ($t[$k] != 0)
     {
      return 0;
     }
   }
  return 1;
}

sub arctan($)
{
  my $s;
  my $n;

  $s = shift;

  $t[0] = 1;
  divide($s);
  add();
  $n = 1;
  do
  {
    mul($n);
    divide($s * $s);
    divide($n += 2);
    if (int(($n-1) / 2) % 2 == 0)
     {
      add();
     }
    else
     {
      subt();
     }
  } while (tiszero() == 0);
}

sub printp()
{
  my $n;

  printf("pi = %d.", $p[0]);
  for ($n = 1; $n < $q; $n++)
   {
    printf("%d", $p[$n]);
   }
  printf("\n");
}

# Check the arguments
if ($#ARGV < 0)
 {
  printf("usage: $0 <precision>\n");
 }
else
 {
  my $n;
  for ($n = 1; $n < $q; $n++)
   {
    $p[$n] = 0;
    $t[$n] = 0;
   }
  $q = $ARGV[0];

  # Compute PI
  arctan(2);
  arctan(3);
  mul4();
  printp();
 }

exit 0;

