Calculate PI and E to Arbitrary Precision

Once upon a time, in a land far away, I wanted to compute pi to several thousand digits. I thought it would be trivial for computers. Not so! Hence, I wrote a Basic program to do it. I have since written it in Pascal (here) and as a Java applet (where?)

Pi = 3.14159 26535 89793 23846 26433
       83279 50288 41971 69399 37510
       58209 74944 59230 78164 06286
       20899 86280 34825 34211 70679
       82148 08651 32823 06647 09384
       46095 50582 23172 53594 08128
       48111 74502 84102 70193 85211
       05559 64462 29489 54930 38196

PI.pas - main program for PI

   1: Program Pi;  {Calculate the number Pi = 3.14159265358979...}
   2:  
   3: uses
   4:   MemTypes, QuickDraw, OSIntf;
   5: 
   6: {$I TIMER.P}              {For elapsed time}
   7: 
   8: Const
   9:   MaxAcc     = 1000;      {Maximum of 1000 digits}
  10:   NBytes     = 501;       {= MaxAcc/2 + 1 extra for rounding}
  11: 
  12: Var
  13:   Accuracy   : Integer;   {Number of digits to report}
  14:   BAcc       : Integer;   {= Accuracy/2 + 1}
  15:   Divisor    : Integer;
  16:   FirstTerm  : Integer;   {First non-zero cell in Term}
  17:   I          : Integer;
  18:   Minus      : Boolean;   {Whether to add or subtract Term}
  19:   Prt        : Text;
  20: 
  21: {$I ARITH.P}              {Extended Arithmetic routines}
  22: 
  23: Begin {======================================================== Main Program}
  24:   ReWrite(Prt, 'PI.PRT');
  25: 
  26:   WriteLn;
  27:   Write('How many digits of accuracy? -->');
  28:   ReadLn(Accuracy);
  29:   If Accuracy > MaxAcc then Begin
  30:     WriteLn('Sorry; maximum accuracy is set to ', MaxAcc, ' digits.');
  31:     Halt;
  32:   End;
  33:   WriteLn(Prt);
  34:   WriteLn(Prt,'Pi to ', Accuracy,' digits of accuracy.');
  35:   BAcc := (Accuracy div 2) + 1;
  36: 
  37: {Print a brief description}
  38: 
  39:   WriteLn(Prt);
  40:   WriteLn(Prt, '  Pi = 16*ArcTan(1/5) - 4*ArcTan(1/239)');
  41:   WriteLn(Prt);
  42:   WriteLn(Prt, '                       1      1       1       1        ');
  43:   WriteLn(Prt, '  Where ArcTan(1/x) = --- - ----- + ----- - ----- + ...');
  44:   WriteLn(Prt, '                       x    3*x^3   5*x^5   7*x^7      ');
  45: 
  46: {Initialize Pi to 1/5 and Term to 16/5}
  47: 
  48:   For I := 0 to BAcc do Begin
  49:     Result[I] := 0;
  50:     Term[I] := 0;
  51:   End;
  52:   Result[1] := 20;  {20/100 = 1/5}
  53:   Term[0] := 3;
  54:   Term[1] := 20;    {320/100 = 16/5}
  55:   FirstTerm := 0;
  56:   Minus := False;
  57:   Divisor := 1;
  58: 
  59: {Calculate next iteration}
  60: 
  61:   StartTimer;
  62:   Repeat
  63:     Minus := Not Minus;    {Signs alternate}
  64:     Multiply(Divisor);
  65:     Divisor := Divisor + 2;
  66:     Divide(Divisor);
  67:     Divide(25);     {5^2 = 25}
  68:     If Minus
  69:       Then Subtract
  70:       Else Add;
  71:   Until FirstTerm >= BAcc;
  72:   WriteLn(Prt);
  73:   WriteLn(Prt, 'Number of terms in ArcTan(1/5) needed was ', Divisor div 2, '.');
  74: 
  75: {Now subtract 4*ArcTan(1/239)}
  76: 
  77:   FirstTerm := 0;
  78:   Term[0] := 4;
  79:   Minus := False;
  80:   Divisor := 1;
  81:   Divide(239);
  82: 
  83: {Divide term by 239 twice, then by divisor}
  84: 
  85:   Repeat
  86:     Minus := Not Minus;    {Signs alternate}
  87:     If Minus
  88:       Then Subtract
  89:       Else Add;
  90:     Multiply(Divisor);
  91:     Divisor := Divisor + 2;
  92:     Divide(239);
  93:     Divide(239);
  94:     Divide(Divisor);
  95:   Until FirstTerm >= BAcc;
  96:   WriteLn(Prt, 'Number of terms in ArcTan(1/239) needed was ', Divisor div 2, '.');
  97: 
  98: {Write out the answer}
  99: 
 100:   WriteLn(Prt);
 101:   WriteLn(Prt, 'Elapsed time = ', ElapsedTime, ' seconds.');
 102:   WriteLn(Prt);
 103:   Write(Prt, ' Pi = 3.');
 104:   PrintResult;
 105:   
 106:   WriteLn;
 107:   Write('Please Press Return');
 108:   ReadLn;
 109:   Close(Prt);
 110: End.
 111: 

