The Museum of HP Calculators

HP Forum Archive 14

[ Return to Index | Top of Index ]

CALCULATING MANY DIGITS OF PI
Message #1 Posted by Gene Wright on 19 Mar 2004, 9:40 a.m.

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.

      
Re: CALCULATING MANY DIGITS OF PI
Message #2 Posted by HrastProgrammer on 19 Mar 2004, 9:54 a.m.,
in response to message #1 by Gene Wright

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.

            
Email me ...
Message #3 Posted by Gene on 19 Mar 2004, 11:22 a.m.,
in response to message #2 by HrastProgrammer

and I'll forward your email to Palmer.

Thanks! Gene

            
Re: CALCULATING MANY DIGITS OF PI
Message #4 Posted by Palmer O. Hanson, Jr. on 20 Mar 2004, 8:00 p.m.,
in response to message #2 by HrastProgrammer

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.

      
Re: CALCULATING MANY DIGITS OF PI
Message #5 Posted by David Smith on 19 Mar 2004, 12:25 p.m.,
in response to message #1 by Gene Wright

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).

            
Re: CALCULATING MANY DIGITS OF PI
Message #6 Posted by Nelson M. Sicuro (Brazil) on 19 Mar 2004, 1:39 p.m.,
in response to message #5 by David Smith

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

                  
Re: CALCULATING MANY DIGITS OF PI
Message #7 Posted by Gordon Dyer on 19 Mar 2004, 3:14 p.m.,
in response to message #6 by Nelson M. Sicuro (Brazil)

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

                        
Routine to calculate PI using only integer math [LONG, Pascal code]
Message #8 Posted by Nelson M. Sicuro (Brazil) on 19 Mar 2004, 5:43 p.m.,
in response to message #7 by Gordon Dyer

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

            
Approximating PI on a IBM 1130
Message #9 Posted by Andrés C. Rodríguez (Argentina) on 19 Mar 2004, 4:38 p.m.,
in response to message #5 by David Smith

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)

                  
Re: Approximating PI on a IBM 1130
Message #10 Posted by David Smith on 19 Mar 2004, 6:18 p.m.,
in response to message #9 by Andrés C. Rodríguez (Argentina)

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

                        
Re: Approximating PI on a IBM 1130
Message #11 Posted by Andrés C. Rodríguez (Argentina) on 19 Mar 2004, 7:15 p.m.,
in response to message #10 by David Smith

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)

                              
Re: Approximating PI on a IBM 1130
Message #12 Posted by David Smith on 20 Mar 2004, 12:41 p.m.,
in response to message #11 by Andrés C. Rodríguez (Argentina)

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.

                                    
Re: Approximating PI on a IBM 1130
Message #13 Posted by Andrés C. Rodríguez (Argentina) on 20 Mar 2004, 5:02 p.m.,
in response to message #12 by David Smith

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

      
Re: CALCULATING MANY DIGITS OF PI
Message #14 Posted by Gordon Dyer on 19 Mar 2004, 3:30 p.m.,
in response to message #1 by Gene Wright

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!!)

            
Re: CALCULATING MANY DIGITS OF PI - HP-42S
Message #15 Posted by Gordon Dyer on 19 Mar 2004, 4:27 p.m.,
in response to message #14 by Gordon Dyer

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.

      
PI approximation by the Monte Carlo method (long, examples)
Message #16 Posted by Andrés C. Rodríguez (Argentina) on 19 Mar 2004, 4:30 p.m.,
in response to message #1 by Gene Wright

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.

      
Re: CALCULATING MANY DIGITS OF PI
Message #17 Posted by Olivier De Smet on 20 Mar 2004, 7:02 a.m.,
in response to message #1 by Gene Wright

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

            
Re: CALCULATING MANY DIGITS OF PI
Message #18 Posted by Katie on 20 Mar 2004, 11:10 a.m.,
in response to message #17 by Olivier De Smet

Olivier,

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

-Katie

                  
Re: CALCULATING MANY DIGITS OF PI
Message #19 Posted by Olivier De Smet on 20 Mar 2004, 12:02 p.m.,
in response to message #18 by Katie

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

                        
Re: CALCULATING MANY DIGITS OF PI
Message #20 Posted by Katie on 20 Mar 2004, 9:37 p.m.,
in response to message #19 by Olivier De Smet

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

            
Re: CALCULATING MANY DIGITS OF PI
Message #21 Posted by Eddie Shore on 20 Mar 2004, 2:29 p.m.,
in response to message #17 by Olivier De Smet

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

            
Hmm ;)
Message #22 Posted by Tizedes Csaba [Hungary] on 21 Mar 2004, 2:21 p.m.,
in response to message #17 by Olivier De Smet

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

                  
Little challenge - I'll provocative
Message #23 Posted by Tizedes Csaba [Hungary] on 23 Mar 2004, 3:00 a.m.,
in response to message #22 by Tizedes Csaba [Hungary]

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

                        
Re: Little challenge - I'll provocative
Message #24 Posted by Olivier De Smet on 23 Mar 2004, 4:17 a.m.,
in response to message #23 by Tizedes Csaba [Hungary]

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.

                              
Re: Little challenge HP-71B
Message #25 Posted by Gordon Dyer on 23 Mar 2004, 8:05 p.m.,
in response to message #24 by Olivier De Smet

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

                                    
Re: Little challenge HP-71B
Message #26 Posted by Karl Schneider on 24 Mar 2004, 12:00 a.m.,
in response to message #25 by Gordon Dyer

Gordon --

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

Thanks,

-- Karl S.

                                          
Re: Little challenge HP-71B
Message #27 Posted by Gordon Dyer on 26 Mar 2004, 4:46 p.m.,
in response to message #26 by Karl Schneider

No Math ROM required...

                                    
Re: Little challenge HP-71B
Message #28 Posted by Katie on 24 Mar 2004, 12:53 p.m.,
in response to message #25 by Gordon Dyer

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.

                                          
Re: Little challenge HP-71B
Message #29 Posted by Palmer O. Hanson, Jr. on 26 Mar 2004, 7:53 p.m.,
in response to message #28 by Katie

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.

      
Re: CALCULATING MANY DIGITS OF PI
Message #30 Posted by Katie on 26 Mar 2004, 1:35 a.m.,
in response to message #1 by Gene Wright

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

            
Re: CALCULATING MANY DIGITS OF PI
Message #31 Posted by David Smith on 26 Mar 2004, 12:59 p.m.,
in response to message #30 by Katie

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

                  
Re: Favorite Programming Commands
Message #32 Posted by Paul Brogger on 26 Mar 2004, 1:05 p.m.,
in response to message #31 by David Smith

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

                        
Re: Favorite Programming Commands
Message #33 Posted by Ed Look on 26 Mar 2004, 10:30 p.m.,
in response to message #32 by Paul Brogger

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!


[ Return to Index | Top of Index ]

Go back to the main exhibit hall