Gerson's Pi Program

[EdS2]

Did anyone else have a go at de-obfuscating and running Gerson's short pi-calculating program from the festive greetings thread?

(12-25-2020 11:45 PM)Gerson W. Barbosa Wrote:  

I'd be interested to know how it works!

I made a plain copy in BBC Basic, which you can run in your browser here.

And here it is below anti-spoiler scroll filler:

Code:


.
..
.
...
.
..
.
....
.
..
.
...
.
..
.

10 INPUT n
20 c=0
30 d=8*n+4
40 e=d-8
50 m=4*n*n - 1
60 w=2
70 FOR i=1 TO n
80 w=w+w/m
90 c=m/(c+d)
100 m=m-e
110 e=e-8
120 NEXT
130 PRINT w*(2/(c+d-1)+1)

[Allen]

Quote:I'd be interested to know how it works!

Code:


50 m=4*n*n - 1
60 w=2
...
80 w=w+w/m

Looks like a nicely compact convergence based on Wallis' Product

[Gerson W. Barbosa]

(12-28-2020 11:20 AM)EdS2 Wrote:  I'd be interested to know how it works!

Hello, Ed!

Thank you very much for your interest! That’s an MSX BASIC listing. The MSX computer was very popular in Brazil and England in the mid 80’s (but not in the USA). That was my second computer. I liked its 14 significant digits, much better than my first computer, an Apple II clone.

The program is based on this formula:

\(\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 )\)


That’s the third Wallis-Wasicki formula, as I called it, rather jokingly. Notice the reuse of the denominators in the Wallis product as the numerators in the continued fraction.

Here’s the thread in which I presented it, in case you missed it:

https://www.hpmuseum.org/forum/post-1394...#pid139434

Best regards,

Gerson.

[EdS2]

Thanks! I will surely have seen that post before, but evidently failed to recognise the approach. I think I might blame my ever-increasing age.

[Gerson W. Barbosa]

(12-29-2020 09:15 AM)EdS2 Wrote:  Thanks! I will surely have seen that post before, but evidently failed to recognise the approach. I think I might blame my ever-increasing age.

My pleasure! Time goes forward for us all. I have to update my avatar so that it reflects my current age, but I don't remember how I did that twenty years ago.

Actually, this is the post where I first presented that formula, along with an explanation how I found it (not a proof, though).

https://www.hpmuseum.org/forum/post-1345...#pid134502

It appears to produce 25*n/12 correct digits (I have tested the algorithm to 1000 digits only).

The previous program is optimized for speed as it avoids two extra multiplications per iteration. For the few digits provided by most computer languages, it should be better to use this even more compact program:

Code:


10 INPUT N
15 C=0
20 D=8*N+4
25 W=2
30 FOR I=N TO 1 STEP -1
35 T=4*I*I-1
40 W=W+W/T
45 C=T/(C+D)
50 NEXT I
55 PRINT W*(2/(C+D-1)+1)

[Gerson W. Barbosa]

Here is an arbitrary precision version for the MSX computer. It shouldn't be difficult to convert it to the HP-BASIC (HP-75C and HP-71B), but the base should be lowered to 1E6 or 1E5.

Code:


10 REM  ***  PI BY WALLIS-WASICKI FORMULA  ***
20 CLS: KEYOFF: WIDTH 40
30 DEF FNFRAC(X) = X-INT(X)
40 DEF FNRMDR(X,Y)=X-Y*INT(X/Y)
50 DEFINT D,E,G-N
60 DIM A(50), B(50), C(100), Q(50), X(50), Y(100), W(50)
70 LINE INPUT "Number of decimal digits: ";N$: ND=VAL(N$): PRINT
80 B=10000000!: N=ND\7+2: K=ND*12\25+1
90 TIME=0
100 D=8*K+4: E=D-8: G=4*K*K-1
110 W(1)=4: B(1)=0
120 FOR I=2 TO N
130   W(I)=0: B(I)=0
140 NEXT I
150 FOR IT=1 TO K
160   FOR I=0 TO N
170     Q(I)=W(I)
180   NEXT I
190   GOSUB 860
200   GOSUB 930
210   B(1)=B(1)+D
220   GOSUB 420
230   GOSUB 1000
240   FOR I=1 TO N: B(I)=C(I): NEXT I
250   G=G-E
260   E=E-8
270 NEXT IT
280 B(1)=B(1)+D-1
290 GOSUB 420
300 FOR I=1 TO N
310   A(I)=W(I): X(I)=C(I)
320 NEXT I
330 X(2)=X(2)+B/2
340 GOSUB 670
350 TM = TIME/60
360 FOR I=1 TO N
370   PRINT USING"########";C(I);
380 NEXT I
390 PRINT: PRINT: PRINT USING"####.##";TM;: PRINT "s" : PRINT
400 END
410 REM  ***   LINV   ***
420 T=1/(B(1)+B(2)/B)
430 X(1)=INT(T): X(2)=INT(B*FNFRAC(T))
440 M=10*LOG(ND)/LOG(1000)-1
450 FOR L=1 TO M
460   FOR I=1 TO N
470     A(I)=B(I)
480   NEXT I
490   GOSUB 670
500   BR%=0
510   FOR I=N TO 2 STEP -1
520     T=FNRMDR(B-BR%-C(I),B)
530     BR%=BR% OR SGN(C(I))
540     Y(I)=T
550   NEXT I
560   Y(1)=2-BR%-C(1)
570   FOR I=1 TO N
580     A(I)=Y(I)
590   NEXT I
600   GOSUB 670
610   FOR I=1 TO N
620     X(I)=C(I)
630   NEXT I
640 NEXT L
650 RETURN
660 REM  ***   LMULT   ***
670 FOR I=1 TO N
680   C(I)=0
690 NEXT I
700 FOR I=N TO 1 STEP -1
710   TP=0: AC=0
720   FOR J=N-I+2 TO -1 STEP -1
730     CT=(TP+A(I)*X(J+1))/B
740     AC=FNFRAC(CT)*B
750     C(I+J)=C(I+J)+AC
760     TP=INT(CT)
770   NEXT J
780 NEXT I
790 FOR I=N+2 TO 2 STEP -1
800   T=C(I)/B
810   C(I)=FNFRAC(T)*B
820   C(I-1)=C(I-1)+INT(T)
830 NEXT I
840 RETURN
850 REM  ***   DIVIDE   ***
860 Y=0
870 FOR I= 1 TO N
880   X=W(I)+Y*B
890   W(I)=INT(X/G): Y=FNRMDR(X,G)
900 NEXT I
910 RETURN
920 REM  ***   ADD   ***
930 Y=0
940 FOR I=N TO 1 STEP -1
950   X=W(I)+Q(I)+Y
960   W(I)=FNRMDR(X,B): Y=INT(X/B)
970 NEXT I
980 RETURN
990 REM  ***   MULT   ***
1000 Y=0
1010 FOR I=N TO 1 STEP -1
1020   X=C(I)*G+Y
1030   C(I)=FNRMDR(X,B): Y=INT(X/B)
1040 NEXT I
1050 RETURN

Because full-precision multiplications and inversions are required in the main loop, this is not meant to be fast. It took eight minutes and eleven seconds on an emulator running ten times as fast compared to the real machine for just the first 100 digits. The program may return a few extra digits, but only the ones you ask for can be trusted.


Number of decimal digits: 100

       3 1415926 5358979 3238462 6433832
 7950288 4197169 3993751  582097 4944592  
 3078164  628620 8998628  348253 4211706  
 7982103

512.45s

Edited to fix a typo in line 720 and to update the new timings.

[EdS2]

Very nice! I'd like to port it to BBC Basic, where we have 4 byte integers and 5 byte floats. But it's not working... any ideas?

Minor notes:

The DEFINT are ignored, because in BBC Basic we need to decorate integer variables with % suffix. I was hoping that would merely be a speed advantage as the 5 byte floats are very nearly as good as the 4 byte ints for accuracy and range.

I took the ! off the end of the B=100000000 but I don't know what it means.

I changed the \ operators to DIV, guessing that this is what's meant.

I changed the printing to print an integer rounded version of the number.

Here's the link.

(Edit: typo for 4 byte ints)

[Gerson W. Barbosa]

(03-01-2022 09:39 AM)EdS2 Wrote:  Very nice! I'd like to port it to BBC Basic, where we have 4 byte integers and 5 byte floats. But it's not working... any ideas?

Minor notes:

The DEFINT are ignored, because in BBC Basic we need to decorate integer variables with % suffix. I was hoping that would merely be a speed advantage as the 5 byte floats are very nearly as good as the 5 byte ints for accuracy and range.

I took the ! off the end of the B=100000000 but I don't know what it means.

I changed the \ operators to DIV, guessing that this is what's meant.

B stands for base. Default variable types in MSX are double precision (8 bytes, 14 significant digits), which allows for bases as high as 10000000. With 9-digit precision, you should choose base = 10000 or perhaps 100000. The \ operator is the integer division, as you have correctly guessed. The exclamation sign says it's a single precision number.

Try changing lines 80 and 440:

80 B = 10000: N=ND DIV 4 + 2: K=ND*12 DIV 25 + 1

440 M=10*LOG(ND)/LOG(1000)

For 30 digits the output should be

3 1415 9265 3589 7932 3846 2643 3832 7941


P.S.: For a much faster program, you might want to try Katie Wasserman's program for the HP-71B:

CALCULATING MANY DIGITS OF PI

[EdS2]

Hmm, thanks for the suggested edits, but unfortunately something isn't quite right...
12 digits gives 3 1416 448 2514 8212

(Thanks also the link back to the earlier thread. I'll see about porting that too.)

(Sorry my link to owlet is not working well: something about url-escaping. My posted text works, the rendered link doesn't. If you open a draft which quotes my post you'll get the good link text.)

[Gerson W. Barbosa]

(03-02-2022 11:17 AM)EdS2 Wrote:  Hmm, thanks for the suggested edits, but unfortunately something isn't quite right...
12 digits gives 3 1416 448 2514 8212

I have a working version for the BeebEm emulator, but I don't know how to get a listing. Anyway, basically you should delete line 40 and use MOD instead:

520 T=(B-BR%-C(I)) MOD B

890 W(I)=INT(X/G): Y=X MOD G

960 W(I)=X MOD B: Y=INT(X/B)

1030 C(I)=X MOD B: Y=INT(X/B)

Also, lines 30, 740 and 740 should be changed to

30 DEF FNFRAC(X) = X+5E-6 - INT(X+5E-6)

740 AC=INT(B*FNFRAC(CT))

810 C(I)=INT(B*FNFRAC(T))

I've used only real variables, including the one in line 520, but you can change the types of the variables accordingly for speed and memory saving. On the MSX I cannot define AC and C() to integer, because its integer type is only 2-byte long. For the same reason MOD cannot be used on the MSX for bases greater than 10000.

For 12 digits you should get 3 1415 9265 3589 8142

[EdS2]

Thanks! With your changes, my suboptimal port produces the 12 correct digits as it should, in 53 seconds.

For reference, my simplistic and inefficient port of Katie's program produced 92 correct digits in 71 seconds.

[Gerson W. Barbosa]

(03-04-2022 04:50 PM)EdS2 Wrote:  With your changes, my suboptimal port produces the 12 correct digits as it should, in 53 seconds.

For reference, my simplistic and inefficient port of Katie's program produced 92 correct digits in 71 seconds.

As stated earlier, the WWF is not intended to be fast. According to an article by the Borwein Brothers in Scientific American (February 1988), not even 100 years of computing on a supercomputer programmed to evaluate the Wallis Product or the Gregory’s Series would yield 100 digits of π. So, getting 100 digits in 100 minutes on an 8-bit personal computer of that era using the Wallis Product as a basis is an unexpected surprise.

[Gerson W. Barbosa]

(03-04-2022 04:50 PM)EdS2 Wrote:  Thanks! With your changes, my suboptimal port produces the 12 correct digits as it should, in 53 seconds.

When porting the program to the HP-75C I noticed an error in line 720.

720 FOR J=N-1+2 TO -1 STEP -1

The first 1 should be an I:

720 FOR J=N-I+2 TO -1 STEP -1

Please correct it if you haven’t done so yet. This should make the program a little less slow.
I am away from home and have brought only the HP-75C along. Also, I can’t make your links work. Sorry for the mistake!

[Valentin Albillo]

.
Hi, Gerson,

(03-06-2022 01:25 AM)Gerson W. Barbosa Wrote:  According to an article by the Borwein Brothers in Scientific American (February 1988), not even 100 years of computing on a supercomputer programmed to evaluate the Wallis Product or the Gregory’s Series would yield 100 digits of π.

Perhaps, but then again there's this:

Speeding it up !

Regards.
V.
.
.

[Gerson W. Barbosa]

(03-08-2022 05:42 PM)Valentin Albillo Wrote:  

(03-06-2022 01:25 AM)Gerson W. Barbosa Wrote:  According to an article by the Borwein Brothers in Scientific American (February 1988), not even 100 years of computing on a supercomputer programmed to evaluate the Wallis Product or the Gregory’s Series would yield 100 digits of π.

Perhaps, but then again there's this:

Speeding it up !

Hello, Valentin,

Yes, I remember that particular HP-15C Mini-challenge of yours. It was linked by Thomas Klemm in the beginning of a thread I started. As the thread progressed, I ended up with what could have been another solution to your challenge, which I didn’t participate in:

https://www.hpmuseum.org/forum/post-9194.html#pid9194

This is a correction term to the Gregory’s Series in continued-fraction form, similar to the correction factor to the Wallis Product.

Best regards,

Gerson.

P.S.: This is a closer quotation of the Scientific American Article:

Quote:WALLIS PRODUCT (1665)

Code:

π   2×2   4×4   6×6   8×8
– = ––– × ––– × ––– × ––– × …
2   1×3   3×5   5×7   7×9

GREGORY’S SERIES (1671)

Code:


π       1   1   1
– = 1 - – + – - – + …
4       3   5   7

TERMS OF MATHEMATICAL SEQUENCES can be summed or multiplied to yield values for pi (divided by a constant) or its reciprocal. The first two sequences, discovered respectively by the mathematicians John Wallis and James Gregory, are probably among the best-known, but they are practically useless for computational purposes. Not even 100 years of computing on a supercomputer programmed to add or multiply the terms of either sequence would yield 100 digits of pi. The formula discovered by John Machin made the calculation of pi feasible, since calculus allows the inverse tangent (arc tangent) of a number, x, to be expressed in terms of a sequence whose sum converges more rapidly to the value of the arc tangent the smaller x is.

[Ángel Martin]

Hi Gerson (and All),

Obviously pi continues to be a top-subject on the forum, with fascinating contributions being added every few months - and that happening for years already. Your latest ones have made me decide to put together a PI/E ROM as an attempt to gather in a common place (i.e. an HP-41 Module) some of the material the forum members have been producing. Since my math prowess is limited and won't allow for more substantial contributions here, this ROM is my tribute to the forum activity on the pi subject.

Obviously the selection is biased by my own preferences, but hopefully I didn't miss the more original or intriguing ones. Besides your latest articles, the sources used obviously include Valentín's wondrous challenges and frequent articles, a treasure trove on the subject. The hyperlinks and forum archives have also been very helpful to compile the material, of course.

To add some value to the project I've made MCODE versions of some algorithms, making them faster and slightly more accurate (due to the 13-digit routines) - though this is a negligible feature in the transcendental pi-world...

Anyway, be as limited as it may, the ROM should be ready in a few more days - ideally on the 14th. but no assurances are made ;-)

Thanks for the contributions, a lot of fun came from working on it.
Best,
ÁM

[EdS2]

(03-08-2022 05:17 PM)Gerson W. Barbosa Wrote:  Please correct it if you haven’t done so yet. This should make the program a little less slow.

Thanks! There is indeed a speedup. I'm now going through to change to integer variables to see what difference that can make.

Quote:Also, I can’t make your links work.

Sorry about that - I will publish some updated ones, perhaps on the 14th!

(03-09-2022 08:12 AM)Ángel Martin Wrote:  Hi Gerson (and All),

Obviously pi continues to be a top-subject on the forum, with fascinating contributions being added every few months - and that happening for years already. Your latest ones have made me decide to put together a PI/E ROM...

Excellent news! (I hope you find a spigot approach to include: such a thing will produce digits successively from left to right.)

(03-08-2022 05:42 PM)Valentin Albillo Wrote:  ... there's this:

Speeding it up !

Thanks for the link, Valentin, and again to Gerson for all the ideas and links.

[Gerson W. Barbosa]

(03-10-2022 07:47 AM)EdS2 Wrote:  

(03-08-2022 05:17 PM)Gerson W. Barbosa Wrote:  Please correct it if you haven’t done so yet. This should make the program a little less slow.

Thanks! There is indeed a speedup. I'm now going through to change to integer variables to see what difference that can make.

Thank you! I’ve updated the listing (will update the timings later).
Times on the HP-75C:
100 digits: 8083.154 s ( ~ 135 min )
152 digits: 24386.212 s ( ~ 406 min )

[Valentin Albillo]

.
Hi, EdS2,

EdS2 Wrote:Excellent news! (I hope you find a spigot approach to include: such a thing will produce digits successively from left to right.)
[...]
Thanks for the link, Valentin, [...]

You're welcome. I know you like spigot algorithms for \(\pi\), for instance you ported my multiprecision BASIC spigot program for a SHARP pocket computer to the BBC micro, which ran flawlesly (and faster when using integer variables throughout), while I also ported it to the HP-71B, which could produce many thousands of digits (one at a time, mind you) if you had enough RAM and time.

Curiously enough, using INTEGER variables on the 71B results in a slower program than using REAL variables.

As for what Ángel might include in his forthcoming PI/E ROM, who knows ...

Best regards.
V.

[EdS2]

Indeed I do like a pi program, and a spigot program, and anything multiprecision.

About porting to BBC Basic - in the case of your spigot, a very naive initial port would produce 50 digits in 118 seconds. With a change to integer variables, it becomes 66 seconds, and then by changing / to DIV, it becomes 54 seconds. We can even save another second by setting D% to 10 and using that instead of the literal in a couple of places.

So by recoding in the fastest idiom for BBC Basic, we get a 2x speedup. (The above timings are for Basic 2 on a 2MHz 6502 - later Basics were faster, and of course CPU clock speeds also improved.)

Other changes which can be beneficial:
- using only single-letter variables, especially A% to Z% which are special-cased
- removing the loop variable from NEXT
- using FOR, REPEAT, PROC and FN instead of line numbers
- inlining functions and subroutines
- removing blank spaces and consolidating multiple lines

Of course, some of these re-codings make the program a bit less clear, and one gets diminishing returns. It's a pity that all the % signs for integer variables are a bit ugly especially if you're not used to them: BBC Basic has no DEF INT.

In the case of Gerson's program, the version below produces 12 digits in 14 seconds, compared to 43 seconds for the initial port. This version produces 36 digits in 186s. (Edit: and 100 digits in 2510s.)

I was quite pleased to do without both FN definitions and especially the 5E-6 which I didn't understand or feel comfortable with. I could use DIV and MOD in a natural way.

Code:


10 REM  ***  PI BY WALLIS-WASICKI FORMULA  ***
20 REM CLS: REM KEYOFF: WIDTH 40
60 DIM A%(50), B%(50), C%(100), Q%(50), X%(50), Y%(100), W%(50)
70 INPUT "Number of decimal digits: ";N$: ND=VAL(N$): PRINT
80 B% = 10000: N%=ND DIV 4 + 2: K%=ND*12 DIV 25 + 1
85 M%=10*LOG(ND)/LOG(1000)
90 TIME=0
100 D%=8*K%+4: E%=D%-8: G%=4*K%*K%-1
110 W%(1)=4: B%(1)=0
120 FOR I%=2 TO N%
130   W%(I%)=0: B%(I%)=0
140 NEXT
150 FOR IT=1 TO K%
160   FOR I%=0 TO N%
170     Q%(I%)=W%(I%)
180   NEXT
190   GOSUB 860
200   GOSUB 930
210   B%(1)=B%(1)+D%
220   GOSUB 420
230   GOSUB 1000
240   FOR I%=1 TO N%: B%(I%)=C%(I%): NEXT
250   G%=G%-E%
260   E%=E%-8
270 NEXT
280 B%(1)=B%(1)+D%-1
290 GOSUB 420
300 FOR I%=1 TO N%
310   A%(I%)=W%(I%): X%(I%)=C%(I%)
320 NEXT
330 X%(2)=X%(2)+B%/2
340 GOSUB 670
350 TM = TIME/100
360 FOR I%=1 TO N%
370   @%=8:PRINT ;" ";INT(.5+C%(I%));
380 NEXT
390 PRINT: PRINT: PRINT "Time: ";TM;: PRINT "s" : PRINT
400 END
410 REM  ***   LINV   ***
420 T=1/(B%(1)+B%(2)/B%)
430 X%(1)=T: X%(2)=(T*B%)MODB%
450 FOR L%=1 TO M%
460   FOR I%=1 TO N%
470     A%(I%)=B%(I%)
480   NEXT
490   GOSUB 670
500   R%=0
510   FOR I%=N% TO 2 STEP -1
520     T%=(B%-R%-C%(I%))MOD B%
530     R%=R% OR SGN(C%(I%))
540     Y%(I%)=T%
550   NEXT
560   Y%(1)=2-R%-C%(1)
570   FOR I%=1 TO N%
580     A%(I%)=Y%(I%)
590   NEXT
600   GOSUB 670
610   FOR I%=1 TO N%
620     X%(I%)=C%(I%)
630   NEXT
640 NEXT
650 RETURN
660 REM  ***   LMULT   ***
670 FOR I%=1 TO N%
680   C%(I%)=0
690 NEXT
700 FOR I%=N% TO 1 STEP -1
710   P%=0
720   FOR J%=N%-I%+2 TO -1 STEP -1
730     T%=(P%+A%(I%)*X%(J%+1))
750     C%(I%+J%)=C%(I%+J%)+T%MODB%
760     P%=T%DIVB%
770   NEXT
780 NEXT
790 FOR I%=N%+2 TO 2 STEP -1
800   T%=C%(I%)
810   C%(I%)=T%MODB%
820   C%(I%-1)=C%(I%-1)+T% DIV B%
830 NEXT
840 RETURN
850 REM  ***   DIVIDE   ***
860 Y%=0
870 FOR I%= 1 TO N%
880   X%=W%(I%)+Y%*B%
890   W%(I%)=X%DIVG%: Y%=X% MODG%
900 NEXT
910 RETURN
920 REM  ***   ADD   ***
930 Y%=0
940 FOR I%=N% TO 1 STEP -1
950   X%=W%(I%)+Q%(I%)+Y%
960   W%(I%)=X%MODB%: Y%=X%DIVB%
970 NEXT
980 RETURN
990 REM  ***   MULT   ***
1000 Y%=0
1010 FOR I%=N% TO 1 STEP -1
1020   X%=C%(I%)*G%+Y%
1030   C%(I%)=X%MODB%: Y%=X%DIVB%
1040 NEXT
1050 RETURN