E.pas - main program for E

   1: Program E;  {Calculate the number E = 2.718281828459045...}
   2: 
   3: uses
   4:   MemTypes, QuickDraw, OSIntf;
   5: 
   6: {$I TIMER.P}              {For elapsed time}
   7: 
   8: Const
   9:   MaxAcc     = 1000;      {Maximum of 1000 digits}
  10:   NBytes     = 501;       {= MaxAcc/2 + 1 extra for rounding}
  11: 
  12: Var
  13:   Accuracy   : Integer;   {Number of digits to report}
  14:   BAcc       : Integer;   {= Accuracy/2 + 1}
  15:   Divisor    : Integer;
  16:   FirstTerm  : Integer;   {First non-zero cell in Term}
  17:   I          : Integer;
  18:   Prt        : Text;
  19: 
  20: {$I ARITH.P}              {Extended Arithmetic routines}
  21: 
  22: Begin {======================================================== Main Program}
  23:   ReWrite(Prt, 'E.PRT');
  24: 
  25:   WriteLn;
  26:   Write('How many digits of accuracy? -->');
  27:   ReadLn(Accuracy);
  28:   If Accuracy > MaxAcc then Begin
  29:     WriteLn('Sorry; maximum accuracy is set to ', MaxAcc, ' digits.');
  30:     Halt;
  31:   End;
  32:   WriteLn(Prt);
  33:   WriteLn(Prt,'E to ', Accuracy,' digits of accuracy.');
  34:   BAcc := (Accuracy div 2) + 1;
  35: 
  36: {Print a brief description}
  37: 
  38:   WriteLn(Prt);
  39:   WriteLn(Prt, '       1     1     1     1     1       ');
  40:   WriteLn(Prt, '  E = --- + --- + --- + --- + --- + ...');
  41:   WriteLn(Prt, '       0!    1!    2!    3!    4!');
  42: 
  43: {Initialize E to 0 and Term to 1}
  44: 
  45:   For I := 0 to BAcc do Begin
  46:     Result[I] := 0;
  47:     Term[I] := 0;
  48:   End;
  49:   Term[0] := 1;
  50:   FirstTerm := 0;
  51:   Divisor := 1;
  52: 
  53: {Calculate next iteration}
  54: 
  55:   StartTimer;
  56:   Repeat
  57:     Divisor := Divisor + 1;
  58:     Divide(Divisor);
  59:     Add;
  60:   Until FirstTerm >= BAcc;
  61:   WriteLn(Prt);
  62:   WriteLn(Prt, 'Highest Divisor needed was ', Divisor div 2, '.');
  63: 
  64: {Write out the answer}
  65: 
  66:   WriteLn(Prt);
  67:   WriteLn(Prt, 'Elapsed time = ', ElapsedTime, ' seconds.');
  68:   WriteLn(Prt);
  69:   Write(Prt, '  E = 2.');
  70:   PrintResult;
  71:   
  72:   WriteLn;
  73:   Write('Please Press Return');
  74:   ReadLn;
  75:   Close(Prt);
  76: End.
  77: 

