The Museum of HP CalculatorsHP Articles Forum
[Return to the Index ]
[ Previous | Next ]
Many Digits of Pi (HP32sii, HP-16c, HP-12c, HP-12cp, HP-19c, HP-30b, HP-15c)
Posted by Katie Wasserman on 17 Oct 2008, 11:56 p.m.
Here are several fairly simple programs for calculating many digits of pi on various HP calculators. My goal in writing these is to see how many digits of pi can be calculated and saved in memory on calculators of limited memory (under 1K bytes).
XEQ P - results are in (A) - (Z) 8 digits per register no leading zeros - takes about 40 minutes.based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits
P01 | LBL P P02 | CLVARS P03 | CLSIGMA
P04 | 28 ; Set up to store the iteration counter in register n P05 | STO i
P06 | 26 ; calculate number of iterations needed = log-base-2(10) * N P07 | 8 P08 | x P09 | 2 P10 | LOG P11 | / P12 | IP P13 | STO(i) ; Store in n
P14 | 2 P15 | STO A ; starting value P16 | x P17 | 3 P18 | + P19 | SIGMA+ ; the SIGMAx register is used to store (iteration_counter x 2 + 1)
B01 | LBL B ; main loop B02 | 1.026 ; set up for register loop B03 | STO i B04 | 0 ; initial carry
C01 | LBL C C02 | 1E8 ; 100,000,000 this determines the number of digits per register C03 | x
C04 | n C05 | RCLx(i) C06 | + ; multiply and add carry
C07 | ENTER C08 | ENTER C09 | SIGMAx C10 | / C11 | IP ; divide by (iteration_counter x 2 + 1) C12 | STO(i)
C13 | SIGMAx C14 | x C15 | - ; remainder is the carry
C16 | ISG i C17 | GTO C
C18 | 2 C19 | STO+ A C20 | SIGMA- ; decrement iteration_counter by 1; and (2 x iteration_counter + 1) by 2 C21 | x>0? C22 | GTO B
C23 | 26 ; set up for register loop to carry overflows up one register, C24 | STO i ; register value can never be more than 2 x 10^8
E01 | LBL E E02 | RCL(i) E03 | 1E8 E04 | - E05 | x<=0? E06 | GTO F ; no carry E07 | STO(i) E08 | 1 E09 | STO- i ; these 3 instructions add the carry forward E10 | STO+(i) E11 | STO+ i E12 | GTO E ; this isn't really needed since the overflow can never be more than 1E8, ; but it's nice to have if other algorithms are used
F01 | LBL F F02 | DSE i F03 | GTO E F04 | RTN
Checksums:
P=28.5 CK=76B9 B=14.0 CK=334E C=44.0 CK=DFE2 E=26.0 CK=C6AE F=06.0 CK=BF13
GSB A - results are in (0) - (F) 13 digits per register no leading zeros - takes about 17 hours.
based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits
001 | LBL A | 43,22, A
002 | FLOAT 0 | 42,45, 0 ; Set 56 bit word size, 2's compliment 003 | DEC | 24 004 | CLR REG | 42 34 005 | 2s COMP | 42 2
006 | 6 | 6 ; 663 iterations = 195 digits, convergence is about 1 bit per iteration 007 | 6 | 6 008 | 3 | 3
009 | LBL B | 43,22, B 010 | STO .0 | 44 .0 ; register .0 is the iteration_counter
011 | 0 | 0 ; initial carry
012 | STO I | 44 32 ; register I is the register pointer (0 - F)
013 | LBL C | 43,22, C 014 | GSB 8 | 21 8 ; get 10,000,000,000,000 015 | X | 20 ; Carry x 10^13
016 | RCL .0 | 45 .0 017 | RCL (i) | 45 31 018 | X | 20 019 | + | 40 ; multiply and add carry
020 | Enter | 36 021 | Enter | 36
022 | RCL .0 | 45 .0 023 | SL | 42 A 024 | 1 | 1 025 | + | 40 026 | / | 10 ; divide by (iteration_counter x 2 + 1)
027 | STO (i) | 44 31 ; store result
028 | R down | 33 029 | LST X | 43 36 030 | RMD | 42 9 ; remainder is the carry
031 | X <> I | 42 22 ; increment register pointer and check for end of loop 032 | 1 | 1 033 | 5 | 5 034 | X = Y | 43 49 035 | GTO D | 22 D 036 | R down | 33 037 | 1 | 1 038 | + | 40 039 | X <> I | 42 22 040 | GTO C | 22 C
041 | LBL D | 43,22, D ; decrement iteration_counter and check for end 042 | 2 | 2 043 | STO 0 | 44 0 044 | RCL .0 | 45 .0 045 | 1 | 1 046 | - | 30 047 | X > 0 | 43 30 048 | GTO B | 22 B
049 | 1 | 1 ; loop through registers to carry any overflows up to 050 | 5 | 5 ; next register. Overflow can never be more than 2 x 10^13. 051 | STO I | 44 32 052 | 0 | 0 053 | LBL E | 43,22, E 054 | GSB 7 | 21 7 ; this is weird, but saves one byte 055 | GSB 8 | 21 8 056 | - | 30 057 | X < 0 | 43 2 058 | GTO F | 22 F 059 | STO (i) | 44 31 060 | 1 | 1 061 | GTO 9 | 22 9 062 | LBL F | 43,22, F 063 | 0 | 0 064 | LBL 9 | 43,22, 9 065 | DSZ | 43 23 066 | GTO E | 22 E 067 | LBL 7 | 43,22, 7 068 | RCL (i) | 45 31 069 | + | 40 070 | STO (i) | 44 31 071 | RTN | 43 21
072 | LBL 8 | 43,22, 8 ; put 10^13 on the stack 073 | 1 | 1 074 | 0 | 0 075 | 0 | 0 076 | 0 | 0 077 | Enter | 36 078 | X | 20 079 | Enter | 36 080 | X | 20 081 | 1 | 1 082 | 0 | 0 083 | X | 20 084 | RTN | 43 21
CLEAR REG R/S - results are in (0) - (.0) 7 digits per register no leading zeros - takes about 90 minutes.based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits
01 | EEX | 26 02 | 7 | 7 03 | STO FV | 44 15 04 | 2 | 2 05 | 4 | 4 06 | 0 | 0 07 | STO i | 44 12 08 | GTO 42 | 43,33 42 09 | 0 | 0 10 | STO n | 44 11 11 | STO PMT | 44 14 12 | RCL PMT | 45 14 13 | RCL FV | 45 15 14 | x | 20 15 | RCL CFj | 45,43 14 16 | RCL i | 45 12 17 | x | 20 18 | + | 40 19 | ENTER | 36 20 | ENTER | 36 21 | RCL i | 45 12 22 | 2 | 2 23 | x | 20 24 | 1 | 1 25 | + | 40 26 | STO PV | 44 13 27 | / | 10 28 | INTG | 43 25 29 | CFj | 43 14 30 | RCL PV | 45 13 31 | x | 20 32 | - | 30 33 | STO PMT | 44 14 34 | 1 | 1 35 | 0 | 0 36 | RCL n | 45 11 37 | 1 | 1 38 | + | 40 39 | STO n | 44 11 40 | x<=y | 43 34 41 | GTO 12 | 43,33 12 42 | 2 | 2 43 | STO 0 | 44 0 44 | RCL i | 45 12 45 | 1 | 1 46 | - | 30 47 | STO i | 44 12 48 | x=0 | 43 35 49 | GTO 51 | 43,33 51 50 | GTO 09 | 43,33 09 51 | 1 | 1 52 | 0 | 0 53 | STO n | 44 11 54 | RCL i | 45 12 55 | RCL CFj | 45,43 14 56 | + | 40 57 | RCL FV | 45 15 58 | / | 10 59 | INTG | 43 25 60 | STO i | 44 12 61 | LASTx | 43 36 62 | FRAC | 43 24 63 | RCL FV | 45 15 64 | x | 20 65 | CFj | 43 14 66 | RCL CFj | 45,43 14 67 | RCL n | 45 11 68 | 0 | 0 69 | x<=y | 43 34 70 | GTO 54 | 43,33 54 71 | GTO 00 | 43,33 00
R/S - results are in (0) - (79) 8 digits per register no leading zeros - takes about 45 hours.This program will show register (0) when it ends and start a continuous loop thru all the registers, advanced each time R/S is pressed. IMPORTANT NOTE: If the calculator times-out and shuts down before you get a chance to press R/S you must do a .GTO 081 after turning the calculator back on. Failure to do this will cause the program to re-run wiping out your result. This happens because the calculator resets the program pointer to 000 when it turns on.
based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits
001 | RPN | 42 16 002 | CLEAR FIN | 42 34 003 | 0 | 0 004 | CFj | 43 14 005 | 7 | 7 006 | 9 | 9 007 | RCL n | 45 11 008 | x<=y | 43 34 009 | GTO 003 | 43,33,003 010 | (F) 0 | 42 0 011 | 2 | 2 012 | 1 | 1 013 | 0 | 0 014 | 0 | 0 015 | STO i | 44 12 016 | EEX | 26 017 | 8 | 8 018 | STO FV | 44 15 019 | 0 | 0 020 | STO n | 44 11 021 | STO PMT | 44 14 022 | RCL PMT | 45 14 023 | RCL FV | 45 15 024 | x | 20 025 | RCL CFj | 45,43 14 026 | RCL i | 45 12 027 | x | 20 028 | + | 40 029 | ENTER | 36 030 | ENTER | 36 031 | RCL i | 45 12 032 | 2 | 2 033 | x | 20 034 | 1 | 1 035 | + | 40 036 | STO PV | 44 13 037 | / | 10 038 | INTG | 43 25 039 | CFj | 43 14 040 | RCL PV | 45 13 041 | x | 20 042 | - | 30 043 | STO PMT | 44 14 044 | 7 | 7 045 | 9 | 9 046 | RCL n | 45 11 047 | 1 | 1 048 | + | 40 049 | STO n | 44 11 050 | x<=y | 43 34 051 | GTO 022 | 43,33,022 052 | 2 | 2 053 | STO 0 | 44 0 054 | RCL i | 45 12 055 | 1 | 1 056 | - | 30 057 | STO i | 44 12 058 | x=0 | 43 35 059 | GTO 061 | 43,33,061 060 | GTO 019 | 43,33,019 061 | 7 | 7 062 | 9 | 9 063 | STO n | 44 11 064 | RCL i | 45 12 065 | RCL CFj | 45,43 14 066 | + | 40 067 | RCL FV | 45 15 068 | / | 10 069 | INTG | 43 25 070 | STO i | 44 12 071 | LASTx | 43 40 072 | FRAC | 43 24 073 | RCL FV | 45 15 074 | x | 20 075 | CFj | 43 14 076 | RCL CFj | 45,43 14 077 | RCL n | 45 11 078 | 0 | 0 079 | x<=y | 43 34 080 | GTO 064 | 43,33,064 081 | 0 | 0 082 | STO n | 44 11 083 | RCL CFj | 45,43 14 084 | R/S | 31 085 | 7 | 7 086 | 9 | 9 087 | RCL n | 45 11 088 | 2 | 2 089 | + | 40 090 | x<=y | 43 34 091 | GTO 082 | 43,33,082 092 | GTO 081 | 43,33,081
GSB 0 - results are printed are in (2) - (29), 6 digits per register no leading zeros - takes about 8.5 hours. Also, a progress countdown is given showing the iteration counter. After a few minutes you'll see 53 print than 52, 51, ... every 9.5 minutes or so.
based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits
01 | *LBL 0 | 25 14 00 02 | FIX0 | 16 13 00 03 | CLRG | 16 23 04 | 5 | 05 05 | 3 | 03 06 | 8 | 08 ; 538 iterations = 162 digits, convergence is about 1 bit per iteration 07 | STO1 | 45 01
08 | *LBL1 | 25 14 01 ; top of main loop
09 | RCL1 | 55 01 10 | 1 | 01 11 | 0 | 00 12 | / | 61 13 | INT | 16 52 14 | LSTX | 16 63 15 | X=Y? | 16 61 16 | PRTX | 65 ; print status every 10 iterations
17 | 2 | 02 18 | STO0 | 45 00 19 | 0 | 00 ; initial value of carry
20 | *LBL2 | 25 14 02 ; top of register loop 21 | EEX | 23 22 | 6 | 06 23 | x | 51 24 | RCLi | 55 12 25 | RCL1 | 55 01 26 | x | 51 27 | + | 41 28 | ENT^ | 21 29 | ENT^ | 21 30 | ENT^ | 21 31 | RCL1 | 55 01 32 | 2 | 02 33 | x | 51 34 | 1 | 01 35 | + | 41 36 | / | 61 37 | LSTX | 16 63 38 | X<>Y | 11 39 | INT | 16 52 40 | STOi | 45 12 41 | x | 51 42 | - | 31 43 | ISZ | 25 55 44 | RCL0 | 55 00 45 | 3 | 03 46 | 0 | 00 47 | - | 31 48 | X=0? | 25 61 49 | GTO3 | 14 03 50 | Rv | 12 51 | GTO2 | 14 02
52 | *LBL3 | 25 14 03 53 | 2 | 02 54 | STO2 | 45 02 55 | RCL1 | 55 01 56 | 1 | 01 57 | - | 31 58 | STO1 | 45 01 59 | X>0? | 25 41 60 | GTO1 | 14 01
61 | 2 | 02 ; loop through registers and carry overflow up 62 | 9 | 09 63 | STO0 | 45 00 64 | 0 | 00
65 | *LBL4 | 25 14 04 66 | ST+i | 45 41 12 67 | RCLi | 55 12 68 | EEX | 23 69 | 6 | 06 70 | / | 61 71 | INT | 16 52 72 | STO1 | 45 01 73 | EEX | 23 74 | 6 | 06 75 | x | 51 76 | ST-i | 45 31 12 77 | RCL1 | 55 01 78 | DSZ | 25 45 79 | GTO4 | 14 04
80 | 2 | 02 ; print results in registers (2) through (29) 81 | STO0 | 45 00
82 | *LBL5 | 25 14 05 83 | RCLi | 55 12 84 | PRTX | 65 85 | 1 | 01 86 | ST+0 | 45 41 00 87 | RCL0 | 55 00 88 | 3 | 03 89 | 0 | 00 90 | X>Y? | 16 41 91 | GOT5 | 14 05 92 | RTN | 25 13 93 | R/S | 64
Put 99 (maximum number, you can use less) on the stack and run.
- results are in Data(X) and Data(Y), 7 digits per register no leading zeros - takes about 39 minutes (no other programs in the calc) - a progress countdown is given showing the iteration counter
based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits 1 STO 9 -- Save the number of registers to fill
2 RCL Mode -- save the Mode 3 STO 5 4 3 5 0 6 0 7 STO Mode -- Set RPN and FIX 00
8 Reset -- Clear the CF registers to make maximum room in the Stats Registers 9 DOWN 10 DOWN 11 INPUT 12 INPUT 13 Reset 14 DOWN 15 DOWN 16 DOWN 17 INPUT 18 INPUT
19 RCL 9 -- Figure out how many iterations will be needed based on the number of registers to fill 20 7 -- 7 digits per register 21 * 22 2 23 * 24 2 25 LN 26 1 27 0 28 LN 29 / 30 / - # digits x log-base-2(10) are the iterations needed (and we count down by 2)
31 MATH -- take the IP 32 UP 33 UP 34 INPUT 35 STO 1 -- and save in 1
36 LBL 80 - top of main loop 37 RCL 9 38 EEX 39 3 40 / 41 STO 0 -- set up number of register loop for ISG 42 2 -- initial value of term 43 STO Data 44 0 -- initial carry 45 INPUT
46 LBL 81 -- top of loop through the Stat registers 47 EEX 48 7 49 * 50 RCL DATA 51 RCL x 1 52 + -- numerator is Data(i)*n + carry * 10^7 53 STO 2 54 RCL 1 55 Disp0 56 2 57 * 58 1 59 + -- denominator is 2n+1 60 STO 3 61 / 62 MATH -- take the IP 63 UP 64 UP 65 INPUT 66 STO Data -- IP((Data(i)*n+carry*10^7)/(2n+1)) --> Data(i) 67 RCL x 3 68 RCL 2 69 SWAP 70 - -- carry into next register 71 ISG 0 -- check for end of Stat register loop 72 GTO 81
73 0 -- next 4 lines are needed because of the way the loop ends, we need "2" in Data(0) 74 STO 0 75 2 76 STO Data 77 DSE 1 -- check for end of term loop 78 GTO 80
79 RCL 9 -- This is the final "fix up" loop, just one pass through the registers to adjust the overflows 80 STO 0 81 0 -- initial value of carry 82 INPUT 83 LBL 82 84 RCL Data 85 + 86 STO Data -- Add carry to register 87 EEX 88 7 89 ?< 90 GT 83 -- no overflow so skip over adjustment 91 RCL Data 92 EEX 93 7 94 - -- back out overflow 95 STO Data 96 1 -- carry into next register 97 INPUT 98 GTO 84 99 LBL 83 100 0 -- no carry 101 INPUT 102 LBL 84 103 DSE 0 -- decrement register pointer and loop 104 GTO 82 105 RCL Data -- add carry to Data(0) 106 + 107 STO Data
108 RCL 5 -- done, restore mode settings 109 STO Mode 110 Data -- leave the user in Data review mode to see the results 111 Stop
--------------------------------------------- 157 Bytes 014 Checksum
GSB 0 - results in (2) - (53), 6 digits per register no leading zeros. It takes about 43 hours to run on an original 15C or 16.1 minutes on a 15C LE.based on: ------------------------------------------------------------------------- pi = 2 x sigma[i,1,1,log2(10)*N,i/(2*i+1)] = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + 4/9 * (2 + ......
where N = number of decimal digits
01 | *LBL 0 | 42,21, 0 02 | 5 | 5 03 | 3 | 3 04 | DIM (i) | 42,23,24 05 | FIX 0 | 42, 7, 0 06 | CLR REG | 42 34 07 | 1 | 1 08 | 0 | 0 09 | 1 | 1 10 | 7 | 7 ; 1017 iterations = 306 digits, convergence is about 1 bit per iteration 11 | STO 1 | 44 1
12 | *LBL 1 | 42,21, 1 ; top of main loop 13 | 2 | 2 14 | STO I | 44 25 15 | 0 | 0 ; initial value of carry
16 | *LBL 2 | 42,21, 2 ; top of register loop 17 | EEX | 26 18 | 6 | 6 19 | x | 20 20 | RCL (i) | 45 24 21 | RCL 1 | 45 1 22 | x | 20 23 | + | 40 24 | ENTER | 36 25 | ENTER | 36 26 | ENTER | 36 27 | RCL 1 | 45 1 28 | 2 | 2 29 | x | 20 30 | 1 | 1 31 | + | 40 32 | / | 10 33 | LST x | 43 36 34 | X<>Y | 34 35 | INT | 43 44 36 | STO (i) | 44 24 37 | x | 20 38 | - | 30 39 | 1 | 1 40 | STO+I | 44,40,25 41 | Rv | 33 42 | RCL I | 45 25 43 | 5 | 5 44 | 4 | 4 45 | - | 30 46 | X=0? | 43 20 47 | GTO 3 | 22 3 48 | Rv | 33 49 | GTO 2 | 22 2
50 | *LBL 3 | 42,21, 3 51 | 2 | 2 52 | STO 2 | 44 2 53 | RCL 1 | 45 1 54 | 1 | 1 55 | - | 30 56 | STO 1 | 44 1 57 | X>0? | 43,30, 1 (TEST 1) 58 | GTO 1 | 22 1
59 | 5 | 5 ; loop through registers and carry overflow up 60 | 3 | 3 61 | STO I | 44 25 62 | 0 | 0
63 | *LBL 4 | 42,21, 4 64 | STO+(i) | 44,40,24 65 | RCL (i) | 45 24 66 | EEX | 26 67 | 6 | 6 68 | / | 10 69 | INT | 43 44 70 | STO 1 | 44 1 71 | EEX | 26 72 | 6 | 6 73 | x | 20 74 | STO-(i) | 44,30,24 75 | RCL 1 | 45 1 76 | DSE I | 42, 5,25 77 | GTO 4 | 22 4
78 | RCL 2 | 45 2 79 | RTN | 43 32
Edited: 20 Jan 2012, 3:18 p.m.
Go back to the main exhibit hall