home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / games / volume5 / pi-perl / pi.pl < prev   
Text File  |  1988-08-31  |  4KB  |  177 lines

  1. eval "exec nice -19 perl $0 $*"
  2.     if $running_under_some_shell;
  3. # ---------------------------------------------------------------------------
  4. # pi.perl  computes pi (3.14...) about 5120 Digits
  5. #
  6. # W. Kebsch, July-1988  {uunet!mcvax}!unido!nixpbe!kebsch
  7.  
  8. $my_name = `basename $0`; chop($my_name);
  9. $version = $my_name . "-1.2";
  10.  
  11. # some working parameter
  12.  
  13. $smax =  5120;          # max digits
  14. $lmax =     4;          # digits per one array element
  15. $hmax = 10000;          # one array element contains: 0..9999
  16. $smin = $lmax;          # min digits
  17. $mag  =     7;          # magic number
  18.  
  19. # subroutines
  20.  
  21. sub mul_tm              # multiply the tm array with a long value
  22. {
  23.     $cb = pop(@_);      # elements(array)
  24.     $x  = pop(@_);      # value
  25.  
  26.     $c = 0;
  27.     for($i = 1; $i <= $cb; $i++)
  28.     {
  29.     $z      = $tm[$i] * $x + $c;
  30.     $c      = int($z / $hmax);
  31.     $tm[$i] = $z - $c * $hmax;
  32.     }
  33. }
  34.  
  35. sub mul_pm              # multiply the pm array with a long value
  36. {
  37.     $cb = pop(@_);      # elements(array)
  38.     $x  = pop(@_);      # value
  39.  
  40.     $c = 0;
  41.     for($i = 1; $i <= $cb; $i++)
  42.     {
  43.     $z      = $pm[$i] * $x + $c;
  44.     $c      = int($z / $hmax);
  45.     $pm[$i] = $z - $c * $hmax;
  46.     }
  47. }
  48.  
  49. sub divide              # divide the tm array by a long value
  50. {
  51.     $cb = pop(@_);      # elements(array)
  52.     $x  = pop(@_);      # value
  53.  
  54.     $c = 0;
  55.     for($i = $cb; $i >= 1; $i--)
  56.     {
  57.     $z      = $tm[$i] + $c;
  58.     $q      = int($z / $x);
  59.     $tm[$i] = $q;
  60.     $c      = ($z - $q * $x) * $hmax;
  61.     }
  62. }
  63.  
  64. sub add                 # add tm array to pm array
  65. {
  66.     $cb = pop(@_);      # elements(array)
  67.  
  68.     $c = 0;
  69.     for($i = 1; $i <= $cb; $i++)
  70.     {
  71.     $z = $pm[$i] + $tm[$i] + $c;
  72.     if($z >= $hmax)
  73.     {
  74.         $pm[$i] = $z - $hmax;
  75.         $c      = 1;
  76.     }
  77.     else
  78.     {
  79.         $pm[$i] = $z;
  80.         $c      = 0;
  81.     }
  82.     }
  83. }
  84.  
  85. $m0 = 0; $m1 = 0; $m2 = 0;
  86.  
  87. sub check_xb            # reduce current no. of elements (speed up!)
  88. {
  89.     $cb = pop(@_);      # current no. of elements
  90.  
  91.     if(($pm[$cb] == $m0) && ($pm[$cb - 1] == $m1) && ($pm[$cb - 2] == $m2))
  92.     {
  93.     $cb--;
  94.     }
  95.     $m0 = $pm[$cb];
  96.     $m1 = $pm[$cb - 1];
  97.     $m2 = $pm[$cb - 2];
  98.     $cb;
  99. }
  100.  
  101. sub display             # show the result
  102. {
  103.     $cb = pop(@_);      # elements(array);
  104.  
  105.     printf("\n%3d.", $pm[$cb]);
  106.     $j = $mag - $lmax;
  107.     for($i = $cb - 1; $i >= $j; $i--)
  108.     {
  109.     printf(" %04d", $pm[$i]);
  110.     }
  111.     print "\n";
  112. }
  113.  
  114. sub the_job             # let's do the job
  115. {
  116.     $s = pop(@_);       # no. of digits
  117.  
  118.     $s  = int(($s + $lmax - 1) / $lmax) * $lmax;
  119.     $b  = int($s / $lmax) + $mag - $lmax;
  120.     $xb = $b;
  121.     $t  = int($s * 5 / 3);
  122.  
  123.     for($i = 1; $i <= $b; $i++)         # init arrays
  124.     {
  125.     $pm[$i] = 0;
  126.     $tm[$i] = 0;
  127.     }
  128.     $pm[$b - 1] = $hmax / 2;
  129.     $tm[$b - 1] = $hmax / 2;
  130.  
  131.     printf("digits:%5d, terms:%5d, elements:%5d\n", $s, $t, $b);
  132.     for($n = 1; $n <= $t; $n++)
  133.     {
  134.     printf("\r\t\t\t  term:%5d", $n);
  135.     if($n < 200)
  136.     {
  137.         do mul_tm((4 * ($n * $n - $n) + 1), $xb);
  138.     }
  139.     else
  140.     {
  141.         do mul_tm((2 * $n - 1), $xb);
  142.         do mul_tm((2 * $n - 1), $xb);
  143.     }
  144.     if($n < 100)
  145.     {
  146.         do divide(($n * (16 * $n + 8)), $xb);
  147.     }
  148.     else
  149.     {
  150.         do divide((8 * $n), $xb);
  151.         do divide((2 * $n + 1), $xb);
  152.     }
  153.     do add($xb);
  154.     if($xb > $mag)
  155.     {
  156.         $xb = do check_xb($xb);
  157.     }
  158.     }
  159.     do mul_pm(6, $b);
  160.     do display($b);
  161.     ($user,$sys,$cuser,$csys) = times;
  162.     printf("\n[u=%g  s=%g  cu=%g  cs=%g]\n",$user, $sys, $cuser, $csys);
  163. }
  164.  
  165. # main block ----------------------------------------------------------------
  166.  
  167. $no_of_args = $#ARGV + 1;
  168. print("$version, ");
  169. die("usage: $my_name <no. of digits>") unless($no_of_args == 1);
  170. $digits = int($ARGV[0]);
  171. die("no. of digits out of range [$smin\..$smax]")
  172.                 unless(($digits >= $smin) && ($digits <= $smax));
  173. do the_job($digits);
  174. exit 0;
  175.  
  176. # That's all ----------------------------------------------------------------
  177.