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.
Search
Member List
Calendar
Help
