Gerson's Pi Program

[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