' ---------------------------------------------------------------------
' Program to integrate a user-defined function f(x) from x1 to x2 by
' the QANC8 subroutine with control of absolute and relative precisions
' ---------------------------------------------------------------------
' SAMPLE RUN:
'
'  (Integrate function cos(x) - 2 sin(x) from x=0 to x=1,
'   with a precision <= 1e-10).
'
'  Integral Value  = -.7792440345582405
'
'  Estimated error =  9.72394750795836D-17
'
'  Error code = 0
'
'  Number of function evaluations = 33
'
' ---------------------------------------------------------------------
' Reference: From Numath Library By Tuan Dang Trong in Fortran 77
' [BIBLI 18].  
'
'                              Basic Release 1.0 By J-P Moreau, Paris
'                                        (www.jpmoreau.fr)
' ---------------------------------------------------------------------
'PROGRAM TQANC8

DEFDBL A-H, O-Z
DEFINT I-N

  A = 0#
  B = 1#
  AERR = .000000001#
  RERR = .0000000001#

  GOSUB 1000 'CALL QANC8(A,B,AERR,RERR,RES,XERR,NBF,FLG)

  CLS
  PRINT
  PRINT " Integral Value  = "; RES
  PRINT
  PRINT " Estimated error = "; XERR
  PRINT
  PRINT " Error Code ="; FLG
  PRINT
  PRINT " Number of function evaluations ="; NBF
  PRINT
     
END

500 'FUNCTION FCT(X)
  FCT = COS(X) - 2# * SIN(X)
RETURN


' ----------------------------------------------------------------------
1000  'SUBROUTINE QANC8 (A,B,AERR,RERR,RES,XERR,NBF,FLG)

'     INTEGRATE A REAL FUNCTION FCT(X) FROM X=A TO X=B, WITH
'     GIVEN ABSOLUTE AND RELATIVE PRECISIONS, AERR, RERR.

'     INPUTS:
'     FCT     EXTERNAL USER-DEFINED FUNCTION FOR ANY X VALUE
'             IN INTERVAL (A,B)
'     A,B     LIMITS OF INTERVAL
'     AERR,RERR   RESPECTIVELY ABSOLUTE ERROR AND RELATIVE ERROR
'                 REQUIRED BY USER

'     OUTPUTS:
'     RES     VALUE OF INTEGRAL
'     XERR    ESTIMATED ERROR
'     NBF     NUMBER OF NECESSARY FCT(X) EVALUATIONS
'     FLG     INDICATOR
'             = 0.0       CORRECT RESULT
'             = NNN.RRR   NO CONVERGENCE DU TO A SINGULARITY.
'             THE SINGULAR POINT ABCISSA IS GIVEN BY FORMULA:
'             XS = B-.RRR*(B-A)

'     REFERENCE :
'     FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL
'     COMPUTATIONS. PRENTICE-HALL, INC.
' -----------------------------------------------------------------------
      DIM QR(31), F(16), X(16), FS(8, 30), XS(8, 30)
      LMIN = 1
      LMAX = 30
      LOUT = 6
      NMAX = 5000
      NFIN = NMAX - 8 * (LMAX - LOUT + 2 ^ (LOUT + 1))
      W0 = 39560 / 14175#
      W1 = 235520 / 14175#
      W2 = -37120 / 14175#
      W3 = 419840 / 14175#
      W4 = -181600 / 14175#
      FLG = 0#
      RES = 0#
      COR = 0#
      XERR = 0#
      SUM = 0#
      NBF = 0
      IF A = B THEN RETURN
      L = 0
      NIM = 1
      X0 = A
      X(16) = B
      QP = 0#
      X = X0: GOSUB 500: F0 = FCT
      PAS1 = (B - A) / 16#
      X(8) = (X0 + X(16)) * .5#
      X(4) = (X0 + X(8)) * .5
      X(12) = (X(8) + X(16)) * .5#
      X(2) = (X0 + X(4)) * .5#
      X(6) = (X(4) + X(8)) * .5#
      X(10) = (X(8) + X(12)) * .5#
      X(14) = (X(12) + X(16)) * .5#
      FOR J = 2 TO 16 STEP 2
	X = X(J): GOSUB 500: F(J) = FCT
      NEXT J
      NBF = 9
30    X(1) = (X0 + X(2)) * .5#
      X = X(1): GOSUB 500: F(1) = FCT
      FOR J = 3 TO 15 STEP 2
	X(J) = (X(J - 1) + X(J + 1)) * .5#
	X = X(J): GOSUB 500: F(J) = FCT
      NEXT J
      NBF = NBF + 8
      PAS = (X(16) - X0) / 16#
      QL = (W0 * (F0 + F(8)) + W1 * (F(1) + F(7)) + W2 * (F(2) + F(6)) + W3 * (F(3) + F(5)) + W4 * F(4)) * PAS
      QR(L + 1) = (W0 * (F(8) + F(16)) + W1 * (F(9) + F(15)) + W2 * (F(10) + F(14)) + W3 * (F(11) + F(13)) + W4 * F(12)) * PAS
      QN = QL + QR(L + 1)
      QD = QN - QP
      SUM = SUM + QD
      ERR1 = ABS(QD) / 1023#
      IF AERR > RERR * ABS(SUM) THEN
	TEMP1 = AERR
      ELSE
	TEMP1 = RERR * ABS(SUM)
      END IF
      TOL1 = TEMP1 * (PAS / PAS1)
      IF L < LMIN THEN GOTO 50
      IF L >= LMAX THEN GOTO 62
      IF NBF > NFIN THEN GOTO 60
      IF ERR1 <= TOL1 THEN GOTO 70
50    NIM = 2 * NIM
      L = L + 1
      FOR I = 1 TO 8
	FS(I, L) = F(I + 8)
	XS(I, L) = X(I + 8)
      NEXT I
      QP = QL
      FOR I = 1 TO 8
	F(18 - 2 * I) = F(9 - I)
	X(18 - 2 * I) = X(9 - I)
      NEXT I
      GOTO 30
60    NFIN = 2 * NFIN
      LMAX = LOUT
      FLG = FLG + (B - X0) / (B - A)
      GOTO 70
62    FLG = FLG + 1!
70    RES = RES + QN
      XERR = XERR + ERR1
      COR = COR + QD / 1023!
72    IF NIM / 2 = INT(NIM / 2) THEN GOTO 75
      NIM = NIM / 2
      L = L - 1
      GOTO 72
75    NIM = NIM + 1
      IF L <= 0 THEN GOTO 80
      QP = QR(L)
      X0 = X(16)
      F0 = F(16)
      FOR I = 1 TO 8
	F(2 * I) = FS(I, L)
	X(2 * I) = XS(I, L)
      NEXT I
      GOTO 30
80    RES = RES + COR
      IF ABS(XERR) = 0# THEN RETURN
82    TEMP = ABS(RES) + XERR
      IF TEMP <> ABS(RES) THEN RETURN
      XERR = 2# * XERR
      GOTO 82
RETURN

' end of file tqanc8.bas