/* --------------------------------------------------------------------
! 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  = -7.792440345582e-002
!
!  Estimated error =  1.041851518710e-017
!
!  Error code = 0.000000
!
!  Number of function evaluations = 33
!
! ---------------------------------------------------------------------
! Reference: From Numath Library By Tuan Dang Trong in Fortran 77
!            [BIBLI 18].
!
!                               C++ Release 1.0 By J-P Moreau, Paris
!                                          (www.jpmoreau.fr)
! ------------------------------------------------------------------- */
#include <stdio.h>
#include <math.h>

    double AERROR,CODE,ERROR,RERROR,X1,X2,VAL;
    int NBRF;

    double FCT(double X) {
      return (cos(X)-2*sin(X));
    }

    double MAX(double a,double b) {
      if (a>=b) return a;
      else return b;
    }

// ---------------------------------------------------------------------
    void QANC8 (double A,double B,double AERR,double RERR,
		double *RES, double *ERR, int *NBF, double *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
!     ERR     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.
! --------------------------------------------------------------------- */
//  Labels: e30, e50, e60, e62, e70, e72, e75, e80, e82;
      double QR[32];
      double F[17], X[17];
      double FS[9][31], XS[9][31];
      double COR,F0,PAS,PAS1,QP,SUM,W0,W1,W2,W3,W4,X0;
      double ERR1,QD,QL,QN,TEMP,TOL1;
      int I,J,L,LMIN,LMAX,LOUT,NFIN,NIM,NMAX;

      LMIN = 1;
      LMAX = 30;
      LOUT = 6;
      NMAX = 5000;
      NFIN = (int) floor(NMAX-8*(LMAX-LOUT+pow(2,LOUT+1)));
      W0   =   3956.0/14175.0;
      W1   =  23552.0/14175.0;
      W2   =  -3712.0/14175.0;
      W3   =  41984.0/14175.0;
      W4   = -18160.0/14175.0;
      FLG  = 0;
      *RES  = 0;
      COR  = 0;
      *ERR  = 0;
      SUM  = 0;
      *NBF  = 0;
      if (A==B) return;
      L = 0;
      NIM = 1;
      X0  = A;
      X[16] = (double) B;
      QP  = 0;
      F0   = FCT(X0);
      PAS1 = (B-A)/16.0;
      X[8] = (X0+X[16])*0.5;
      X[4] = (X0+X[8])*0.5;
      X[12]= (X[8]+X[16])*0.5;
      X[2] = (X0+X[4])*0.5;
      X[6] = (X[4]+X[8])*0.5;
      X[10]= (X[8]+X[12])*0.5;
      X[14]= (X[12]+X[16])*0.5;
      J=2;
      while (J<=16) {
        F[J] = FCT(X[J]);
        J+=2;
      }
      *NBF = 9;
e30:  X[1] = (X0+X[2])*0.5;
      F[1] = FCT(X[1]);
      J=3;
      while (J<=15) {
        X[J] = (X[J-1]+X[J+1])*0.5;
        F[J] = FCT(X[J]);
        J+=2;
      }
      *NBF = *NBF+8;
      PAS = (X[16]-X0)/16.0;
      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 = fabs(QD)/1023.0;
      TOL1 = MAX(AERR,RERR*fabs(SUM))*(PAS/PAS1);
      if (L < LMIN) goto e50;
      if (L >= LMAX) goto e62;
      if (*NBF > NFIN) goto e60;
      if (ERR1 <= TOL1) goto e70;
e50:  NIM *= 2;
      L++;
      for (I = 1; I<=8; I++) {
        FS[I][L] = F[I+8];
        XS[I][L] = X[I+8];
      }
      QP = QL;
      for (I = 1; I<=8; I++) {
        F[18-2*I] = F[9-I];
        X[18-2*I] = X[9-I];
      }
      goto e30;

e60:  NFIN *= 2;
      LMAX = LOUT;
      *FLG = *FLG + (B-X0)/(B-A);
      goto e70;
e62:  *FLG = *FLG + 1.0;
e70:  *RES = *RES + QN;
      *ERR = *ERR + ERR1;
      COR = COR + QD/1023.0;
e72:  if (NIM==2*(NIM/2)) goto e75;
      NIM /= 2;
      L--;
      goto e72;
e75:  NIM++;
      if (L <= 0) goto e80;
      QP = QR[L];
      X0 = X[16];
      F0 = F[16];
      for (I = 1; I<=8; I++) {
        F[2*I] = FS[I][L];
        X[2*I] = XS[I][L];
      }
      goto e30;

e80:  *RES = *RES + COR;
      if (*ERR==0) return;
e82:  TEMP = fabs(*RES) + *ERR;
      if (TEMP!=fabs(*RES)) return;
      *ERR *= 2.0;
      goto e82;
}


void main()  {
  
  X1=0.0;
  X2=1.0;
  AERROR=1e-9;
  RERROR=1e-10;

  QANC8(X1,X2,AERROR,RERROR,&VAL,&ERROR,&NBRF,&CODE);

  printf("\n Integral Value  = %20.12e\n\n", VAL);
  printf(" Estimated error = %20.12e\n\n", ERROR);
  printf(" Error Code = %f\n\n", CODE);
  printf(" Number of function evaluations = %d\n\n", NBRF);
     
}

// end of file tqanc8.cpp