timer.p - timer routine

   1: Var
   2:   StartTime    : LongInt;
   3: 
   4: Procedure StartTimer;
   5: Begin {====================================================== StartTimer}
   6:   GetDateTime(StartTime);
   7: End {StartTimer};
   8: 
   9: 
  10: Function ElapsedTime : Integer;
  11: Var
  12:   StopTime     : LongInt;
  13: Begin {===================================================== ElapsedTime}
  14:   GetDateTime(StopTime);
  15:   ElapsedTime := StopTime - StartTime;
  16: End {ElapsedTime};
  17: 

arith.p - extended precision

   1: Type
   2:   Byte           = 0..255;
   3:   BigNumber      = Array[0..NBytes] of Byte;
   4: 
   5: 
   6: Var
   7:   Result         : BigNumber;
   8:   Term           : BigNumber;
   9: 
  10: 
  11: Procedure Add;   {Add Term into Result}
  12: Var
  13:   Carry          : Integer;
  14:   I              : Integer;
  15:   NextNumber     : Integer;
  16:   Pos            : Integer;
  17: Begin {============================================================= Add}
  18:   Carry := 0;
  19:   For I := BAcc downto FirstTerm do Begin
  20:     NextNumber := Carry + Result[I] + Term[I];
  21:     Result[I] := NextNumber mod 100;
  22:     Carry := NextNumber div 100;
  23:   End;
  24:   Pos := FirstTerm;
  25:   While Carry <> 0 do Begin
  26:     Pos := Pos - 1;
  27:     NextNumber := Carry + Result[Pos];
  28:     Result[Pos] := NextNumber mod 100;
  29:     Carry := NextNumber div 100;
  30:   End;
  31: End {Add};
  32: 
  33: 
  34: Procedure Subtract;   {Subtract Term from Result}
  35: Var
  36:   I              : Integer;
  37:   NextNumber     : Integer;
  38:   Pos            : Integer;
  39: Begin {======================================================== Subtract}
  40:   For I := BAcc downto FirstTerm do Begin
  41:     NextNumber := Result[I] - Term[I];
  42:     If NextNumber >= 0 then
  43:       Result[I] := NextNumber
  44:     Else Begin {have to borrow}
  45:       Result[I] := NextNumber + 100;
  46:       Pos := I - 1;
  47:       While Result[Pos] = 0 do Begin
  48:         Result[Pos] := 99;
  49:         Pos := Pos - 1;
  50:       End;
  51:       Result[Pos] := Result[Pos] - 1;
  52:     End;
  53:   End;
  54: End {Subtract};
  55: 
  56: 
  57: Procedure Multiply(K:LongInt);   {Multiply Term by K}
  58: Var
  59:   Carry          : Integer;
  60:   I              : Integer;
  61:   NextNumber     : LongInt;
  62: Begin {======================================================== Multiply}
  63:   Carry := 0;
  64:   For I := BAcc downto FirstTerm do Begin
  65:     NextNumber := K * Term[I] + Carry;
  66:     Term[I] := NextNumber mod 100;
  67:     Carry := NextNumber div 100;
  68:   End;
  69:   If Carry <> 0 then Begin
  70:     FirstTerm := FirstTerm - 1;
  71:     Term[FirstTerm] := Carry;
  72:   End;
  73: End {Multiply};
  74: 
  75: 
  76: Procedure Divide(K:Integer);   {Divide Term by K}
  77: Var
  78:   Hundred        : LongInt;
  79:   I              : Integer;
  80:   NextNumber     : LongInt;
  81:   Remainder      : Integer;
  82: Begin {========================================================== Divide}
  83:   Hundred := 100;
  84:   Remainder := 0;
  85:   For I := FirstTerm to BAcc do Begin
  86:     NextNumber := Hundred * Remainder + Term[I];
  87:     Term[I] := NextNumber div K;
  88:     Remainder := NextNumber mod K;
  89:   End;
  90:   While Term[FirstTerm] = 0 do
  91:     FirstTerm := FirstTerm + 1;
  92: End {Divide};
  93: 
  94: 
  95: Procedure PrintResult;  {Print the result, nicely}
  96: Var
  97:   I              : Integer;
  98:   Pos            : Integer;
  99: Begin {===================================================== PrintResult}
 100:   Pos := 1;
 101:   For I := 1 to Accuracy do Begin
 102:     If Odd(I) then
 103:       Write(Prt, Result[Pos] div 10)
 104:     Else Begin
 105:       Write(Prt, Result[Pos] mod 10);
 106:       Pos := Pos + 1;
 107:     End;
 108:     If I mod 5 = 0 then Write(Prt, ' ');
 109:     If I mod 50 = 0 then Begin
 110:       WriteLn(Prt);
 111:       Write(Prt, '        ');
 112:     End;
 113:   End;
 114:   If Accuracy mod 50 <> 0 then WriteLn(Prt);
 115: End {PrintResult};
 116: 

