begin
comment ****************************************************************
*** ***
*** This calculation of PI uses the Bailey-Borwein-Plouffe ***
*** algorithm, whereby PI = Sum [K=0 to Inf] (1/16)**K ***
*** *( 4/(8*K+1) -2/(8*K+4) -1/(8K+5) -1/(8*K+6) ). ***
*** At each iteration, the four arrays BFrac1 to BFrac4 are ***
*** populated with the binary expansions of the fractional ***
*** terms, then these are combined, and finally shifted by an***
*** amount corresponding to the factor (1/16)**K, upon which ***
*** this is added to the accumulating total. Hence we get ***
*** an answer of gradually increasing precision, viewable ***
*** by setting bit 2 of the console word generator. ***
*** ***
*** Each iteration increases the precision by slightly ***
*** better than 4 binary bits, so we can estimate the decimal***
*** precision correspondingly at each iteration, or otherwise***
*** by simply comparing results of consecutive iterations. ***
*** ***
*** Diagnostics Setting BLine causes binary dump of some ***
*** =========== arrays - to see all remove the Comments ***
*** Setting One on the word generator causes a ***
*** trace of significant operations ***
*** Setting Two on the word generator shows the ***
*** decimal accumulation at EVERY iteration ***
*** ***
*** Author: Bob Firth ***
*** Written: June,July 2010 ***
*** Contact: General_Factotum@hotmail.co.uk ***
*** ***
*** Comment: The machine code of the 803/503 is not ***
*** particularly suited to multi-word precision ***
*** arithmetic, mainly because of the absence of ***
*** double-length add and subtract. Nevertheless ***
*** it is possible to program around this by ***
*** detecting and handling overflow while processing***
*** the lower order word. ***
*** ***
****************************************************************;
procedure PWRDS(Ident,AStart,ASize);
value AStart,ASize;
integer AStart,ASize;
string Ident;
comment *** Dumps binary representation of elements of array ***
*** starting at AStart, length 1+ASize ***
****************************************************************;
begin integer BLine,NoBits,NoCHs,Count,Word;
switch S ≔ PACE,LOOP,NEXT,EXIT,GETOUT;
BLine ≔ 24288;
Elliott(7,0,0 ,0,0,3,BLine );
Elliott(4,2,GETOUT,0,0,0,0 );
for Count ≔ step 1 until ASize do begin
Elliott(7,0,0 ,0,0,3,BLine );
Elliott(4,2,EXIT ,0,0,0,0 );
PRINT(“\n”,Ident,sameline,digits(3),Count,“) ”);
NOBITS ≔ 39;
Elliott(3,0,AStart,0,0,4,Count);
Elliott(2,0,4 ,1,3,0,0 );
Elliott(2,0,Word ,0,0,0,0 );
SPACE: PRINT(“ ”);
NOCHS ≔ 2;
LOOP: ELLIOTT(3,2,NoBits,0,4,2,EXIT );
if WORD < 0 then PRINT(“1”);
if WORD ⩾ 0 then PRINT(“0”);
ELLIOTT(3,0,Word ,0,5,5,1 );
ELLIOTT(2,0,Word ,0,4,3,NEXT );
NEXT: ELLIOTT(3,2,NoCHs ,0,4,2,SPACE );
ELLIOTT(4,0,LOOP ,0,0,0,0 );
EXIT: end;
GETOUT:
end PWRDS;
comment ================================================================;
procedure INTDV(DVDend,DIVsor,AStart,ASize);
value DVDend,DIVsor,AStart,ASize;
integer DVDend,DIVsor,AStart,ASize;
comment *** Divides integer word DVDend by integer word DIVsor, ***
*** result is stored in array ARRAY[0:ASize] ***
*** where whole part of quotient goes to ARRAY[0] ***
*** and fractional part populates rest of array. ***
*** Array position and length given by AStart and ASize+1 ***
****************************************************************;
begin integer Maxpos,One,Remain,Count,Point;
switch S ≔ OTRAC,LOOP;
Maxpos ≔ 74877906943;
One ≔ ;
Elliott(7,0,0 ,0,0,3,One );
Elliott(4,2,NOTRAC,0,0,0,0 );
PRINT(“Dividing ”,SAMELINE,DIGITS(1),DVDend,“ by ”,DIGITS(5),DIVsor);
NOTRAC: Elliott(2,6,Remain,0,2,6,Count );
Elliott(3,0,AStart,0,2,0,Point );
Elliott(3,0,DVDend,0,2,0,4 );
Elliott(5,0,38 ,0,0,0,0 );
LOOP: Elliott(3,0,Remain,0,5,6,DIVsor);
Elliott(0,0,Point ,1,2,0,0 );
Elliott(5,2,DIVsor,0,5,7,0 );
Elliott(0,7,4 ,0,0,3,Maxpos);
Elliott(2,0,Remain,0,2,6,4 );
Elliott(2,2,Point ,0,0,0,0 );
Elliott(3,2,Count ,0,0,5,ASize );
Elliott(5,1,0 ,0,4,1,LOOP );
end INTDV;
comment ================================================================;
procedure COMBI(T1Strt,T2Strt,T3Strt,T4Strt,TTStrt,Max);
value T1Strt,T2Strt,T3Strt,T4Strt,TTStrt,Max;
integer T1Strt,T2Strt,T3Strt,T4Strt,TTStrt,Max;
comment *** Take (binary fraction at) T1Strt, subtract that at T2Strt***
*** and at T3Strt and at T4Strt, and put result at TTStrt ***
****************************************************************;
begin integer Maxpos,One,Count,Carry;
switch S ≔ OOP,UFLO1,SKIP1,UFLO2,SKIP2,UFLO3,SKIP3,UFLO4,SKIP4;
Maxpos ≔ 74877906943;
One ≔ ;
T1Strt ≔ 1Strt+Max;
T2Strt ≔ 2Strt+Max;
T3Strt ≔ 3Strt+Max;
T4Strt ≔ 4Strt+Max;
TTStrt ≔ TStrt+Max;
Elliott(3,0,Max ,0,2,0,Count );
Elliott(2,6,Carry ,0,2,6,4 );
LOOP: Elliott(0,0,T1Strt,1,3,0,0 );
Elliott(0,0,T2Strt,1,0,5,0 );
Elliott(4,1,UFLO1 ,0,4,0,SKIP1 );
UFLO1: Elliott(0,3,MaxPos,0,2,2,4 );
SKIP1: Elliott(0,0,T3Strt,1,0,5,0 );
Elliott(4,1,UFLO2 ,0,4,0,SKIP2 );
UFLO2: Elliott(0,3,MaxPos,0,2,2,4 );
SKIP2: Elliott(0,0,T4Strt,1,0,5,0 );
Elliott(4,1,UFLO3 ,0,4,0,SKIP3 );
UFLO3: Elliott(0,3,MaxPos,0,2,2,4 );
SKIP3: Elliott(0,0,0 ,0,0,5,Carry );
Elliott(4,1,UFLO4 ,0,4,0,SKIP4 );
UFLO4: Elliott(0,3,MaxPos,0,2,2,4 );
SKIP4: Elliott(0,0,TTStrt,1,2,0,0 );
Elliott(3,6,4 ,0,2,0,Carry );
Elliott(3,0,One ,0,2,7,T1Strt);
Elliott(2,7,T2Strt,0,2,7,T3Strt);
Elliott(2,7,T4Strt,0,2,7,TTStrt);
Elliott(0,0,0 ,0,3,7,Count );
Elliott(0,1,0 ,0,4,1,LOOP );
end COMBI;
comment ================================================================;
procedure SHIFT(Expone,TTStrt,SHStrt,Max);
value Expone,TTStrt,SHStrt,Max;
integer Expone,TTStrt,SHStrt,Max;
comment *** Shift contents of array at TTStrt down by Expone*4 bits ***
*** and place result at array SHStrt, both length Max ***
****************************************************************;
begin integer Maxpos,One,Count,ShiftB,ShiftL,ShiftR,ShiftW,Source;
switch S ≔ ONE,NOTRAC,MORE,SHIF1,NONE1,NEXT1,SHIF2,NONE2,NEXT2,NOSHIF,EXIT;
Maxpos ≔ 74877906943;
One ≔ ;
Count ≔ xpone ÷ 19;
ShiftB ≔ Expone - Count×19 )×4;
ShiftW ≔ ount ×2;
ShiftL ≔ hiftR ≔ ource ≔ ;
if ShiftB ⩾ 60 then begin
ShiftW ≔ hiftW +1;
ShiftL ≔ 6 - ShiftB;
goto DONE;
end;
if ShiftB ⩾ 40 then begin
ShiftW ≔ hiftW +1;
ShiftR ≔ hiftB - 38;
Source ≔ ;
goto DONE;
end;
if ShiftB ⩾ 20 then begin
ShiftL ≔ 8 - ShiftB;
goto DONE;
end;
if ShiftB ⩾ 0 then begin
ShiftR ≔ hiftB;
Source ≔ ;
end;
DONE: Elliott(7,0,0 ,0,0,3,One );
Elliott(4,2,NOTRAC,0,0,0,0 );
PRINT(“\nShift Parameters”,SAMELINE,DIGITS(3),Count,ShiftB,ShiftL,ShiftR,ShiftW,Source);
NOTRAC: Elliott(0,0,0 ,0,2,6,Count );
MORE: Elliott(3,0,Count ,0,0,5,ShiftW);
Elliott(4,1,NONE1 ,0,0,4,TTStrt);
Elliott(2,0,4 ,1,3,0,0 );
SHIF1: Elliott(5,0,38 ,0,4,0,NEXT1 );
NONE1: Elliott(0,6,0 ,0,4,0,SHIF1 );
NEXT1: Elliott(3,0,Count ,0,0,5,ShiftW);
Elliott(0,5,One ,0,4,1,NONE2 );
Elliott(0,4,TTStrt,0,0,0,0 );
Elliott(2,0,4 ,1,3,0,0 );
SHIF2: Elliott(4,0,NEXT2 ,0,0,0,0 );
NONE2: Elliott(0,6,0 ,0,0,0,0 );
NEXT2: Elliott(0,0,ShiftL,1,5,4,0 );
Elliott(0,0,ShiftR,1,5,0,0 );
Elliott(0,3,Maxpos,0,2,0,4 );
Elliott(3,0,Source,0,4,2,NOSHIF);
Elliott(5,7,0 ,0,2,0,4 );
NOSHIF: Elliott(3,0,4 ,0,0,0,0 );
Elliott(0,0,SHStrt,1,2,0,0 );
Elliott(2,2,SHStrt,0,3,0,Count );
Elliott(0,4,One ,0,2,0,Count );
Elliott(0,5,Max ,0,0,5,One );
Elliott(4,1,MORE ,0,4,3,EXIT );
EXIT:
end SHIFT;
comment ================================================================;
procedure ACCUM(BStart,BSize,TStart,TSize);
value BStart,BSize,TStart,TSize;
integer BStart,BSize,TStart,TSize;
comment *** Add contents of extended 'integer' at BStart size BSize ***
*** into extended 'integer' at TStart, size TSize. ***
****************************************************************;
begin integer Maxpos,One,Count,Carry;
switch S ≔ OOP,OFLO,SKIP,EXIT;
Maxpos ≔ 74877906943;
One ≔ ;
if BSize NOTEQ TSize then begin
PRINT(“\nIllegal parameters to procedure ACCUM”);
goto EXIT;
end;
TStart ≔ Start+BSize;
BStart ≔ Start+BSize;
Elliott(3,0,BSize ,0,2,0,Count );
Elliott(2,6,Carry ,0,2,6,4 );
LOOP: Elliott(0,0,TStart,1,3,0,0 );
Elliott(0,0,BStart,1,0,4,0 );
Elliott(4,3,OFLO ,0,4,0,SKIP );
OFLO: Elliott(0,3,MaxPos,0,2,2,4 );
SKIP: Elliott(0,4,Carry ,0,0,0,0 );
Elliott(0,0,TStart,1,2,0,0 );
Elliott(3,6,4 ,0,2,0,Carry );
Elliott(3,0,One ,0,2,7,BStart);
Elliott(0,0,0 ,0,2,7,TStart);
Elliott(3,7,Count ,0,0,1,0 );
Elliott(4,1,LOOP ,0,0,0,0 );
EXIT:
end CBTOD;
comment ================================================================;
procedure CBTOD(BStart,BSize,DStart,DSize);
value BStart,BSize,DStart,DSize;
integer BStart,BSize,DStart,DSize;
comment *** Converts binary fraction in 'array' at BStart length ***
*** 1+BSize 'to' Decimal fraction at DStart length 1+DSize ***
****************************************************************;
begin integer Maxpos,Power,One,BCount,DCount,Carry,Item,Point;
switch S ≔ TART,LOOP,OFLO,SKIP,EXIT;
Maxpos ≔ 74877906943;
Power ≔ 0000000000;
One ≔ ;
Elliott(3,0,DSize ,0,2,1,DCount);
START: Elliott(0,6,BStart,1,1,0,0 );
Elliott(0,0,DStart,1,2,0,0 );
Elliott(3,0,DCount,0,4,2,EXIT );
Elliott(2,2,DCount,0,2,6,Carry );
Elliott(3,0,BSize ,0,2,1,BCount);
Elliott(3,0,BStart,0,0,4,BSize );
Elliott(2,0,Point ,0,0,0,0 );
LOOP: Elliott(0,0,Point ,1,3,0,0 );
Elliott(5,2,power ,0,2,0,4 );
Elliott(5,7,0 ,0,0,4,Carry );
Elliott(4,3,OFLO ,0,4,0,SKIP );
OFLO: Elliott(0,3,Maxpos,0,2,2,4 );
SKIP: Elliott(0,0,Point ,1,2,0,0 );
Elliott(3,0,4 ,0,2,0,Carry );
Elliott(3,0,One ,0,2,7,Point );
Elliott(3,2,BCOUNT,0,4,1,LOOP );
Elliott(2,2,DStart,0,4,0,START );
EXIT:
end CBTOD;
comment ================================================================;
procedure PDARY(DStart,DSize);
value DStart,DSize;
integer DStart,DSize;
comment *** PRINT(decimal number at DStart, and fractional component ***
*** at DStart+1, length DSize ***
****************************************************************);
begin integer Cols,Count,Word;
Cols ≔ ;
Elliott(0,0,Dstart,1,3,0,0 );
Elliott(2,0,Word ,0,0,0,0 );
PRINT(“\n”,sameline,Digits(1),Word,“.”);
for Count ≔ step 1 until DSize do begin
Elliott(2,2,Dstart,1,3,0,0 );
Elliott(2,0,Word ,0,0,0,0 );
if ( Count-1 ) = ( ( Count-1 ) ÷ Cols )×Cols then begin
if ( Count-1 ) ⩾ 1 then PRINT(“\n ”); end;
PRINT(sameline,leadzero(“0”),Digits(10),Word);
end;
end PDARY;
comment ================================================================;
integer procedure ESTIM(Expone);
value Expone;
integer Expone;
comment *** Estimate number of decimal places precision available ***
*** after term Expone contributed. ***
****************************************************************;
begin real Temp;
Temp ≔ /(8×Expone+1)-2/(8×Expone+4)-1/(8×Expone+5)-1/(8×Expone+6);
ESTIM ≔ ntier((Expone×LN(16)-LN(Temp))/LN(10)+0·5);
end ESTIM;
comment ================================================================;
procedure TRACE(TEXT);
string TEXT;
comment *** Trace program flow if bit 1 set ***
****************************************************************;
begin integer One;
switch S ≔ OTRAC;
One ≔ ;
Elliott(7,0,0 ,0,0,3,One );
Elliott(4,2,NOTRAC,0,0,0,0 );
PRINT(TEXT);
NOTRAC:
end TRACE;
comment ================================================================;
procedure CALPI(MaxTrm,BMax,DMax);
value MaxTrm,BMax,DMax;
integer MaxTrm,BMax,DMax;
comment *** Perform calculation of PI using algorithm described above***
*** MaxTrm is the number of terms to be calculated ***
*** BMax is words allocated to fractional part of result ***
*** DMax is words allocated (10 digits/word) to decimal ***
*** part of result ***
****************************************************************;
begin integer array BFrac1[0:BMax],BFrac2[0:BMax],BFrac3[0:BMax],BFrac4[0:BMax];
integer array BFracT[0:BMax],BShift[0:BMax],BReslt[0:BMax],BOutpt[0:BMax];
integer array DOutpt[0:DMax];
integer Two,Count,Term,Expone,EPlace,DWords;
switch S ≔ HOW,NOSHOW;
Two ≔ ;
for Count ≔ step 1 until BMax do
BReslt[Count] ≔ ;
for Term ≔ step 1 until MaxTrm do begin
Expone ≔ erm-1;
PRINT(“\nAdding term ”,sameline,digits(4),Term);
INTDV(4,8×Expone+1,ADDRESS(BFrac1),BMax);
comment PWRDS(`BFrac1',ADDRESS(BFrac1),BMax);
INTDV(2,8×Expone+4,ADDRESS(BFrac2),BMax);
comment PWRDS(`BFrac2',ADDRESS(BFrac2),BMax);
INTDV(1,8×Expone+5,ADDRESS(BFrac3),BMax);
comment PWRDS(`BFrac3',ADDRESS(BFrac3),BMax);
INTDV(1,8×Expone+6,ADDRESS(BFrac4),BMax);
comment PWRDS(`BFrac4',ADDRESS(BFrac4),BMax);
TRACE(“\nCombine fractions...”);
COMBI(ADDRESS(BFrac1),ADDRESS(BFrac2),ADDRESS(BFrac3),ADDRESS(BFrac4),ADDRESS(BFracT),BMax);
comment PWRDS(`BFracT',ADDRESS(BFracT),BMax);
comment TRACE(`\nShift down...');
SHIFT(Expone,ADDRESS(BFracT),ADDRESS(BShift),BMax);
PWRDS(“BShift”,ADDRESS(BShift),BMax);
TRACE(“\nAccumulate...”);
ACCUM(ADDRESS(BShift),BMax,ADDRESS(BReslt),BMax);
comment PWRDS(`BReslt',ADDRESS(BReslt),BMax);
if Term ⩾ MaxTrm-2 then goto SHOW;
Elliott(7,0,0 ,0,0,3,Two );
Elliott(4,2,NOSHOW,0,0,0,0 );
SHOW: TRACE(“\nConvert to Decimal...”);
for Count ≔ step 1 until BMax do BOutpt[Count] ≔ Reslt[Count];
CBTOD(ADDRESS(BOutpt),BMax,ADDRESS(DOutpt),DMax);
EPlace ≔ STIM(Expone);
DWords ≔ EPlace+15) ÷ 10;
if DWords GR DMax then DWords ≔ Max;
PRINT(SAMELINE,“ yields value for PI:-”);
PDARY(ADDRESS(DOutpt),DWords);
PRINT(“\nEstimated accurate to”,sameline,digits(4),EPlace,“ digits\n”);
NOSHOW: end;
end CALPI;
comment ================================================================
================================================================
================================================================;
begin integer N,BMax,DMax,Expone,MaxTrm;
real DTB;
switch S ≔ OOP1,LOOP2,ONYERBIKE;
Reader(2);
Punch(3);
Digits(4);
comment The most efficient use of array space is when both Binary and Decimal arrays
hold the same precision, which occurs when DMax = BMax * ( LN(2**39)/LN(10**10))
= BMax * 1.17402
To achieve an accuracy of D decimal digits, this is guaranteed by
using number of iterations given by: Maxtrm = D * ( LN(10)/LN(16) )
= D * 0.830482
where DMax decimal array entries hold 10*DMax digits
For a more accurate estimate, keep increasing Expone until
( Expone*LN(16)-LN(4/(8*Expone+1)-2/(8*Expone+4)-1/(8*Expone+5)-1/(8*Expone+6)) )/LN(10) >D
and then set MaxTrm :=Expone+1;
Read N;
if N ⩽ 0 then begin
PRINT(“Invalid number”);
goto ONYERBIKE;
end;
DMax ≔ N+9 ) ÷ 10;
BMax ≔ ntier( DMax×(10×LN(10))/(39×LN(2)) )+2;
PRINT(“\nRequired decimal places ”,sameline,N);
PRINT(“\nArray size, Decimal storage”,sameline,DMax);
PRINT(“\nArray size, Binary storage ”,sameline,BMax);
PRINT(“\nOverall array space reqd. ”,sameline,DMax+8×BMax);
Expone ≔ ;
LOOP1: Expone ≔ xpone+10;
if ESTIM(Expone) ⩽ (N-1) then goto LOOP1;
LOOP2: Expone ≔ xpone-1;
if ESTIM(Expone) ⩾ (N+1) then goto LOOP2;
Expone ≔ xpone+1;
Maxtrm ≔ xpone+1;
PRINT(“\nRequired number iterations ”,sameline,Maxtrm);
PRINT(“\fCalculation of PI using”,sameline,MaxTrm,“ iterations”);
PRINT(“\n=======================================”,“\n”);
CALPI(MaxTrm,BMax,DMax);
ONYERBIKE:
end;
comment ================================================================;
end;