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