Contributor: DAVID ADAMSON
(* PROGRAM PI5.PAS
By David Adamson
October 23, 1994
This it the third version of a program for computing Pi up to 2150
places.
The first version was written in Turbo Pascal 3. The second was also
in Turbo Pascal, but with all the math procedures converted to inline
for maximum program executing speed.
This version is rewritten to allow compilation with Turbo Pascal 5 and
up.
The inline code was scrapped for version 3. The sacrifice in speed is
not a problem, because it improvement in today's computers overshadows
the processing efficiency loss.
Here are my results using a 33MHz 386.
Pi v2 Pi v3
100 places = 0.1 sec 100 places = 0.3 sec
1000 places = 11 sec 1000 places = 32 sec
2135 places = 48 sec 2135 places = 150 sec
----------------------------------------------------------------------
THIS IS THE ORIGINAL HEADER:
PROGRAM PI1.PAS
TANDY 2000
BY CINO HILLIARD
APRIL 6 , 1986
( AN IMPROVED VERSION OF MY ORIGIONAL PI.PAS IN THE DL)
THIS PROGRAM COMPUTES THE DIGITS OF PI USING THE ARCTANGENT FORMULA
( TERM 1 ) ( TERM 2 )
(1) PI/4 = 4 ARCTAN 1/5 - ARCTAN 1/239
IN CONJUNCTION WITH THE GREGORY SERIES
N
(2) ARCTAN X = SUM [ (-1)^M*(2M + 1)^-1*X^(2M+1) ] (APPROXIMATELY)
M=0
SUBSTITUTING INTO (1) AND (2) THE FIRST FEW VALUES OF N, REARRANGING,
SIMPLYIFING, AND NESTING WE HAVE,
PI= 3.2 + 1/25[-3.2/3 + 1/25[3.2/5 + 1/25[-3.2/7 + ...].].] ( TERM 1 )
-1/239[4 +1/239^2[-4/3 +1/239^2[4/5 +1/239^2[-4/7 +...].].] ( TERM 2 )
USING THE LONG DIVISION ALGORITHM AND SOME TRICKS I DISCOVERED, THIS (
NESTED ) INFINITE SERIES CAN BE USED TO CALCULATE PI TO A LARGE NUMBER
OF DECIMAL PLACES IN A REASONABLE AMOUNT OF TIME. A TIME FUNCTION IS
INCLUDED TO SHOW HOW SLOW THINGS GET WHEN N IS LARGE.
IMPROVEMENTS CAN BE MADE BY INCREASING THE NUMBER OF DIGITS EACH ARRAY
ELEMENT HOLDS AND BY CHANGING THE DATA TYPE FROM INTEGER TO REAL FOR
SELECTED VARIABLES. OF COURSE THE ADDED NUMBER OF DIGITS THESE CHANGES
PRODUCE WILL COST MUCH MUCH MORE TIME. AH INDEED, 'TIS NO FREE LUNCH!
HOWEVER, SINCE TERM 1 AND TERM 2 ARE COMPUTED SEPERATELY AND SINCE THE
ARRAYS ARE STEP BY STEP UPDATED, THE PROGRAM DOES LEND ITSELF TO
PARALLEL OR NON VON ( WHAT A COINCIDENCE - SEE BELOW ) PROCESSING.
FOR EXAMPLE, LET COMPUTER 1 PERFORM TERM 1 AND COMPUTER 2 PERFORM
TERM 2. MOREOVER, LET SEVERAL COMPUTERS SHARE IN THE CALCULATION OF
EACH OF THE INDIVIDUAL TERMS. HOWEVER, TO KEEP EACH COMPUTER EQUALLY
BUSY, A LOGARITHMIC TYPE OF ADJUSTMENT MUST BE MADE TO DECIDE ON THE
NUMBER OF TERMS TO BE ASSIGNED TO EACH COMPUTER. SINCE THE HIGHER
POWER TERMS REQUIRE MUCH LONGER TO COMPUTE, THE COMPUTERS ASSIGNED TO
HIGHER POWERS MUST BE GIVEN FEWER TERMS TO DO.
A LITTLE HISTORY
----------------
IN AUGUST, 1949, PROFESSOR JOHN VON NEUMANN USED FORMULAS (1) AND (2)
TO CALCULATE PI TO 2035 DECIMAL PLACES ON THE ENIAC COMPUTER. THE
EFFORT WAS MADE TO DETERMINE IF THE DIGITS CONFORMED TO SOME TYPE OF
PATTERN OR IF THEY WERE RANDOMLY DISTRIBUTED. THE CALCULATION WAS
COMPLETED OVER THE LABOR DAY WEEKEND WITH THE COMBINED EFFORTS OF FOUR
ENIAC STAFF MEMBERS WORKING IN EIGHT-HOUR SHIFTS TO ENSURE CONTINUOUS
OPERATION OF THE ENIAC. THE CALCULATION (INCLUDING CARD HANDLING TIME)
TOOK APPROXIMATELY 70 HOUR.
THE CONCLUSION WAS AS SUSPECTED - THE DIGITS APPEARED TO BE RANDOM!
SOME YEARS AGO I REQUESTED INFORMATION ON PI FROM THE ENCYCLOPEDIA
BRITANNICA RESEARCH SERVICE. I RECEIVED A REPORT GIVING THE ABOVE
HISTORICAL ACCOUNT PLUS A LISTING OF THE 2035 DIGITS.
A LITTLE COMMENT
----------------
USING MY T2K, PI1.PAS COMPUTES PI TO 2035 PLACES IN 15 MIN 31.97
SEC.
GEE WHIZ! WHAT EARTHLY USE CAN BE MADE OF THIS? I PLAN TO MAKE
DESIGNS AND COLOR PATTERNS USING THE DIGIT VALUES AS THE BASIS. MAYBE
SOMETHING LIKE FRACTALS WOULD BE A GOOD START. THIS IS WHY COMPUTATION
AND SPEED ARE IMPORTANT AS I DO NOT WISH (NOR TRUST) TO KEY IN ENTRIES
FROM A PUBLISHED LIST.
CINO HILLIARD
[72756,672]
---------------------------------------------------------------------------
MODIFICATIONS MADE FOR VERSION 2:
ADDED INLINE CODE TO IMPROVE THE SPEED. USING THE ORIGINAL SOURCE ON
MY TANDY-1000 ( WITH A V20 CPU ) , COMPUTING PI TO 2035 PLACES RAN IN
39 MIN 44.92 SEC.. AFTER ADDING THE INLINE CODE IT RAN IN 10 MIN 25.16
SEC..
ADDED A TIME FUNCTION FROM 'TECH JOURNAL FEB 85'. THIS USES A MSDOS
INTERUPT, SO IT SHOULD RUN ON ANY MSDOS COMPATIBLE COMPUTER.
CHUCK WHITE
[75006,3677]
*)
PROGRAM PI5;
Uses DOS, CRT;
Type
TimeString = string[12];
VAR K,I,I2,J,M,N,Q,V,R,D,Z : INTEGER;
A,P,T : ARRAY[0..5000] OF integer;
TI,T2 : STRING[20];
function Time: TimeString;
var
Hour, Minute, Second, Sec100 : word;
Hr, Min, Sec, Hun : string[2];
begin
gettime(Hour, Minute, Second, Sec100);
str(Hour:2, Hr);
str(Minute:2, Min);
str(Second:2, Sec);
str(Sec100:2, Hun);
if Hr[1] = ' ' then Hr[1] := '0';
if Min[1] = ' ' then Min[1] := '0';
if Sec[1] = ' ' then Sec[1] := '0';
if Hun[1] = ' ' then Hun[1] := '0';
time := Hr+ ':'+ Min+ ':'+ Sec+ '.'+ Hun
end;
PROCEDURE DIV32IA; { DIVIDE 3.2 BY I AND STORE IN ARRAY A }
BEGIN
Q:=3 DIV I;
R:=3 MOD I;
A[0]:=Q;
V:=R*10+2;
Q:=V DIV I;
R:=V MOD I;
A[1]:=Q;
FOR J:=2 TO M DO
BEGIN
V:=R*10;
Q:=V DIV I;
R:=V MOD I;
A[J]:=Q;
END;
END;
PROCEDURE DIVA(D:INTEGER); { DIVIDE A BY SPECIFIED INTEGER FROM }
BEGIN { PROCEDURE COMPUTE; }
R:=0;
FOR J:=0 TO M DO
BEGIN
V:= R*10+A[J];
Q:= V DIV D;
R:= V MOD D;
A[J]:=Q;
END;
END;
PROCEDURE DIV4IA; { DIVIDE 4 BY I AND STORE IN A }
BEGIN
Q:=4 DIV I;
R:=4 MOD I;
A[0]:=Q;
FOR J:=1 TO M DO
BEGIN
V:= R*10;
Q:= V DIV I;
R:= V MOD I;
A[J]:=Q;
END;
END;
PROCEDURE DIV32IP; { DIVIDE 3.2 BY I AND STORE IN P }
BEGIN
Q:=3 DIV I;
R:=3 MOD I;
P[0]:=Q;
V:=R*10+2;
Q:=V DIV I;
R:=V MOD I;
P[1]:=Q;
FOR J:=2 TO M DO
BEGIN
V:= R*10;
Q:= V DIV I;
R:= V MOD I;
P[J]:=Q;
END;
END;
PROCEDURE DIV4IP; { DIVIDE 4 BY I AND STORE IN P }
BEGIN
Q:=4 DIV I;
R:=4 MOD I;
P[0]:=Q;
FOR J:=1 TO M DO
BEGIN
V:= R*10;
Q:= V DIV I;
R:= V MOD I;
P[J]:=Q;
END;
END;
PROCEDURE SUBA; { SUBTRACT A FROM P AND STORE IN A }
BEGIN
FOR J:=0 TO M DO
A[J]:=P[J]-A[J];
END;
PROCEDURE SUB32A; { SUBTRACT A FROM 3.2 AND STORE IN T }
BEGIN
T[0]:=3;
T[1]:=1;
FOR J:=2 TO M DO
T[J]:=9-A[J];
END;
PROCEDURE SUB4A; { SUBTRACT A FROM 4 AND STORE IN A }
BEGIN
A[0]:=3;
FOR J:=1 TO M DO
A[J]:=9-A[J];
END;
PROCEDURE SUBT; { SUBTRACT TERM2 FROM TERM 1 AND STORE IN A }
BEGIN
FOR J:=M DOWNTO 1 DO
BEGIN
A[J]:=T[J]-A[J];
WHILE A[J]<0 DO
BEGIN
A[J]:=A[J]+10; { ADJUST FOR NEGATIVE }
A[J-1]:=A[J-1]+1;
END;
WHILE A[J]>9 DO { ADJUST A FOR EXCESS CARRY }
BEGIN
A[J]:=A[J]-10;
A[J-1]:=A[J-1]-1;
END;
END;
END;
PROCEDURE COMPUTE; { COMPUTE THE NESTED SERIES FOR I2 ITERATIONS }
BEGIN
I:=I2; { COMPUTE TERM 1 TO I2 PLACES }
DIV32IA;
DIVA(25);
WHILE I>3 DO
BEGIN
I:=I-2;
DIV32IP;
SUBA;
DIVA(25);
END;
SUB32A;
I:=I2; { COMPUTE TERM 2 TO I2 PLACES }
DIV4IA;
DIVA(239);
DIVA(239);
WHILE I>3 DO
BEGIN
I:=I-2;
DIV4IP;
SUBA;
DIVA(239);
DIVA(239);
END;
SUB4A;
DIVA(239);
SUBT; { SUBTRACT TERM 2 FROM TERM 1 AND STORE IN A }
T2:=TIME; { SET END OF COMPUTATION TIME }
END;
PROCEDURE STORE; { SAVE PI AND TIME TO DISK FILES }
VAR BUFF : FILE OF INTEGER;
F : TEXT;
BEGIN
ASSIGN(BUFF,'PI.DTA');
ASSIGN(F,'TIME.DTA');
REWRITE(BUFF);
REWRITE(F);
FOR J:=0 TO M-2 DO
WRITE(BUFF,A[J]);
WRITE(F,TI); {Note, does not produce an ascii readable}
WRITE(F,T2); {output. It is for use by PILIST - code }
{for PILIST is at the end after PI5}
CLOSE(BUFF);
CLOSE(F);
END;
procedure compute_time;
var
h1,h2,m1,m2,s1,s2,hun1,hun2,code,x: integer;
s,temp: string[11];
begin
val(copy(ti,1,2),h1,code);
val(copy(t2,1,2),h2,code);
val(copy(ti,4,2),m1,code);
val(copy(t2,4,2),m2,code);
val(copy(ti,7,2),s1,code);
val(copy(t2,7,2),s2,code);
val(copy(ti,10,2),hun1,code);
val(copy(t2,10,2),hun2,code);
if hun2