Calculation of pi on many machines

[Benoit Maag]

Not a new subject by a long shot, the attached PDF lists programs to calculate pi with many decimals on many machines and benchmarking their relative speeds (see graph at the end). This is an update of a previous article, adding the 41CL, WP 34S, DM 42, DM 41X to the benchmark. It is incredible what these little machines can do...

Includes:
- HP-41C, 41CL
- DM42, DM41X
- WP 34S
- HP-15C-LE
- HP-65 (yes, it is possible!)
- HP-29C
- HP-67
- Casio fx-6300G and fx-7000GA

And if you have an HP-71B, Valentin has a 6 line program (!) at:

https://albillo.hpcalc.org/programs/HP%2...20time.pdf

[rprosperi]

Thanks for all your time doing the testing, documenting it all and sharing it here Benoit.

There's LOTS of real useful and interesting stuff just in the chart alone. Among them:

a. the 65 is faster the 67
b. the 41 is only a tad faster than the earlier machines and the 29 is faster still
c. the surprisingly high performance of the early Casio machines
d. the amazing speed of the DM42
e. the quite different slope of the 65 and 29 (gotta go read that in detail)

Looking forward to reading it in more detail.

Thanks!

[Gene]

The 29 is slower than the 41. That's the 29-LP with the upgrade.

The 41 has no time shown for 42 digits so the graph point is an interpolation assuming a linear time - hard to know if the HP-65 is faster without a direct time for the 41.

[Gerson W. Barbosa]

7 seconds for 4900 digits on the iPhone 7 running Ron Knapp’s program (Free42).

Code:


00 { 629-Byte Prgm }
01▸LBL "PI"
02 CLRG
03 CF 29
04 FIX 00
05 9
06 +
07 10
08 ÷
09 IP
10 2
11 ×
12 5
13 +
14 ARCL IND ST X
15 5
16 +
17 2
18 ÷
19 STO 02
20 1.07
21 ×
22 5
23 -
24 RND
25 4
26 ×
27 1
28 -
29 4ᴇ-5
30 +
31 STO 01
32▸LBL 00
33 RCL 01
34 IP
35 ENTER
36 ENTER
37 ENTER
38 4.3
39 ÷
40 IP
41 1
42 X>Y?
43 X<>Y
44 R↓
45 +/-
46 RCL 02
47 +
48 1ᴇ3
49 ÷
50 4
51 +
52 STO 00
53 STO 03
54 R↓
55 239
56 X↑2
57 ×
58 X<>Y
59 2
60 -
61 STO× ST Z
62 -
63 0
64 X<>Y
65 4
66 ×
67 XEQ 03
68 239
69 X↑2
70 X↑2
71 ENTER
72 ENTER
73 RCL 03
74 STO 00
75 CLX
76 DSE 01
77 GTO 04
78 X<>Y
79 239
80 ÷
81 X<>Y
82 XEQ 06
83 GTO 08
84▸LBL 04
85 XEQ 06
86 GTO 00
87▸LBL 02
88 1ᴇ5
89 ×
90 ENTER
91 R↑
92 STO÷ ST Z
93 MOD
94 X<>Y
95 IP
96 1ᴇ5
97 STO× ST Z
98 ×
99 X<>Y
100▸LBL 03
101 ENTER
102 R↑
103 STO÷ ST Z
104 MOD
105 LASTX
106 X<>Y
107 R↑
108 R↑
109 IP
110 +
111 1ᴇ10
112 -
113 RCL IND 00
114 +
115 X<0?
116 GTO 04
117 1
118 STO- 00
119 STO+ IND 00
120 STO+ 00
121 CLX
122 +
123 GTO 05
124▸LBL 04
125 1ᴇ10
126 +
127▸LBL 05
128 STO IND 00
129 CLX
130 +
131 ISG 00
132 GTO 02
133 RTN
134▸LBL 06
135 1ᴇ5
136 ×
137 ENTER
138 R↑
139 STO÷ ST Z
140 MOD
141 X<>Y
142 IP
143 1ᴇ5
144 X>Y?
145 GTO 07
146 ÷
147 IP
148 1
149 STO- 00
150 X<>Y
151 STO+ IND 00
152 R↓
153 STO+ 00
154 CLX
155 LASTX
156 FP
157 1ᴇ5
158 ×
159 LASTX
160▸LBL 07
161 ×
162 X<> IND 00
163 LASTX
164 ÷
165 IP
166 STO+ ST Y
167 X<> ST L
168 FP
169 X<>Y
170 1ᴇ5
171 STO× ST Z
172 ×
173 ENTER
174 R↑
175 STO÷ ST Z
176 MOD
177 LASTX
178 R↓
179 X<>Y
180 IP
181 STO+ IND 00
182 CLX
183 +
184 +
185 ISG 00
186 GTO 06
187 RTN
188▸LBL 08
189 SF 29
190 RCL 02
191 RCL 03
192 IP
193 -
194 STO+ 02
195 1.2
196 ×
197 0.9
198 -
199 RND
200 12
201 ×
202 1
203 -
204 6ᴇ-5
205 +
206 STO 01
207▸LBL 09
208 RCL 01
209 IP
210 ENTER
211 14.3
212 ÷
213 IP
214 1
215 X>Y?
216 X<>Y
217 R↓
218 +/-
219 RCL 02
220 +
221 LASTX
222 2
223 STO÷ ST Y
224 +
225 X<>Y
226 1ᴇ3
227 ÷
228 +
229 STO 00
230 STO 03
231 SF 00
232▸LBL 10
233 RCL 01
234 IP
235 ENTER
236 ENTER
237 ENTER
238 2
239 -
240 STO× ST Z
241 25
242 STO× ST Z
243 X<>Y
244 ×
245 2
246 STO- ST L
247 CLX
248 LASTX
249 STO× ST T
250 STO- ST Y
251 R↓
252 ×
253 R↑
254 STO+ ST T
255 X↑2
256 R↑
257 +
258 +
259 FS? 00
260 GTO 04
261 5ᴇ6
262 ×
263 0.75
264 GTO 05
265▸LBL 04
266 DSE 01
267 40
268 ×
269 6
270▸LBL 05
271 STO÷ ST Z
272 CLX
273 X<>Y
274 XEQ 12
275 RCL 03
276 STO 00
277 FS?C 00
278 GTO 10
279 125
280 X↑2
281 X↑2
282 ENTER
283 ENTER
284 0
285 DSE 01
286 GTO 07
287 3
288 STO× ST Z
289 STO× ST T
290 R↓
291 XEQ 06
292 GTO 13
293▸LBL 07
294 XEQ 06
295 GTO 09
296▸LBL 11
297 ENTER
298 R↑
299 STO÷ ST Z
300 MOD
301 X<>Y
302 IP
303 X<>Y
304 1ᴇ5
305 STO× ST Z
306 ×
307▸LBL 12
308 ENTER
309 R↑
310 STO÷ ST Z
311 MOD
312 X<>Y
313 LASTX
314 R↓
315 IP
316 STO+ ST Z
317 CLX
318 +
319 1ᴇ5
320 ×
321 X<>Y
322 RCL IND 00
323 -
324 X<0?
325 X>0?
326 GTO 04
327 1
328 STO- 00
329 STO- IND 00
330 STO+ 00
331 CLX
332 1ᴇ10
333 +
334▸LBL 04
335 STO IND 00
336 R↓
337 ISG 00
338 GTO 11
339 RTN
340▸LBL 13
341 RCL 03
342 FP
343 LASTX
344 IP
345 1
346 -
347 STO 00
348 1ᴇ3
349 STO× ST Z
350 ÷
351 STO 01
352 +
353 STO 03
354 CLX
355▸LBL 14
356 RCL IND 03
357 X<>Y
358 -
359 RCL IND 00
360 -
361 0
362 X<>Y
363 X<0?
364 X>0?
365 GTO 05
366 1ᴇ10
367 +
368 1
369 X<>Y
370▸LBL 05
371 STO IND 00
372 R↓
373 DSE 00
374 DSE 03
375 GTO 14
376▸LBL 01
377 RCL 01
378 4
379 +
380 STO 00
381 VIEW IND 00
382 PSE
383 ISG 00
384 FIX 05
385 BEEP
386▸LBL 15
387 RCL IND 00
388 1ᴇ5
389 ÷
390 VIEW ST X
391 ISG 00
392 GTO 15
393 RTN
394 GTO 01
395 END

[rprosperi]

(11-22-2020 11:02 PM)Gene Wrote:  The 29 is slower than the 41. That's the 29-LP with the upgrade.

Thanks Gene. That explains why I did a double-take about the 29 being faster, I knew I had seen it slower, but then I also saw it faster and thought I recalled wrong. Well, 29C vs. 29LP, so it was wrong, kinda...

[Benoit Maag]

Bob,

To answer some of your questions:

a. The comparison is not quite apples to apples. For most of the machines (including the HP-67), the data is stored as #####.##### (5 + 5 digits) and the calculation is performed on 5 digits at a time. This means each registered has to be split into 5 digits which are processed separately. But this is not possible on the HP-65 because there is not enough program memory and no indirect addressing. So the data is stored as 0.####### (7 digits), there is no need to split registers and this is why it is faster than the HP-67. It would be interesting to see how the 67 performs with the 65 program

c. I was surprised by the Casio as well, especially the fox-6300G which is a cheap machine. I think the speed is attributable to the very simple programming language

d. The DM42 is super fast because the data is stored as 17+17 digits and all calculations are made 17 digits at a time. I have not checked the performance with a 5+5 digit program. BTW, for the WP 34S, I tested both 5+5 and 8+8 digit data

e. I am not sure why the slopes would be different. After all, the time is proportional to the square of the number of desired digits so all curves should be parallel

More to come as I would like to run the benchmark on the 32sii, 35s, RPL machines, BASIC sharps and casios, NumWorks (python), TI-58, and maybe the RPN-1250 I presented at HHC 2018, though it is not a bcd machine and has poor calculation precision

(11-22-2020 10:50 PM)rprosperi Wrote:  Thanks for all your time doing the testing, documenting it all and sharing it here Benoit.

There's LOTS of real useful and interesting stuff just in the chart alone. Among them:

a. the 65 is faster the 67
b. the 41 is only a tad faster than the earlier machines and the 29 is faster still
c. the surprisingly high performance of the early Casio machines
d. the amazing speed of the DM42
e. the quite different slope of the 65 and 29 (gotta go read that in detail)

Looking forward to reading it in more detail.

Thanks!

[Paul Dale]

I would expect that the WP 34S in double precision would be faster than in single precision -- 34 digits instead of 16 and most intermediate computations are carried in the high precision regardless.

I somewhat suspect that the greatest accuracy could be achieved in integer mode where there are no wasted bits in registers.


Pauli

[PeterP]

There are also some MCODE versions available for HP41 which you might or might not want to include...

https://www.hpmuseum.org/cgi-sys/cgiwrap...587#147587

Cheers

PeterP

[Gerson W. Barbosa]

Your DM-42 version is about three times faster than Ron Knapp’s (17/5 = 3.4), except that it doesn’t work here. I’ve double-checked the listing and compared it against the 41 version, but there must be something I may have overlooked.

Code:


00 { 188-Byte Prgm }
01▸LBL "PI"
02 ENTER
03 CLRG
04 17
05 +
06 2
07 STO 04
08 LOG
09 ÷
10 IP
11 STO 01
12 R↓
13 34
14 ÷
15 4
16 +
17 STO 00
18 STO 02
19 TIME
20 STO 03
21 0
22▸LBL 00
23 RCL IND 00
24 FP
25 XEQ 01
26 X<> IND 00
27 IP
28 1ᴇ17
29 ÷
30 XEQ 01
31 1ᴇ17
32 ×
33 STO+ IND 00
34 R↓
35 DSE 00
36 RCL 00
37 4
38 X>Y?
39 GTO 02
40 R↓
41 R↓
42 GTO 00
43▸LBL 01
44 RCL 01
45 ×
46 X<>Y
47 1ᴇ17
48 ÷
49 +
50 IP
51 LASTX
52 FP
53 RTN
54▸LBL 02
55 0
56 ISG 00
57 FIX 05
58▸LBL 05
59 1ᴇ17
60 ×
61 RCL IND 00
62 IP
63 +
64 XEQ 03
65 X<> IND 00
66 FP
67 +
68 1ᴇ17
69 ×
70 XEQ 03
71 1ᴇ17
72 ÷
73 STO+ IND 00
74 R↓
75 ISG 00
76 FIX 05
77 RCL 00
78 RCL 02
79 X<Y?
80 GTO 04
81 R↓
82 R↓
83 GTO 05
84▸LBL 03
85 ENTER
86 ENTER
87 RCL 01
88 STO+ ST X
89 ISG ST X
90 FIX 05
91 STO ST T
92 MOD
93 X<>Y
94 R↑
95 ÷
96 IP
97 RTN
98▸LBL 04
99 2
100 STO+ 04
101 DSE 01
102 GTO 06
103 TIME
104 RCL 03
105 HMS-
106 VIEW ST X
107 RTN
108▸LBL 06
109 VIEW 01
110 DSE 00
111 0
112 GTO 00
113 END

Edited to fix a typo and a phrase.

[John Keith]

(11-23-2020 04:51 AM)Benoit Maag Wrote:  More to come as I would like to run the benchmark on the 32sii, 35s, RPL machines, ...

A note as to RPL machines: the LongFloat library for the 49/50 can calculate pi up to 9999 digits. I would still be interested in seeing how a userRPL program would compare in speed.

[Werner]

My current HP-41 compatible version runs in less than a minute (58s) on the DM-42 on USB, producing 1017 correct digits, using
Pi = 16*ATAN(1/5) - 4*ATAN(1/239)
It runs unchanged on a 41, 42 and Free42/DM42, but will use all available precision, ie. 10 digits on a 41, 12 on a 42 and 34 on Free42/DM42.
It is written so that it can use any ATAN formula X*ATAN(N/D), as long as N^2+D^2 is less than 'half-length', ie 10^5 for a 41, 10^6 for a 42S and 10^17 for Free42/DM42.
Using another ATAN formula, I can get the DM42 USB time under 38 seconds (but it won't run on a 41 or 42 any more).

Cheers, Werner

[Gerson W. Barbosa]

(11-23-2020 06:38 PM)John Keith Wrote:  

(11-23-2020 04:51 AM)Benoit Maag Wrote:  More to come as I would like to run the benchmark on the 32sii, 35s, RPL machines, ...

A note as to RPL machines: the LongFloat library for the 49/50 can calculate pi up to 9999 digits. I would still be interested in seeing how a userRPL program would compare in speed.

Code:



DIR
  PiM
  « PUSH RAD -105 CF -3 CF DUP 2 / 5 DUP2 LOG / →NUM CEIL 2 * 1 - R→I ATN1X 
16 * 239 ROT OVER LOG / →NUM CEIL 2 * 1 - R→I ATN1X 4 * - EXPAND FXND DUP
 SIZE R→I ALOG OVER - PICK3 * SWAP IQUOT + →STR DUP HEAD 0 I→R →STR TAIL
 + SWAP TAIL + 1 ROT 2 + SUB POP
  »
  ATN1X
  « -1 UNROT 1 4 ROLLD OVER SQ PICK3 INV 4 ROLLD UNROT 3 SWAP
    FOR i OVER * DUP i * INV 5 PICK * 4 ROLL + UNROT 5 ROLL 4 ROLLD 2
    STEP DROP2 UNROT DROP2
  »
END

50 digits: 24.70 s
100 digits: 136.59 s

( HP 50g )

Unpractical running times henceforward, but then again my implementation might be far from optimal. EXPAND takes too much time for long expressions.

P.S.: If you prefer resorting to an external library, a LongFloat version is available here.

707 digits in 1858 seconds or roughly half an hour.

[John Keith]

On EMU48 on my phone (Galaxy S9+):
<< 4900 DIGITS STO FPI
>>
8.876 seconds.

[Massimo Gnerucci]

On my Core I7 with PiFast 4.3:
2000000 digits in 3.08"

:)

