CALCULATING MANY DIGITS OF PI



#29

Gene: This is a write-up I just received from Palmer Hanson, who edited the TI PPC Notes for many years in the 1980s - early 1990s. Thought it might be interesting. Looks like the race to 1000 digits on a machine prior to 1990 now belongs to TI? Can we beat the 5.81 hours on the TI95 with the HP71B? :-)
-----------------------------------------------

CALCULATING MANY DIGITS OF PI
by Palmer Hanson, editor of TI PPC Notes

Students of calculus may remember that converging series can be used for calculation of trigonometric and exponential functions; e.g.,

sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...

pi = 4(1 - 1/3 + 1/5 - 1/7 + 1/9 - ...

where the series for sin(x) converges rapidly but the series for pi converges very, very slowly. That equation for finding pi is ascribed to Leibniz and is a special case of the series for finding the arc tangent

arctan(x) = x - x^3/3 + x^5/5 - x^7/7 ...

evaluated for x = 1/4 radian. For more rapid convergence pi may be calculated as the sum of two arc tangents

pi = 16*arctan(1/5) - 4*arctan(1/239)

In June 1980 the PPC Calculator Journal (V7N5P9-11) presented an HP-41 program by Ron Knapp which would calculate 1000 digits in 11.5 hours and 1160 digits of pi in 15.25 hours. Editor Richard Nelson wrote to Maurice Swinnen, the editor of TI PPC Notes, saying "... if the TI-59 club wants to accept a real challenge, compute the first 1,000 digits of pi in less than 11 hours, 30 minutes. Printout time of the answer need not be included." The challenge was repeated in the "Calcu-Letter" column of the July 1981 issue of Popular Science. A revised version of Ron's program which would calculate 1000 digits of pi in 8.5 hours was published in late 1981 (V8N6P69 of the PPC Calculator Journal). There was no response in the pages of TI PPC Notes until mid 1982 when Bob Fruit submitted a program which would find 460 digits in 6 hours 18 minutes (V7N4/5P26). The size of the TI-59 memory prevented calculation of additional digits using Bob's method. Bob's program runs in three stages. The first stage enters fast mode and calculates the value of 16*arctan (1/5). The second stage calculates the value of 4*arctan(1/239) and subtracts it from the value found in the first stage. The third stage formats and prints the results. It seemed likely that the TI-59 programmers would not be able to overcome the advantage of the larger memory and faster speed of the HP-41. Bob's submission also included a table of the first 10,000 digits of pi (V7N4/5P28).

In February 1983 TI PPC Notes (V8N1P21) reported that Jovan Puzovic of Yugoslavia had found a way to calculate 1188 digits of pi with a TI-59, but even with his program running in fast mode the solution required eighty hours. The program uses the same equations as that of Bob Fruit and proceeds in three stages. First the arctan(1/5) portion is calculated and stored on magnetic cards. This requires about 62 hours. Second, the arctan(1/239) is calculated and stored on magnetic cards. That requires an additional 18 hours. Finally, a short routine combines the two intermediate solutions. It seemed likely that 1188 digits was about as much as could be obtained considering the memory limitations of the TI-59.

In April 1983 TI PPC Notes (V8N2P4) reported that the French scientific journal Science et Vie had published a program by Renaud de La Taille which could deliver 1287 digits of pi on a TI-59. The solution occupies 13 digits per data register times 99 registers for the total of 1287 digits with only one data register (R00) reserved for dsz control. The program is painfully S - L - O - W ; it runs for 24.55 days! In June 1983 TI PPC Notes (V8N3P8) published the program and instructions for its use. In August 1983 TI PPC Notes (V8N4P21) published a Maurice Swinnen's translation of the Science et Vie article which explained the method of calculation as follows:

"The entire art consists of finding the right formula in order to design a short and simple program, and to use the largest possible number of memory registers. ... ... start with the series

pi/2 = 1 + 1/3 + (1*2)/(3*5) + (1*2*3)/(3*5*7) + (1*2*3*4)/(3*5*7*9) + ...

Thus, pi/2 = ((...((2n/(2n+1) + 2)(n-1)/(2n-1) + 2)(n-2)/(2n-3) + 2)(n-3)/...)1/3 +2

This method avoids the addition of terms in the series to the preceding sums. The resulting program contains only a multiplication and a division. ..."

The same issue of TI PPC Notes presented Palmer Hanson's modification of the Science et Vie program to run in fast mode and complete the calculations in 13.39 days. Palmer's modification also provided automatic calculation of the number of registers and iterations to be used when calculating fewer than 1287 digits. Correspondents of Science et Vie used similar techniques with other programmable calculators and reported the following results

250 digits of pi on an HP-67

576 digits of pi on a TI-58

and, finally, 3600 digits of pi on an HP-41, a calculation which requires four months! The longest continuous calculation for a TI-59 may have been a search for "amicable numbers" by Bob Fruit. V8N5P13 of TI PPC Notes reported that he ran his search program for 2662 hours -- nearly 111 days. If nothing else the successful completion of calculations which extend to more than one hundred days attests to the reliability of the devices used.

In 1986 Robert Prins of the Netherlands published a book "Extended Precision Programs for the TI-59". He included a modification of Jovan Puzovic's program which allowed the user to calculate fewer digits than 1188 and a minor modification of the fast mode version of the Science et Vie program. He also included a TI-66 program which could calculate 494 digits of pi in about 7.5 days

In mid 1987 Hewlett Ladd reported that he had translated the Science et Vie TI-59 program for use with the TI-95. His translation would calculate up to 1573 digits. He did not report run times at that time. His program was resurrected for this historical account. It will deliver 1000 digits of pi in 40.95 hours. .

An HP-32SII program by Katie Wasserman which will calculate 99 digits of pi in 11 minutes appears at hpmuseum.org as Articles Forum No. 281 submitted in June 2002. That program uses a different formula for pi

pi = 20*arctan(1/7) + 8*arctan(3/79)

One might think that a savings in execution time would result because the first term would converge more rapidly; however, the second term will converge more slowly such that the total execution time will be expected to be very similar to that with the equation used by Knapp, Fruit and Puzovic.

In early 2004 Palmer Hanson noted that the five-fold increase in speed of Hewlett Ladd's translation of the TI-59 Science et Vie program for the TI-95 suggested that a conversion of Bob Fruit's TI-59 program for use on the TI-95 might yield results which were faster than the Knapp programs on the HP-41. That turned out to be true. The resulting program will deliver 30 places in 30 seconds, 90 places in 3 minutes 5 seconds, 99 places in 3 minutes 43 seconds, 200 places in 13 minutes 35 seconds, 460 places in 1 hour 8 minutes, 1000 places in 5.81 hours, 1160 places in 7.72 hours and 2070 places in 24.61 hours.

Epilog:

Other discussions of methods for finding or remembering pi were reported in later issues of TI PPC Notes:

V8N5P10 reported that George Thomson had recognized the editor's fascination with pi by addressing letters to "Palmer Pi Hanson" and passing on a mnemonic for remembering the first thirty digits of pi which had appeared on page 309 of the April 1969 issue of Datamation, where one counts the letters in each word to obtain each successive digit of pi.

Now I, even I, would celebrate in rhymes inept,
The great immortal Syracusan, rivaled nevermore,
Who in his wondrous lore passed on before,
Left men his guidance how to circles mensurate.

V10N1P14/15 discussed a method for finding pi through the use of random numbers which appeared in A. K. Dewdney's column "Computer Recreations" in the April 1985 issue of Scientific American. The analogy is to a cannon firing into a square field which includes an inscribed circular pond. If the cannon is fired many times and the hits are distributed uniformly over the area of the square then the ratio of the number of hits in the pond to the total number of hits in the square should approach the ratio of the area of the circle to the area of the square which is pi/4. The simulation on a calculator or a computer is performed by using a random number generator to determine the position of each hit. Clearly, such a method for finding pi depends upon the quality of the random number generator. John Von Neumann's admonition that "Anyone who considers mathematical methods of producing random numbers is, of course, in a state of sin" (Donald Knuth's The Art of Computer Programming, Vol. 2, page 1) may be applicable here.

V10N3P4 noted that approximating pi with the fractions 22/7 or 355/113 was a well known technique. Ron Wagner had obtained another technique known as the Hick's (or maybe Hicks) method at a computer trade show. Divide 2143 by 22 and take the square root two times. On the TI-59 the fraction 22/7 yields three correct digits, the fraction 355/113 yields seven correct digits and the Hick's method yields nine correct digits where the tenth digit is within one of the tenth digit displayed for pi. On the HP-41 the Hick's method yields 3.141592653 which is within one in the tenth digit of the value of pi used in that device.

V10N4P26 reported that Larry Leeds had used a decimal to fraction converter on his Model 100 to search for fractions and fractions followed by square roots which would return the first thirteen digits of pi. He found three fractions: 4272943/1360120, 5419351/1725033 and 61905677/19705189,. He also found three fractions followed by a square root: 3044467/308469, 17007401/1723210 and 26140802/2648617. He did not find any fractions which when combined with two square roots would yield thirteen correct digits of pi. The first thirteen digits of pi 3.14159 26535 89 (obtained by truncation) are NOT the same as the rounded value 3.14159 26535 90 which is obtained in response to the pi key on the TI-59, TI-66 and TI-95. It is more than a little curious that the TI-59 and TI-66, which use truncated thirteen digit arithmetic, store an internal value for pi which is a thirteen digit rounded value. It seems likely that memorizing thirteen digits of pi would be easier than memorizing any of those methods.

The discussion in V10N4P26 also proposed a simple-minded search which would compare the decimal equivalent of fractions against the value of pi (or any other number to be converted) loaded to the accuracy of the machine. Start with the fraction 1/1. If the value of the fraction is greater than the decimal to be converted then add a 1 to the denominator and try again. If the value of the fraction is less than the number to be converted then add a 1 to the numerator and try again. Sample programs for a CC-40 were included. Use of the simple minded search on a TI-95 yields two fractions of interest. 5419351/1725033 = 3.14159 26535 8985 ... which becomes 1E-12 less than the defined value on machines that truncate to thirteen digits such as the TI-59 and TI-66, and which becomes the defined value on machines that round to thirteen digits such as the TI-95 and 6565759/2089946 = 3.14159 26565 9009 ... which becomes the defined thirteen digit value whether truncated or rounded. The same technique on an HP-41 yields fractions such as 104348/33215 and 209051/66543 which will yield the ten digit value used for pi on that device.


#30

Hi Gene,

Do you have this PI calculating program for TI-95? I would like to test it on my TI-95 emulator ...

Best regards.


#31

and I'll forward your email to Palmer.

Thanks!
Gene

#32

I have a pasteup of the PC-324 printout of the program. It includes about 750 steps so I don't want to go through the labor of typing them all in. I haven't figured out how to transmit material like the pasteup through the museum so if you will send an e-mail address I will send it as an attachment.

#33

In high school I wrote a program (FORTRAN) for calculating 1 million digits of pi on an IBM1130. It actually calculated arctan(1) and multiplied it by 4. It ran over Christmas break (around two weeks).


#34

I have a Pascal (Delphi) program that calculates PI with any decimal places (initially with 1000, easy to change) using ONLY INTEGER OPERATIONS!

If there is someone interested, I post here (code not so big).

Best regards,

Nelson


#35

I am interested in your Pascal code for PI. Please post it in here.

Thanks.


#36

There is the routine (is a pascal unit, need to be called from another [main] program).

unit expi;

interface

function CalcPI(numdig:integer): string;

procedure TestUnit;

implementation
{/* +++Date last modified: 05-Jul-1997 */
/*
This program was based on a Pascal program posted in the FIDO 80XXX
assembly conference. That Pascal program had the following comment:
------------------------
This program, which produces the first 1000 digits of PI, using
only integer arithmetic, comes from an article in the "American
Mathematical Monthly", Volume 102, Number 3, March 1995, page
195, by Stanley Rabinowitz and Stan Wagon.
------------------------
My C implementation is placed into the public domain, by its author
Carey Bloodworth, on August 22, 1996.

I have not seen the original article, only that Pascal implementation,
but based on a discussion in the 80XXX conference, I believe the
program should be accurate to at least 32 million digits when using
long unsigned integers. Using only 16 bit integers causes the
variables to overflow after a few hundred digits.
*/
}

var
Digits:string;
PDig:pchar;
pos:integer;

procedure SetupDigits(len:integer);
begin
SetLength(Digits,len+2);
Digits[1]:='3';
Digits[2]:='.';
pos:=3;
Digits[pos]:=#0;
PDig:=pchar(Digits);
end;

procedure OutDig(dig:integer);
begin
Digits[pos]:=char(ord('0')+dig);
inc(pos);
end;


type
PCardinalArray=^TCardinalArray;
TCardinalArray=array[0..0] of cardinal;

function Inner(pi:PCardinalArray; len:cardinal):cardinal;
var
i,x,k:cardinal;
j:cardinal;
begin
j:=10;
result:=0;
i:=len;
while i>0 do // converted for to while
begin
x:=j*pi[i]+result*i; // converted 10 to a variable 'j' (PII optimization only)
k:=i+i-1;
result:=x div k;
pi[i]:=x-result*k; // eliminated mod (which is a division)
dec(i);
end;
end;

function CalcPI(numdig:integer): string;
var
i,j,nines,predigit:cardinal;
len,q,qi:cardinal;
pi:array of cardinal;
begin
len := (numdig*10) div 3;
SetLength(pi,len+1);
for i:=1 to len do
pi[i]:=2;
nines := 0;
predigit := 0;
SetupDigits(len);
for j:=0 to numdig do
begin
qi:=Inner(pointer(pi),len);
q:=qi div 10;
pi[1]:=qi-q*10;
if (q=9) then
inc(nines)
else if q=10 then
begin
OutDig(1+predigit);
while nines>0 do
begin
OutDig(0);
dec(nines);
end;
predigit:=0;
end
else
begin
if j>1 then
OutDig(predigit);
while nines>0 do
begin
OutDig(9);
dec(nines);
end;
predigit:=q;
end;
end;
OutDig(predigit);
Result := Digits;
end;

procedure TestUnit;
begin
CalcPI(1000);
end;

end.

Best regards,

Nelson

#37

Coincidence!. In 1978, in my first programming course (EE 3rd year), I wrote a FORTRAN program to calculate PI approximations, running on a IBM 1130.

Oh!, I had to use punched cards for the 1130, but, at the same time, I was doing similar things on the HP-25 (see long post below)


#38

//JOB
//FOR
*IOCS(1132PRINTER,CARD,KEYBOARD,PLOTTER,DISK)
//


#39

Exactly the same, but as far as I recall, our IOCS had no plotter... And we were not allowed to program anything using the console (100% batch, cards with data followed the program cards). All output was hardcopy.

Main memory was 16 K, and the large rotating disks (looked like car wheels) held a mere hundred K... (less than 1 MBy, I mean)


#40

Actually the disks bigger than 100K... at least big enough to store 1 million digits. Our machine only had 8K of memory. It could still run Fortran though. A 24 pass compiler. It could even run ECAP.


#41

You are right, I think they held "a mere hundreds of K", but I made a mistake omitting the "s" without noticing the error. Thank you

#42

A much faster way is by Gauss-Legendre.



Here is the algorithm:

define numbers a,b,c such that at the initial state:

a=x=1

b=1/sqrt(2)

c=1/4



then iterate:

y=a

a=(a+b)/2

b=sqrt(b*y)

c=c-x(a-y)*(a-y)

x=2x



such that pi=(a+b)*(a+b)/4c



This algorithm doubles the numbers of pi every iteration!!

The algorithm has second order convergent nature. Then if you want to calculate up to n digits, aniteration count of the order log2(n) is sufficient.

E.g. 19 times for 1 million decimal digits, 31 times for 3.2 billon decimal digits.

(This algorithm is used for the world-record on calculating pi!!)


#43

Here is an HP-42S program for the Gauss-Legendre method.

Registers used are:

a= reg 01, b= reg 02, c=reg 03, x=reg 04, y= reg 05



PI Program for HP-42S

uses Gauss-Legendre

Converges to 11 decimal places in 3 iterations!



01 LBL "PI-GL"

02 1

03 STO 01

04 2

05 SQRT

06 1/X

07 STO 02

08 4

09 1/X

10 STO 03

11 1

12 STO 04

13 LBL 01

14 RCL 01

15 STO 05

16 RCL 02

17 +

18 2

19 /

20 STO 01

21 RCL 02

22 RCL 05

23 x

24 SQRT

25 STO 02

26 RCL 01

27 RCL 05

28 -

29 X^2

30 RCL 04

31 x

32 +/- (CHS)

33 RCL 03

34 +

35 STO 03

36 RCL 04

37 2

38 x

39 STO 04

40 RCL 01

41 RCL 02

42 +

43 X^2

44 RCL 03

45 4

46 x

47 /

48 PSE

49 GTO 01

50 .END.

#44

Quote:

V10N1P14/15 discussed a method for finding pi through the use of random numbers which appeared in A. K. Dewdney's column "Computer Recreations" in the April 1985 issue of Scientific American. The analogy is to a cannon firing into a square field which includes an inscribed circular pond. If the cannon is fired many times and the hits are distributed uniformly over the area of the square then the ratio of the number of hits in the pond to the total number of hits in the square should approach the ratio of the area of the circle to the area of the square which is pi/4. The simulation on a calculator or a computer is performed by using a random number generator to determine the position of each hit. Clearly, such a method for finding pi depends upon the quality of the random number generator. John Von Neumann's admonition that "Anyone who considers mathematical methods of producing random numbers is, of course, in a state of sin" (Donald Knuth's The Art of Computer Programming, Vol. 2, page 1) may be applicable here.

Such problem and a proposed solution appeared in the documentation of the TI-56 (IIRC, a keystroke programmable, contemporary of the TI-52), circa 1977. The TI-56 was (more or less) in the HP-25 class; and the TI-52, with its magnetic cards, was TI's not-so-good answer to the HP-65.

The problem was used to illustrate the MonteCarlo method, and the usage of random numbers.

Below you will find a couple of HP-25 ports of the solution (just from my old memories, not for actual usage). This congruential random number generator is supposed to pass the spectral test as per Knuth Vol. II, variations of it appeared in HP-25 and HP-41 manuals.

The first version:

01 RCL 0; random number seed

02 RCL 1; contains 9821

03 x

04 RCL 2; contains 0.211327

05 +

06 FRAC; obtain first random number, also used as seed

07 STO 0; saved as new seed

08 RCL 0; the first random number stays in Y register

09 RCL 1

10 x

11 RCL 2

12 +

13 FRAC; obtain second random number

14 STO 0; saved as new seed

15 R=>P; converts the random numbers in X and Y in a "vector"

16 INT; a cheap manner to normalize results and separate "0" from "1" i.e.: "in" or "out" of the pond

17 E+; sums the "out" shots in R7, increments R3 as iteration counter

18 GTO 01; loop back, if desired RCL 3 and PAUSE can be inserted here to show progress

After running the program for some time, the MEAN function will give an approximation to (1-(PI/4)). The summation registers should be cleared before each run.

The second (not very different) version

01 RCL 0; random number seed

02 RCL 1; contains 9821

03 x

04 RCL 2; contains 0.211327

05 +

06 FRAC; obtain first random number, also used as seed

07 ENTER; is not really necessary to save the new seed,

08 ENTER; and the first random number is also kept in the Y register

09 RCL 1

10 x

11 RCL 2

12 +

13 FRAC; obtain second random number

14 STO 0; saved as new seed

15 R=>P; converts the random numbers in X and Y in a "vector"

16 1

17 STO + 3; increment iteration counter

18 X>Y?; in this case...

19 STO + 7; ... we count the "in" shots!

20 GTO 01; loop back, if desired RCL 3 and PAUSE can be inserted here to show progress

In this case, after running the program for some time, the MEAN function will give an approximation to (PI/4). The R3 and R7 registers should be cleared before each run.

The small differences between the programs appear in steps 07-08 and from step 16 on.

Edited: 20 Mar 2004, 10:24 a.m.

#45

Just for fun an adapted program from a 32Sii one for HP49G+
(using infinite integer):

%%HP: T(3)A(R)F(.);
\<< -105. CF \-> N
\<< 280000 N ALOG * DUP 'R' STO 2 \-> y
\<<
DO y * y 1 + IQUOT 50 IQUOT DUP 'R' STO+ 2 'y' STO+
UNTIL DUP 0 ==
END DROP
\>> 30336 N ALOG * DUP 'R' STO+ 2 \-> y
\<<
DO 9 * y * 6250 IQUOT y 1 + IQUOT DUP 'R' STO+ 2 'y' STO+
UNTIL DUP 0 ==
END DROP
\>>
\>> R
\>>

gives 1000 decimals of PI in 231 sec.

3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198937624


#46

Olivier,

Nice work! RPL and infinite integer arithmetic sure makes that algorithm a lot more understandable.

-Katie


#47

Katie,

I just remembered where I found the original 32SII program,
in the Articles Forum (99 Digits of PI on an HP 32SII) and just saw that you were the author of it ... :)

Thanks again for this nice algorithm

Olivier


#48

Olivier,

I'd love to take credit for the formula, but it's one that Euler came up with 250 years ago using the idea that Machin suggested 50 years earlier (see "A History of Pi"). I just did the tweaking to get it into an easy to express algorithm.

-Katie

#49

Very impressive everybody. I would like to learn how to do that.

#50

Okay, this is 'VERY' beautiful, but I think the TI95s haven't got 'infinite integer'-long variables... ;)

So, what about if we want to use traditional variable formats...?

Csaba


#51

So, that is seems to me, the calculating 1000 digits of PI is too big bite anybody here, with no-infinite-long variables...! ;)

Let's make a challenge, and try we to push down the TI's record!

The youngest calc is usable let the HP48 series!

Good working!

Csaba


#52

The same prog with no infinite-integer :

N 'PPI' warning N is the number of digits divided by 10 (i.e. 100 for ~1000 digits)

Only valid for hp49g(+) and not very optimized
(NIP, IREMAINDER, ... )

I'll try an hp48gx one later...

Olivier

Ln
« DUP 1000000000. MOD SWAP 1000000000. / IP OBJ-> 0. SWAP ->LIST NIP ADD »

Ld
« OVER SIZE UNROT
-> l d
« 0. SWAP 1. SWAP
FOR j 'l' j DUP2 GET 4. ROLL + d DUP2 6. ROLLD 6. ROLLD / IP PUT IREMAINDER 1000000000. *
NEXT DROP l
»
»

Lm
« * Ln
»

Lp
« ADD Ln
»

PPI
« -> N
« 280000000. 0. N 1. - NDUPN 1. + ->LIST 'R' STO R 2. -> y
«
DO y Lm y 1. + Ld 50. Ld DUP R Lp 'R' STO 2. 'y' STO+
UNTIL DUP \SigmaLIST DUP ->STR y ->STR " " + SWAP + 1. DISP 0. ==
END DROP
» 30336000. 0. N 1. - NDUPN 1. + ->LIST DUP R Lp 'R' STO 2. -> y
«
DO 9. Lm 125. Ld y Lm y 1. + Ld 50. Ld DUP R Lp 'R' STO 2. 'y' STO+
UNTIL DUP \SigmaLIST DUP ->STR y ->STR " " + SWAP + 1. DISP 0. ==
END DROP
»
» R
»


Edited: 23 Mar 2004, 7:14 a.m.


#53

Here is an 'infinite PI prog for the HP-71B.

I inherited this on some mag cards I bought but I haven't deciphered the algorithm yet.

The value of PI is stored into a text file called 'PI'



HP-71B Program 'PIVAL':

10 ON ERROR GOSUB 'ERR'

20 DELAY 0

30 DESTROY ALL

40 L=11 @ H=10^L

50 INPUT "No. of DP? ";D0

60 T1=TIME @ T9=T1

70 B=INT(D0/L)+2 @ T0=1.66*D0

80 DIM P(B),T(B)

90 T(B-1)=H/2 @ P(B-1)=H/2

100 FOR N=1 TO T0

110 T8=INT((TIME-T9)*(T0-N+1)+TIME)

120 IF N=1 THEN 170

130 D9=INT(T8/(24*3600))

140 T8=MOD(T8,24*3600)

155 IMAGE "TERM:"ZZZZ,XDD"D ",ZZ":",ZZ":",ZZ

160 DISP USING 155;T0-N,D9,INT(T8/3600),MOD(T8/60,60),MOD(T8,60)

170 T9=TIME @ X=N+N-1 @ GOSUB 380 @ GOSUB 380

180 X=8*N @ GOSUB 430

190 X=N+N+1 @ GOSUB 430

200 GOSUB 480 @ NEXT N

210 C=0 @ FOR i=1 TO B

220 P(I)=P(I)*6+C

230 C=INT(P(I)/H)

240 P(I)=P(I)-C*H

250 NEXT I

260 'SHOW':

270 PURGE PI @ CREATE TEXT PI @ ASSIGN #1 TO PI

280 ! PRINT #1;STR$(P(B))&"."

290 FOR I=B-1 TO 1 STEP -1

300 A$="000000000000"&STR$(P(I))

310 PRINT #1;A$[LEN(A$)+1-L]

320 NEXT I

330 DELAY .5

340 OFF TIMER #1

350 T1=TIME-T1 @ T0$=TIME$ @ OFF

360 DISP "CPU TIME=";T1-60*O0

370 OFF TIMER #1 @ END

380 C=0 @ FOR I=1 TO B

390 T(I)=T(I)*X+C

400 C=INT(T(I)/H)

410 T(I)=T(I)-C*H

420 NEXT I @ RETURN

430 C=0 @ FOR I=B TO 1 STEP -1

440 Z=T(I)+C @ C=0

450 Q=INT(Z/X) @ T(I)=Q

460 C=H*(Z-Q*X)

470 NEXT I @ RETURN

480 C=0 @ FOR I=1 TO B

490 P(I)=P(I)+T(I)+C @ C=0

500 IF P(I)<H THEN 520

510 P(I)=P(I)-H @ C=1

520 NEXT I @ RETURN

530 'ERR': T0$=TIME$

540 OFF

550 DISP "ERROR:";ERRM$;" in line";ERRL;"at ";T0$

560 PAUSE

570 RETURN



#54

Gordon --

Does this program require a Math ROM? I don't see any "high-level math" stuff in it...

Thanks,

-- Karl S.


#55

No Math ROM required...

#56

A little investigation reveals that this is using the series:

             1     1   9       1     9     25                  prod((2n-1)^2)   
pi/6= --- + --- ----- + --- ------ ------- + ..... + ------------------
2 2 8*1*3 2 8*1*3 8*2*5 prod(2*n*8*(2n+1))

which is :

pi/6 = arcsin(1/2)

Simply to program but not fast at converging.

Also, I think that this program will produce the wrong results after a short while because the 71b stores 12 BCD digits which will overflow when the program when n>=50. This can be easily fixed by lowering the value of L (in line 40) to 7 or 8, but you'll get fewer digits per array element.

Edited: 24 Mar 2004, 11:25 p.m.


#57

When I converted Bob Fruit's TI-59 program for the TI-95 I found that it did not obtain correct results when the number of decimal places(N)exceeded 770. Bob's program used ten digit registers. I fixed the problem by shifting to nine digit registers when N>770. This preserves the faster calculations when N<771 but some confusion would result when some outputs appeared in ten digit groups and others in nine digit groups. I added a routine to convert the nine digit groups to ten digit groups for output. This made it easier to check results since my table of 10,000 places of pi is in ten digit groups.

#58

Here's the same algorithm that I used on the 32SII but for the 71b. It's not quite as fast as the TI-95 implementation taking 7+ hours for 1000 digits. But it is reasonably memory efficient and will produce up to 100,000 digits of pi if you have enough memory, time and battery power! This could certainly be optimized for speed and probably beat the TI-95 record with little trouble, but I wrote this for clarity so that other's might do their own modifications:

1 ! Compute PI in integer array P() 5 digits per element
2 ! up to 100000 digits (or whatever memory allows)
3 !
4 ! Uses the formula: pi = 20*arctan(1/7) + 8*arctan(3/79)
5 !
6 ! and arctan(x) = (y/x) * (1 + 2/3*y + 2/3*4/5*y*y + 2/3*4/5*6/7*y*y*y + .....)
7 ! where y = x*x/(1+ x*x)
8 !
9 !
10 DESTROY ALL ! clear all variables
20 T1=TIME ! note the time
30 F=100000 ! 5 digits per array element
40 M=IP(MEM/3/2)-100 ! number of 3 byte array elements available
50 DISP "MAX DIGITS=";M*5; ! prompt the user for the number of digits
60 INPUT M1
70 N=M1 DIV 5
80 INTEGER P(N),S(N) ! need 2 arrays, P=pi, S=partial sum

90 FOR I=1 TO N @ S(I)=0 @ NEXT I @ S(1)=28000 ! first term starting value is .28000.....
100 GOSUB 'ADD' ! add this into P()
110 FOR T=2 TO INF STEP 2 ! loop until S() is zero
120 D=T @ GOSUB 'MULTIPLY' ! mutiply S() by T
130 D=50*(T+1) @ GOSUB 'DIVIDE' ! dividie S() by 50*(T+1)
140 IF Z=0 THEN 'TERM2' ! Check for end of loop
150 GOSUB 'ADD' ! add S() to P()
160 NEXT T

170 'TERM2': ! second term starting value is .030336000.....
180 FOR I=1 TO N @ S(I)=0 @ NEXT I @ S(1)=3033 @ S(2)=60000
190 GOSUB 'ADD' ! add this to P()
200 FOR T=2 TO INF STEP 2 ! loop unitl S() is zero
210 D=T*9 @ GOSUB 'MULTIPLY' ! multiply S() by T*9
220 D=T+1 @ GOSUB 'DIVIDE' ! divide S() by 6250*(T+1)
230 D=6250 @ GOSUB 'DIVIDE' ! (in 2 steps to avoid overflows)
240 IF Z=0 THEN 'FINISH' ! check for end of loop
250 GOSUB 'ADD' ! add S() to P()
260 NEXT T

270 'FINISH':
280 T2=TIME-T1 @ DISP T2;"SECONDS" ! show the time in seconds that it took
290 END

300 'ADD': ! add S() into P()
310 C=0 ! clear carry
320 FOR I=N TO 1 STEP -1
330 X=P(I)+S(I)+C ! do the add and include the carry
340 P(I)=MOD(X,F) @ C=X DIV F ! adjust for possible overflow and get new carry
350 NEXT I
360 RETURN

370 'DIVIDE': ! divide S() by D
380 C=0 @ Z=0 ! clear carry and zero flag
390 FOR I=1 TO N
400 X=S(I)+C*F ! add in carry shifted by 5 digits
410 S(I)=X DIV D @ C=MOD(X,D) ! compute quotient and remainder (carry)
420 IF S(I) THEN Z=1 ! set the zero flag for S()=0 check
430 NEXT I
440 RETURN

450 'MULTIPLY': ! multiply S() by D
460 C=0 ! clear carry
470 FOR I=N TO 1 STEP -1
480 X=S(I)*D+C ! do the multiply and add in the carry
490 S(I)=MOD(X,F) @ C=X DIV F ! compute the product and the new carry
500 NEXT I
510 RETURN


#59

"DESTROY ALL" ... definitely my favorite programming command.


#60

IIRC, the 6809 assembly language mnemonic for "Sign Extend" was "SEX" . . .


#61

Well, this one's not a programming command, but I saw it today, from a foreign scientific journal: the authors wanted to indicate a variable or unknown content of selenium atoms (Se) in a compound and somehow, in the abstract of the article, the "x" for variable number of atoms of Se didn't get printed as a subscript! And even if it did... !

I didn't share this one with my immediate company as it was also mixed company! And, I'm not making this up!


Possibly Related Threads...
Thread Author Replies Views Last Post
  [OT] Mathematica free for Raspberry PI BruceH 32 3,769 11-23-2013, 05:24 AM
Last Post: Nick_S
  Computing pi with the PC-1300S Kiyoshi Akima 0 499 11-17-2013, 12:24 AM
Last Post: Kiyoshi Akima
  [HP Prime] Calculating Prandtl Number with Units (bug found in USIMPLIFY) Timothy Roche 1 604 11-13-2013, 04:07 PM
Last Post: cyrille de Brébisson
  Calculating Pi LHH 9 1,323 09-27-2013, 10:50 PM
Last Post: Gerson W. Barbosa
  Visualization of pi Bruce Bergman 13 1,727 08-17-2013, 05:00 PM
Last Post: Howard Owen
  OT: Happy Pi Day! Eddie W. Shore 13 1,818 03-22-2013, 10:44 AM
Last Post: Les Koller
  Totally OT ... Pi Day for my car Maximilian Hohmann 18 2,293 03-10-2013, 01:15 PM
Last Post: chris smith
  [WP34S] A funny bug in Pi (prod) Eduardo Duenez 3 758 01-28-2013, 03:41 AM
Last Post: Walter B
  28S Pi Functionality Matt Fegenbush 3 749 10-17-2012, 02:15 AM
Last Post: Nick_S
  Obtaining More Decimal Digits (50g) Eddie W. Shore 3 753 09-13-2012, 05:38 PM
Last Post: Gilles Carpentier

Forum Jump: