Happy Pi Day!

[Gerson W. Barbosa]

The MSX BASIC code runs with little modification in QBASIC, mainly because the former was also developed by Microsoft. Thus, by running it in non-emulated DOS we can achieve much greater speed even in my old desktop computer from fourteen years ago: 1002 decimal places in "only" 248 seconds (compiled) or 307 seconds (interpreted). BTW, 276 decimal places in 6 seconds (2 seconds when compiled). The program works also on Turbo BASIC, albeit about 30% slower.

Code:


10 REM PI BY WALLIS-WASICKI FORMULA - QBASIC version
20 CLS
30 MAXFILES = 1: OPEN "PI.DAT" FOR OUTPUT AS #1
60 DEFINT D-E, H-N, V
65 DEFDBL A-C, F-G, O-U, W-Z
67 DEF FNFRAC (X) = X - INT(X)
68 DEF FNRMDR (X, Y) = X - Y * INT(X / Y)
70 DIM A(171), B(171), C(342), Q(171), X(171), Y(342), W(171)
80 LINE INPUT "Number of decimal digits: "; N$: ND = VAL(N$): PRINT
90 ND = 6 * (ND \ 6 + SGN(ND MOD 6)): LOCATE 1, 26: PRINT ND: PRINT
100 B = 1000000: N = ND \ 6 + 2: K = ND * 12 \ 25 + 1
110 M = 10 * LOG(ND) / LOG(1000) - 1
130 TM = TIMER
140 D = 8 * K + 4: E = D - 8: G = 4# * K * K - 1
150 W(1) = 4
160 FOR IT = 1 TO K
170   FOR I = 0 TO N
180     Q(I) = W(I)
190   NEXT I
200   GOSUB 990
210   GOSUB 1060
220   B(1) = B(1) + D
230   GOSUB 560
240   GOSUB 1130
250   FOR I = 1 TO N
260     B(I) = C(I)
270   NEXT I
280   G = G - E
290   E = E - 8
300 NEXT IT
310 B(1) = B(1) + D - 1
320 GOSUB 560
330 FOR I = 1 TO N
340   A(I) = W(I): X(I) = C(I)
350 NEXT I
360 X(2) = X(2) + B / 2
370 GOSUB 800
375 TM = TIMER - TM
380 C$ = STR$(C(1)): PRINT TAB(6); : PRINT C$; "."; : PRINT #1, C$
390 V = 0
400 FOR I = 2 TO N - 1
410   V = V + 1
415   IF V MOD 11 = 0 THEN V = 1: PRINT : PRINT "        ";
420   C$ = STR$(C(I)): C$ = RIGHT$(C$, LEN(C$) - 1)
430   IF 6 - LEN(C$) = 0 THEN 450
440     C$ = "0" + C$
445   GOTO 430
450   PRINT C$; " ";
460   PRINT #1, C$
470 NEXT I
480 PRINT #1, STR$(C(N))
490 PRINT : PRINT : PRINT "Runtime:"; : PRINT INT(TM + .5); : PRINT "s": PRINT
500 CLOSE #1
510 END
550 REM  ***   LINV   ***
560 T = 1 / (B(1) + B(2) / B)
570 X(1) = INT(T): X(2) = INT(B * FNFRAC(T) + .01)
580 FOR L = 1 TO M
590   FOR I = 1 TO N
600     A(I) = B(I)
610   NEXT I
620   GOSUB 800
630   BR% = 0
640   FOR I = N TO 2 STEP -1
650     T = FNRMDR(B - BR% - C(I), B)
660     BR% = BR% OR SGN(C(I))
670     Y(I) = T
680   NEXT I
690   Y(1) = 2 - BR% - C(1)
700   FOR I = 1 TO N
710     A(I) = Y(I)
720   NEXT I
730   GOSUB 800
740   FOR I = 1 TO N
750     X(I) = C(I)
760   NEXT I
770 NEXT L
780 RETURN
790 REM  ***   LMULT   ***
800 FOR I = 1 TO N
810   C(I) = 0
820 NEXT I
830 FOR I = N TO 1 STEP -1
840   TP = 0: AC = 0
850   FOR J = N - I + 2 TO -1 STEP -1
860     CT = (TP + A(I) * X(J + 1)) / B
870     AC = INT(B * FNFRAC(CT) + .01)
880     C(I + J) = C(I + J) + AC
890     TP = INT(CT)
900   NEXT J
910 NEXT I
920 FOR I = N + 2 TO 2 STEP -1
930   T = C(I) / B
940   C(I) = INT(B * FNFRAC(T) + .01)
950   C(I - 1) = C(I - 1) + INT(T)
960 NEXT I
970 RETURN
980 REM  ***   DIVIDE   ***
990 Y = 0
1000 FOR I = 1 TO N
1010   X = W(I) + Y * B
1020   W(I) = INT(X / G): Y = FNRMDR(X, G)
1030 NEXT I
1040 RETURN
1050 REM  ***   ADD   ***
1060 Y = 0
1070 FOR I = N TO 1 STEP -1
1080   X = W(I) + Q(I) + Y
1090   W(I) = FNRMDR(X, B): Y = INT(X / B)
1100 NEXT I
1110 RETURN
1120 REM  ***   MULT   ***
1130 Y = 0
1140 FOR I = N TO 1 STEP -1
1150   X = C(I) * G + Y
1160   C(I) = FNRMDR(X, B): Y = INT(X / B)
1170 NEXT I
1180 RETURN