[Gerson W. Barbosa]

On my HP-42S, using my code and my formula:
12 digits in 1.7"

:-)


00 { 55-Byte Prgm }
01▸LBL "π"
02 4
03 0
04 8
05 RCL× ST T
06 RCL+ ST Z
07▸LBL 00
08 R↑
09 RCL+ ST X
10 STO× ST X
11 STO× ST T
12 DSE ST X
13 STO÷ ST T
14 X<>Y
15 STO+ ST Z
16 X<> ST Z
17 STO÷ ST Y
18 X<> ST L
19 R↓
20 X<>Y
21 DSE ST T
22 GTO 00
23 DSE ST X
24 +
25 1/X
26 0.5
27 +
28 ×
29 END


5 XEQ π →

3.14159265358


\(\pi \approx \left ( \frac{4}{3} \times \frac{16}{15}\times \frac{36}{35}\times\frac{64}{63} \times \cdots \times \frac{ 4n ^{2}}{ 4n ^{2}-1}\right ) \left ( 2+\frac{4}{8n+3+\frac{3}{8n+4+\frac{15}{8n+4+ \frac{35}{8n+4 + \frac{63}{\dots\frac{\ddots }{8n+4+\frac{4n^{2}-1}{8n+4}}} }}} } \right )\)

[pier4r]

(11-24-2020 01:49 AM)John Keith Wrote:  On EMU48 on my phone (Galaxy S9+):
<< 4900 DIGITS STO FPI
>>
8.876 seconds.

FPI ? is that a builtin command?

[Gerson W. Barbosa]

(11-24-2020 06:29 PM)pier4r Wrote:  

(11-24-2020 01:49 AM)John Keith Wrote:  On EMU48 on my phone (Galaxy S9+):
<< 4900 DIGITS STO FPI
>>
8.876 seconds.

FPI ? is that a builtin command?

No, it’s from the LongFloat library, the first 1000 digits being hard-coded:

1000 digits in 0.066 s, 1100 in 51.710 s. (Real HP 50g)

[Gerson W. Barbosa]

This is an RPL version of the HP-42S program above:

Code:


« PUSH RAD -105 CF -3 CF DUP 10 * 21 / 1 + CEIL 0 2 ROT 8 OVER * 
PICK3 + UNROT 1
  FOR i i DUP + DUP 1 - UNROT * OVER / SWAP SQ 4 ROLL 4 PICK + / 
UNROT -1
  STEP SQ UNROT 2 / + / EXPAND FXND DUP SIZE R→I ALOG OVER - 
PICK3 * SWAP IQUOT + →STR DUP HEAD 0 I→R →STR TAIL + SWAP TAIL 
+ 1 ROT 2 + SUB POP
»

# 8F62h, 258 bytes

51 decimal digits: 10.502 s
85 decimal digits: 32.557 s
153 decimal digits: 160.952 s
357 decimal digits: 4393.037 s

(HP 50g)

———

PS:

Actually, the RPL program above is based on one of my older formulae. The following program is slower, but it is in accordance with the formula I used in the RPN program:

Code:


« PUSH RAD -105 CF -3 CF R→I DUP 12 * 25 / CEIL DUPDUP 8 * 4 + DUP 
8 - ROT SQ 4 * 1 - 2 0 UNROT 1 7 ROLL
  START OVER 1 + * OVER / UNROT DUP ROT 6 PICK + / UNROT 4 ROLL 
SWAP OVER - UNROT 8 - 4 ROLLD
  NEXT NIP 2 ROT 5 ROLL + 1 - / 1 + * NIP EXPAND FXND DUP SIZE R→I 
ALOG OVER - PICK3 * SWAP IQUOT + →STR DUP HEAD 0 I→R →STR TAIL 
+ SWAP TAIL + 1 ROT 2 + SUB POP
»

# 5CB1h, 301.5 bytes

