'Subject: Calculation of pi with the 1976 Salamin-Brent
'         arithmetic-geometric mean (AGM) algorithm.
'Author : Sjoerd.J.Schaper
'Date   : 11-11-2004
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST A = 0, B = 1, D = 2, S = t3
i = LargeInit(D, "pi_agm")
DIM e(13): e(1) = 2: e(2) = 4: e(3) = 10
e(4) = 21: e(5) = 43: e(6) = 86: e(7) = 174
e(8) = 348: e(9) = 698: e(10) = 1394
e(11) = 2793: e(12) = 5587: e(13) = 11176
CLS ' number of correct digits

begin:
LOCATE 1, 1
PRINT "   +--------------------------+"
PRINT "   |  AGM-iterations [1..12], |"
PRINT "   |  0 to exit.              |"
PRINT "   |                          |"
PRINT "   +--------------------------+"
LOCATE 4, 7: INPUT ; "imx"; g$
IF INSTR(g$, "'") THEN GOTO begin
IF VAL(g$) < 1 OR VAL(g$) > 13 THEN
  PRINT "value out of range": GOTO eind
END IF
LOCATE 7, 2
imx = VAL(g$): tim# = TIMER

Letf A, 10: Powr A, e(imx) + 4 '      initialize arithmetic mean A,
Squ A, S: Rsft S, 1 '                 sum S[0]:= A^2 / 2,
Isqrt S, B '                          geometric mean B
FOR i = 1 TO imx
   Dup A, D: PRINT i;
   Subt D, B: Rsft D, 1 '             D:= (A[n-1] - B[n-1]) / 2
   Squ D, t0: Lsft t0, i
   Subt S, t0 '                       S[n]:= S[n-1] - D^2 * 2^i
   Mult A, B, D
   Add A, B: Rsft A, 1 '              A[n]:= (A[n-1] + B[n-1]) / 2
   Isqrt D, B '                       B[n]:= Sqrt(A[n-1] * B[n-1])
NEXT i
PRINT
Squ A, D: Lsft D, 1 '                 D:= 2 * A[n]^2
Printr D, S, "pi = ", e(imx)
Prints STR$(e(imx)) + " digits", 2

eind:
g$ = STR$(CSNG(TIMER - tim#))
Prints "Timer:" + g$ + " s", 1
CALL Term
SYSTEM

