Gerson's Pi Program

[Gerson W. Barbosa]

(03-11-2022 12:42 PM)EdS2 Wrote:  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.

Very good work! Thank you.
5E-6 was added to avoid errors due to binary representation. For instance, I notice a number that looked like 4482 on the BBC Micro, but its integer part was 4481. When I took its fractional part the result was 0.9999996185. So, INT(4481 + 0.999999 + 0.000005) will return 4482 as it should. Of course, this can be fixed by working with integer variables and functions whenever possible. One hundred digits in 42 minutes is quite good, considering I estimate it will take twice as long on the MSX computer (Z80 @ 3.57 Mhz). Perhaps I should improve mine also.

The multiprecision program is basically equivalent to the obfuscated one in the card, except that the expression for π is 4*W*(1/(C+D-1)+1/2) instead of W*(2/(C+D-1)+1). Also, long divisions are performed by multiplying by the inverse (for example M/(C + D) = M*1/(C + D).
I wrote the LMULT routine (multiply two n-digit numbers giving an n-digit result) many years ago. Usually if I can do something by hand I don't have trouble programming it. But when I tried to do a long division by hand, I discovered I don't remember how to do it anymore. As LDIVIDE was needed, I did 1/y using Newton's method: xₙ = xₙ₋₁(2 - xₙ₋₁.y). The initial guess, x₀, is given with the first few exact digits to ensure a faster convergence (lines 440 and 450 in the HP-75C program below). ADD, MULT and DIVID were taken from Katie Wasserman's program.

HP-75C program:

Code:


10 REM  **   PI BY WALLIS-WASICKI FORMULA   ***           
20 ASSIGN # 1 TO 'PI'
30 OPTION BASE 0
40 INTEGER B0,D,E,G,I,I0,J,K,L,M,N,N0
50 DIM A(55),B(55),C(110),Q(55),X(55),Y(110),W(55)
60 INPUT "Number of decimal digits: "; N$ @ N0=VAL(N$)
70 B=1000000 @ N=N0 DIV 6+2 @ K=N0*12 DIV 25+1
80 T0=TIME
90 D=8*K+4 @ E=D-8 @ G=4*K*K-1
100 FOR I=0 TO N+2
110 W(I)=0 @ A(I)=0 @ B(I)=0 @ C(I)=0 @ X(I)=0
120 NEXT I
130 W(1)=4
140 FOR I0=1 TO K
150 FOR I=0 TO N
160 Q(I)=W(I)
170 NEXT I
180 GOSUB 880
190 GOSUB 950
200 B(1)=B(1)+D
210 GOSUB 440
220 GOSUB 1020
230 FOR I=1 TO N
240 B(I)=C(I)
250 NEXT I
260 G=G-E
270 E=E-8
280 NEXT I0
290 B(1)=B(1)+D-1
300 GOSUB 440
310 FOR I=1 TO N     
320 A(I)=W(I) @ X(I)=C(I)
330 NEXT I
340 X(2)=X(2)+B/2
350 GOSUB 690
360 T0=TIME-T0 @ T$=STR$(T0) @ T$=T$&' s'
370 FOR I=1 TO N
380 PRINT # 1 ; C(I)
390 NEXT I
400 PRINT # 1 ; T$
410 ASSIGN # 1 TO *
420 END
430 REM   ***   LINV   ***
440 T=1/(B(1)+B(2)/B)
450 X(1)=IP(T) @ X(2)=IP(B*FP(T))
460 M=10*LOG(N0)/LOG(1000)-1
470 FOR L=1 TO M               
480 FOR I=1 TO N
490 A(I)=B(I)
500 NEXT I
510 GOSUB 690
520 B0=0
530 FOR I=N TO 2 STEP -1
540 T=MOD(B-B0-C(I),B)
550 B0=B0 OR SGN(C(I))
560 Y(I)=T
570 NEXT I
580 Y(1)=2-B0-C(1)
590 FOR I=1 TO N
600 A(I)=Y(I)
610 NEXT I
620 GOSUB 690       
630 FOR I=1 TO N
640 X(I)=C(I)
650 NEXT I
660 NEXT L
670 RETURN
680 REM  ***   LMULT   ***
690 FOR I=1 TO N
700 C(I)=0
710 NEXT I
720 FOR I=N TO 1 STEP -1
730 T1=0 @ A0=0
740 FOR J=N-I+2 TO -1 STEP -1
750 C0=(T1+A(I)*X(J+1))/B
760 A0=IP(B*FP(C0))
770 C(I+J)=C(I+J)+A0
780 T1=IP(C0)
790 NEXT J
800 NEXT I
810 FOR I=N+2 TO 2 STEP -1
820 T=C(I)/B
830 C(I)=IP(B*FP(T))
840 C(I-1)=C(I-1)+IP(T)
850 NEXT I
860 RETURN
870 REM  ***   DIVIDE   ***
880 Y=0
890 FOR I=1 TO N
900 X=W(I)+Y*B
910 W(I)=IP(X/G) @ Y=MOD(X,G)
920 NEXT I
930 RETURN
940 REM  ***   ADD   ***       
950 Y=0
960 FOR I=N TO 1 STEP -1
970 X=W(I)+Q(I)+Y
980 W(I)=MOD(X,B) @ Y=IP(X/B)
990 NEXT I
1000 RETURN
1010 REM  ***   MULT   ***
1020 Y=0
1030 FOR I=N TO 1 STEP -1
1040 X=C(I)*G+Y
1050 C(I)=MOD(X,B) @ Y=IP(X/B)
1060 NEXT I
1070 RETURN

The digits of π are written to a BASIC file with DATA lines suitable for occasional further processing (like insertion of the missing zeros at left, for example).

1 DATA 3
2 DATA 141592
3 DATA 653589
4 DATA 793238
5 DATA 462643
6 DATA 383279
7 DATA 502884
8 DATA 197169
9 DATA 399375
10 DATA 105820
11 DATA 974944
12 DATA 592307
13 DATA 816406
14 DATA 286208
15 DATA 998628
16 DATA 34825
17 DATA 342117
18 DATA 67967
19 DATA '8083.154 s'

100 correct digits in 135 minutes, but a really competent HP-75C programmer can surely improve on that.

P.S.: Here's a QBASIC version based on you port. The number of digits is forced to the next multiple of 4. All digits are supposed to be exact, no rounding in the last one. Because of the choice for 32-bit integers, the number of digits is limited to 480. One hundred digits in 71 seconds (DOSBox with default settings). 480 digits in about 27 minutes (when compiled). I would like to try it on my HP-200LX, but its screen has gone dark.

Code:


10 REM  ***  PI BY WALLIS-WASICKI FORMULA  ***
30 DEFLNG A-Y
40 DEFSNG Z
60 DIM A(125), B(125), C(250), Q(125), X(125), Y(250), W(125)
70 CLS : LINE INPUT "Number of decimal digits: "; N$: ND = VAL(N$): IF ND > 480 THEN 70
75 ND = 4 * (ND \ 4 + SGN(ND MOD 4)): LOCATE 1, 26: PRINT ND: PRINT
80 B = 10000: N = ND \ 4 + 2: K = ND * 12 \ 25 + 1
85 M = 10 * LOG(ND) / LOG(1000)
90 TM = TIMER
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 = TIMER - TM
355 PRINT " "; STR$(C(1)); ".";
357 V = 0
360 FOR I = 2 TO N - 1
361   V = V + 1
362   IF V MOD 16 = 0 THEN V = 1: PRINT : PRINT "    ";
364   C$ = STR$(C(I)):  C$ = RIGHT$(C$, LEN(C$) - 1)
366   WHILE 4 - LEN(C$) <> 0
367     C$ = "0" + C$
369   WEND
370   PRINT C$; " ";
380 NEXT
390 PRINT : PRINT : PRINT "Time: "; TM; : PRINT "s": PRINT
400 END
410 REM  ***   LINV   ***
420 Z = 1 / (B(1) + B(2) / B)
430 X(1) = Z: X(2) = (Z * B) MOD B
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 MOD B
760     P = T \ B
770   NEXT
780 NEXT
790 FOR I = N + 2 TO 2 STEP -1
800   T = C(I)
810   C(I) = T MOD B
820   C(I - 1) = C(I - 1) + T \ 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 \ G: Y = X MOD G
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 MOD B: Y = X \ B
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 MOD B: Y = X \ B
1040 NEXT
1050 RETURN

[Gerson W. Barbosa]

(03-11-2022 12:42 PM)EdS2 Wrote:  

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

On the MSX LOG is the natural logarithm whereas on the BBC Micro it is the decimal logarithm. This won't change the result for M% in line 85 above, but it can be simplified to

85 M%=10/3*LOG(ND)

Also, matrices are initialized to zero on both the MSX and the BBC micro. Thus, we can eliminate the lines 120 through 140 and edit the line 110 to

110 W%(1)=4

My mistake, especially because more matrices must be initiliazed to zero as you can see from the HP-75C listing above.

It's already Pi Day in countries which use MM/DD date format, so I decided to turn on my old MSX computer after 15 years to check the running time for 100 digits. The 2-byte system variable TIME allows for timings up to 1092 seconds only, but this limitation can be overcome by using an interruption to updated a timing variable every 10 seconds.



Same results on the BlueMSX emulator, in one tenth of the time.

Code:


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

The matrices are dimensioned for 300 digits or so, but G will overflow for ND around 180, unless the line 60 is changed to

60 DEFINT D,E,H-N

In this case, the constant in line 10 below should be changed to 1000.

Code:


10 CLEAR 500: MAXFILES=1: OPEN "PI" FOR INPUT AS #1: INPUT N%: INPUT #1,A$: A$=A$+".":P$=A$
20 IF EOF(1) THEN 60
30 INPUT #1,A$: L%=LEN(A$):Z$=""
40 IF L%<7 THEN FOR I%=1 TO 7-L%: Z$=Z$+"0":NEXT: A$=Z$+A$
50  P$=P$+A$: GOTO 20
60 FOR I%=1 TO N%+2: PRINT MID$(P$,I%,1);: NEXT: CLOSE #1

———

PS. Sure 86 minutes for just the first 100 digits is no big deal, but considering the most powerful supercomputer of the late eighties wouldn’t do it in even 100 years using the plain Wallis Product, this looks great.
The HP-49G calculator, released in 1999, does it in 137.3 seconds. In the picture, the slightly more recent HP 50g computes 100 decimal digits of \(\pi\) in one third of that time using an equivalent User-RPL program:



Edited to fix a couple of elusive typos.

[EdS2]

Very nice to run it on a real machine!