PIE.out - sample output

   1: Calculating Pi = 3.14159 ...
   2: 
   3: Pi to 500 digits of accuracy.
   4: 
   5:   Pi = 16*ArcTan(1/5) - 4*ArcTan(1/239)
   6: 
   7:                        1      1       1       1        
   8:   Where ArcTan(1/x) = --- - ----- + ----- - ----- + ...
   9:                        x    3*x^3   5*x^5   7*x^7      
  10: 
  11: Number of terms in ArcTan(1/5) needed was 356.
  12: Number of terms in ArcTan(1/239) needed was 105.
  13: 
  14: Elapsed time = 135 seconds.
  15: 
  16:  Pi = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 
  17:         58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 
  18:         82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 
  19:         48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 
  20:         44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 
  21:         45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 
  22:         72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 
  23:         78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 
  24:         33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 
  25:         07446 23799 62749 56735 18857 52724 89122 79381 83011 94913 
  26: 
  27: 
  28: Calculating E = 2.71828 ...
  29: 
  30: E to 500 digits of accuracy.
  31: 
  32:        1     1     1     1     1       
  33:   E = --- + --- + --- + --- + --- + ...
  34:        0!    1!    2!    3!    4!
  35: 
  36: Highest Divisor needed was 127.
  37: 
  38: Elapsed time = 28 seconds.
  39: 
  40:   E = 2.71828 18284 59045 23536 02874 71352 66249 77572 47093 69995 
  41:         95749 66967 62772 40766 30353 54759 45713 82178 52516 64274 
  42:         27466 39193 20030 59921 81741 35966 29043 57290 03342 95260 
  43:         59563 07381 32328 62794 34907 63233 82988 07531 95251 01901 
  44:         15738 34187 93070 21540 89149 93488 41675 09244 76146 06680 
  45:         82264 80016 84774 11853 74234 54424 37107 53907 77449 92069 
  46:         55170 27618 38606 26133 13845 83000 75204 49338 26560 29760 
  47:         67371 13200 70932 87091 27443 74704 72306 96977 20931 01416 
  48:         92836 81902 55151 08657 46377 21112 52389 78442 50569 53696 
  49:         77078 54499 69967 94686 44549 05987 93163 68892 30098 79311 
  50: 
  51: 
  52: Steven A. O'Hara
  53: October 1985
  54: 
Email: steve@oharasteve.com