51 decimal digits: 12.344 s
85 decimal digits: 42.373 s
153 decimal digits: 308.168 s

(HP 50g)

1000 decimal digits: 2988.379 s

(iPhone7, emulated HP 49g+, m48+)

Algorithm:

Code:


input digits ; number of decimal digits 
n := ceil(12*digits/25)
c := 0
d := 8*n + 4
e := d - 8
m := 4*n*n - 1
w := 2
i := n
repeat
  w := w*(m + 1)/m
  c := m/(c + d)
  m := m - e
  e := e - 8
  i := i - 1
until i = 0
c := 2/(c + d - 1) + 1
p := c*w
output p

———

PPS:

This is a longer version of the previous program, but it’s significantly faster for larger numbers of decimal digits (and even faster than the first program). Also, as in the previous program, all decimal digits are correct (actual digits, no roundings).

Code:


« PUSH RAD -105 CF -3 CF R→I DUP 12 * 25 / CEIL DUPDUP DUP + 1 + 
OVER 8 * 4 + DUP 8 - 0 5 ROLL SQ 4 * 1 - 5 PICK 1 DUP 9 ROLL 
  START OVER 2 - UNROT OVER 1 - ROT / * PICK3 6 PICK - UNROT 4 ROLL 
5 ROLL 7 PICK + / 4 ROLLD 5 ROLL 8 - 5 ROLLD
  NEXT SQ 7 ROLL * 4 5 ROLL 7 ROLL + 1 - / 2 + * 4 ROLLD 3 DROPN 
EXPAND FXND DUP SIZE R→I ALOG OVER - PICK3 * SWAP IQUOT + →STR 
DUP HEAD 0 I→R →STR TAIL + SWAP TAIL + 1 ROT 2 + SUB POP
»

# E567h, 354 bytes

51 decimal digits: 10.355 s
85 decimal digits: 31.200 s
153 decimal digits: 157.039 s
357 decimal digits: 3820.207 s

(HP 50g)

1000 decimal digits: 2588.102 s

(iPhone7, emulated HP 49g+, m48+)

Algorithm:

Code:


input digits ; number of decimal digits 
n := ceil(12*digits/25)
a := 2*n + 1
b := a
c := 0
d := 8*n + 4
e := d - 8
m := 4*n*n - 1
w := 1
i := n
repeat
  w := w*(a - 1)/a
  c := m/(c + d)
  a := a - 2
  m := m - e
  e := e - 8
  i := i - 1
until i = 0
c := 4/(c + d - 1) + 2
p := b*c*w*w
output p

[Gerson W. Barbosa]

The Science et Vie formula makes for very short programs, but they are too slow.

Code:


« PUSH RAD -105 CF -3 CF 2 OVER 2 LOG / CEIL DUPDUP + 1 + UNROT 1
  FOR i i * SWAP DUP 2 - UNROT / 2 + -1
  STEP NIP EXPAND FXND DUP SIZE R→I ALOG OVER - PICK3 * SWAP IQUOT + 
 →STR DUP HEAD 0 I→R →STR TAIL + SWAP TAIL + 1 ROT 2 + SUB POP
»

# FE89h, 200 bytes.

51 decimal digits: 634.528 s.

Most of the slowness is due to the difficulty of EXPAND to process the long resulting algebraic expression, however. For benchmarking purposes it is better to take a more conventional approach.

[berndpr]

Dear Benoit!

Quote: And if you have an HP-71B, Valentin has a 6 line program (!) at:

This spigot program for pi is not the best and fastet algorithm and needs a lot of storage.
But it's very short.

I have used the same references to make a similar program and it's nearly identical to Valentin's.

But I got another program from a japanese person (Hiro510) which is using different version of pi/4 approximations (like Leibniz, Machin, Klingenstirna, Gauss, Störmer, Takano Kikuo and 2 unknown).
Some are very fast like Kikuo's and some are extremly slow like Leibniz's.

But if you are interested, the calculation of e (Euler constant) with a similar spigot program is much more faster and need much less storage.

Bye
Berndpr