C      ALGORITHM 699 , COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 17, NO. 4, DECEMBER, 1991, PP. 457-461.
      program TEQUAD
c
c     Test the re-implementation of Patterson's QUAD.
c
c     Parameters
c
c NEPS    is the number of tolerances at which to exercise QUAD.
c NKASES  is the number of integrands to use to exercise QUAD.
c PI      is pi.
c
      integer NEPS, NKASES
      parameter (NEPS=3, NKASES=13)
      real PI
      parameter (PI = 3.1415 92653 58979 32384 62643)
c
c     Arguments of QUAD.
c
      real A, B
      real EPSIL
      real F
      external F
      integer ICHECK, K, NPTS
      real RESULT(8)
c
c     Other variables
c
c ANSWRS  is an array of RESULT(K) found by calling QUAD; indexed by
c         the tolerance index.
c AS, BS  are arrays of lower and upper bounds, indexed by KASE, which
c         is in common /COM/.
c EPSILS  are values to assign to EPSIL, indexed by I.
c I       is a loop inductor and subscript.
c ICHEKS  is an array of ICHECK values returned by QUAD; indexed by
c         the tolerance index.
c NPTSS   is an array of NPTS values returned by QUAD; indexed by
c         the tolerance index.
c
      real ANSWRS(NEPS), AS(NKASES), BS(NKASES), EPSILS(NEPS)
      integer I, ICHEKS(NEPS), NPTSS(NEPS)
c
c     Common Variables
c
c KASE    tells the integrand function what case to work on.
c
      integer KASE
      common /COM/ KASE
c
c     Data
c
      data AS /0.0, -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0,
     1         0.0, 0.0/
      data BS /1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 10.0, pi, 1.0,
     1         1.0, 1.0/
      data EPSILS /1.e-3, 1.e-6, 1.e-8/
c
c     Executable Statements
c
      print 10, epsils
10    format (' eps=',5(:'------- ',1pg9.1,' ------ '))
      print 20, (' ',i=1,neps)
20    format (' KASE',5(:' NPTS     RESULT    CHK ',a1))
      do 100 kase = 1, nkases
        a = as(kase)
        b = bs(kase)
        do 50 i = 1, neps
          epsil = epsils(i)
          call quad (a, b, result, k, epsil, npts, icheck, f)
          answrs(i)=result(k)
          nptss(i)=npts
          icheks(i)=icheck
50      continue
        print 70, kase, (nptss(i),answrs(i),icheks(i),i=1,neps)
70      format (i5,5(i4,1pe16.8,i2,3x))
100   continue
      stop
      end
c     =====     Integrand Function     =================================
      real function F (X)
      real X
      real PI
      parameter (PI = 3.1415 92653 58979 32384 62643)
      real T
      integer KASE
      common /COM/ KASE
      go to (1,2,3,4,5,6,7,8,9,10,11,12,13), kase
1     f=sqrt(x)
      return
2     f=0.92*cosh(x)-cos(x)
      return
3     t=x*x
      f=1.0/(t*(t+1.0)+0.9)
      return
4     f=sqrt(x)**3
      return
5     f=1.0/(1.0+x**4)
      return
6     f=1.0/(1.0+0.5*sin(10.0*pi*x))
      return
7     f=x/(exp(x)-1.0)
      return
8     f=sin(100.0*pi*x)/(pi*x)
      return
9     f=50.0/((2500.0*x*x+1.0)*pi)
      return
10    t=x+x
      f=cos(cos(x)+3.0*(sin(x)+cos(3.0*x)+sin(t))+2.0*cos(t))
      return
11    f=log(x)
      return
12    f=4.0*pi*pi*x*sin(20.0*pi*x)*cos(2.0*pi*x)
      return
13    f=1.0/(1.0+(230.0*x-30.0)**2)
      return
      end
c     =====     Revised QUAD     =======================================
      subroutine QUAD (A, B, RESULT, K, EPSIL, NPTS, ICHECK, F)
c
c     This quadrature program uses formulae due to T. N. L. Patterson,
c     Mathematics of computation, Volume 22, 1968, pages 847-856, as
c     modified by F. T. Krogh and W. V. Snyder, ACM Transactions on
c     Mathematical Software 17, 4 (December 1991) pp 457-461.  It is a
c     functional replacement for Algorithm 468, T. N. L. Patterson,
c     Communications of the ACM 16, 11 (November 1973) 694-699.
c
c     *****     Formal Arguments     ***********************************
c
c Input:
c  A, B   Lower and upper limits of integration, respectively.
c  EPSIL  Relative accuracy required.  When the relative difference of
c         two successive formulae does not exceed EPSIL the last formula
c         computed is taken as the result.
c  F      A FUNCTION subprogram that evaluates the integrand at a given
c         abscissa.  F is invoked F(X).  F MUST BE MENTIONED IN AN
c         EXTERNAL STATEMENT IN THE CALLING PROGRAM.
c Output:
c  RESULT This array, which should be declared to have at least 8
c         elements, holds the results obtained by the 1, 3, 7, etc.
c         point formulae.  The number of formulae computed depends on
c         EPSIL.
c  K      RESULT(K) holds the value of the integral to the specified
c         relative accuracy.
c  NPTS   Reports the number of integrand evaluations.
c  ICHECK On exit normally ICHECK=0.  However if convergence to the
c         accuracy requested is not achieved ICHECK=1 on exit.
c
      real A, B, EPSIL, F
      external F
      real RESULT(8)
      integer K, NPTS, ICHECK
c
c     *****     Local Variables     ************************************
c
c ACUM    is the accumulating estimate of the integral.
c DELTA   is B - A.
c DIFF    is 0.5 * DELTA.
c FH, FL  contain subscripts to indicate where to store integrand
c         samples in WORK.
c FNCVAL  is an integrand sample, or a sum of two symmetrically placed
c         integrand samples.
c IP      is a subscript used to index P.
c J       is a subscript and loop inductor.
c J1, J2  are bounds of the indexes at which integrand samples are
c         stored in WORK.
c JH, JL  are bounds for the loop that accumulates new function values
c         into the integral estimate.
c KH, KL  are bounds for indexes at which integrand samples are
c         retrieved from WORK in order to begin applying a quadrature
c         formula.
c KK      is a loop inductor and subscript.
c KX      is a list of bounds of subscripts into KH and KL.  KH is
c         indexed by K.  The set of indexes from which to retrieve
c         integrand samples is given by the set of bounds in KL and KH
c         indexed by KX(K-1)+1 to KX(K) inclusively.
c P       contains the coefficients necessary to the quadrature
c         formulae.  Their organization is described below in the
c         section on DATA statements.
c PACUM   is the previous value of ACUM.
c X       is the distance of the abscissa from the boundary of the
c         region .
c WORK    is space in which to store function values to be used in the
c         next formula.
c
      double precision ACUM
      real DELTA, DIFF
      integer FH(2:8), FL(2:8)
      real FNCVAL
      integer IP
      integer J, J1, J2, JH, JL
      integer KH(11), KK, KL(11), KX(1:8)
      real P(305)
      double precision PACUM
      real WORK(17), X
c
c     *****     Data Statements     ************************************
c
      data FL /2, 3, 5, 9,12,14, 1/
      data FH /2, 4, 8,16,17,17, 0/
      data KL    /1, 1, 1, 1, 1, 3, 5, 9, 5, 9,12/
      data KH    /1, 2, 4, 8,16, 3, 6,17, 5, 9,17/
      data KX /0, 1, 2, 3, 4, 5,       8,      11/
c
c     In the comments below, F(K,I) refers to the function value
c     computed for the I'th node of the K'th formula.  The abscissae and
c     weights are stored in order according to the distance from the
c     boundary of the region, not from the center.  Since we store
c     1 - |abscissa|, the first "node" coefficient for each formula is
c     the smallest.
c
c     Corrections, nodes and weights for the 3-point formula.
c
c     Correction for F(1,1).
      data P(  1)/-.11111111111111111111e+00/
c     Node and weight for F(2,1).
      data P(  2)/+.22540333075851662296e+00/
      data P(  3)/+.55555555555555555556e+00/
c
c     Corrections, nodes and weights for the 7-point formula.
c
c     Corrections for F(1,1) and F(2,1).
      data P(  4)/+.647209421402969791  e-02/
      data P(  5)/-.928968790944433705  e-02/
c     Nodes and weights for F(3,1-2)
      data P(  6)/+.39508731291979716579e-01/
      data P(  7)/+.10465622602646726519e+00/
      data P(  8)/+.56575625065319744200e+00/
      data P(  9)/+.40139741477596222291e+00/
c
c     Corrections, nodes and weights for the 15-point formula.
c
c     Corrections for F(1,1), F(2,1), F(3,1-2).
      data P( 10)/+.5223046896961622    e-04/
      data P( 11)/+.17121030961750000   e-03/
      data P( 12)/-.724830016153892898  e-03/
      data P( 13)/-.7017801099209042    e-04/
c     Nodes and weights for F(4,1-4).
      data P( 14)/+.61680367872449777899e-02/
      data P( 15)/+.17001719629940260339e-01/
      data P( 16)/+.11154076712774300110e+00/
      data P( 17)/+.92927195315124537686e-01/
      data P( 18)/+.37889705326277359705e+00/
      data P( 19)/+.17151190913639138079e+00/
      data P( 20)/+.77661331357103311837e+00/
      data P( 21)/+.21915685840158749640e+00/
c
c     Corrections, nodes and weights for the 31-point formula.
c
c     Corrections for F(1,1), F(2,1), F(3,1-2), F(4,1-4).
      data P( 22)/+.682166534792        e-08/
      data P( 23)/+.12667409859336      e-06/
      data P( 24)/+.59565976367837165   e-05/
      data P( 25)/+.1392330106826       e-07/
      data P( 26)/-.6629407564902392    e-04/
      data P( 27)/-.704395804282302     e-06/
      data P( 28)/-.34518205339241      e-07/
      data P( 29)/-.814486910996        e-08/
c     Nodes and weights for F(5,1-8).
      data P( 30)/+.90187503233240234038e-03/
      data P( 31)/+.25447807915618744154e-02/
      data P( 32)/+.18468850446259893130e-01/
      data P( 33)/+.16446049854387810934e-01/
      data P( 34)/+.70345142570259943330e-01/
      data P( 35)/+.35957103307129322097e-01/
      data P( 36)/+.16327406183113126449e+00/
      data P( 37)/+.56979509494123357412e-01/
      data P( 38)/+.29750379350847292139e+00/
      data P( 39)/+.76879620499003531043e-01/
      data P( 40)/+.46868025635562437602e+00/
      data P( 41)/+.93627109981264473617e-01/
      data P( 42)/+.66886460674202316691e+00/
      data P( 43)/+.10566989358023480974e+00/
      data P( 44)/+.88751105686681337425e+00/
      data P( 45)/+.11195687302095345688e+00/
c
c     Corrections, nodes and weights for the 63-point formula.
c
c     Corrections for F(1,1), F(2,1), F(3,1-2), F(4,1-4), F(5,1-8).
      data P( 46)/+.371583              e-15/
      data P( 47)/+.21237877            e-12/
      data P( 48)/+.10522629388435      e-08/
      data P( 49)/+.1748029             e-14/
      data P( 50)/+.3475718983017160    e-06/
      data P( 51)/+.90312761725         e-11/
      data P( 52)/+.12558916            e-13/
      data P( 53)/+.54591               e-15/
      data P( 54)/-.72338395508691963   e-05/
      data P( 55)/-.169699579757977     e-07/
      data P( 56)/-.854363907155        e-10/
      data P( 57)/-.12281300930         e-11/
      data P( 58)/-.462334825           e-13/
      data P( 59)/-.42244055            e-14/
      data P( 60)/-.88501               e-15/
      data P( 61)/-.40904               e-15/
c     Nodes and weights for F(6,1-16).
      data P( 62)/+.12711187964238806027e-03/
      data P( 63)/+.36322148184553065969e-03/
      data P( 64)/+.27937406277780409196e-02/
      data P( 65)/+.25790497946856882724e-02/
      data P( 66)/+.11315242452570520059e-01/
      data P( 67)/+.61155068221172463397e-02/
      data P( 68)/+.27817125251418203419e-01/
      data P( 69)/+.10498246909621321898e-01/
      data P( 70)/+.53657141626597094849e-01/
      data P( 71)/+.15406750466559497802e-01/
      data P( 72)/+.89628843042995707499e-01/
      data P( 73)/+.20594233915912711149e-01/
      data P( 74)/+.13609206180630952284e+00/
      data P( 75)/+.25869679327214746911e-01/
      data P( 76)/+.19305946804978238813e+00/
      data P( 77)/+.31073551111687964880e-01/
      data P( 78)/+.26024395564730524132e+00/
      data P( 79)/+.36064432780782572640e-01/
      data P( 80)/+.33709033997521940454e+00/
      data P( 81)/+.40715510116944318934e-01/
      data P( 82)/+.42280428994795418516e+00/
      data P( 83)/+.44914531653632197414e-01/
      data P( 84)/+.51638197305415897244e+00/
      data P( 85)/+.48564330406673198716e-01/
      data P( 86)/+.61664067580126965307e+00/
      data P( 87)/+.51583253952048458777e-01/
      data P( 88)/+.72225017797817568492e+00/
      data P( 89)/+.53905499335266063927e-01/
      data P( 90)/+.83176474844779253501e+00/
      data P( 91)/+.55481404356559363988e-01/
      data P( 92)/+.94365568695340721002e+00/
      data P( 93)/+.56277699831254301273e-01/
c
c     Corrections, nodes and weights for the 127-point formula.
c
c     Corrections for F(3,1), F(4,1-2), F(5,1-3), F(6,1-6).
      data P( 94)/+.1041098             e-15/
      data P( 95)/+.249472054598        e-10/
      data P( 96)/+.55                  e-20/
      data P( 97)/+.290412475995385     e-07/
      data P( 98)/+.367282126           e-13/
      data P( 99)/+.5568                e-18/
      data P(100)/-.871176477376972025  e-06/
      data P(101)/-.8147324267441       e-09/
      data P(102)/-.8830920337          e-12/
      data P(103)/-.18018239            e-14/
      data P(104)/-.70528               e-17/
      data P(105)/-.506                 e-19/
c     Nodes and weights for F(7,1-32).
      data P(106)/+.17569645108401419961e-04/
      data P(107)/+.50536095207862517625e-04/
      data P(108)/+.40120032808931675009e-03/
      data P(109)/+.37774664632698466027e-03/
      data P(110)/+.16833646815926074696e-02/
      data P(111)/+.93836984854238150079e-03/
      data P(112)/+.42758953015928114900e-02/
      data P(113)/+.16811428654214699063e-02/
      data P(114)/+.85042788218938676006e-02/
      data P(115)/+.25687649437940203731e-02/
      data P(116)/+.14628500401479628890e-01/
      data P(117)/+.35728927835172996494e-02/
      data P(118)/+.22858485360294285840e-01/
      data P(119)/+.46710503721143217474e-02/
      data P(120)/+.33362148441583432910e-01/
      data P(121)/+.58434498758356395076e-02/
      data P(122)/+.46269993574238863589e-01/
      data P(123)/+.70724899954335554680e-02/
      data P(124)/+.61679602220407116350e-01/
      data P(125)/+.83428387539681577056e-02/
      data P(126)/+.79659974529987579270e-01/
      data P(127)/+.96411777297025366953e-02/
      data P(128)/+.10025510022305996335e+00/
      data P(129)/+.10955733387837901648e-01/
      data P(130)/+.12348658551529473026e+00/
      data P(131)/+.12275830560082770087e-01/
      data P(132)/+.14935550523164972024e+00/
      data P(133)/+.13591571009765546790e-01/
      data P(134)/+.17784374563501959262e+00/
      data P(135)/+.14893641664815182035e-01/
      data P(136)/+.20891506620015163857e+00/
      data P(137)/+.16173218729577719942e-01/
      data P(138)/+.24251603361948636206e+00/
      data P(139)/+.17421930159464173747e-01/
      data P(140)/+.27857691462990108452e+00/
      data P(141)/+.18631848256138790186e-01/
      data P(142)/+.31701256890892077191e+00/
      data P(143)/+.19795495048097499488e-01/
      data P(144)/+.35772335749024048622e+00/
      data P(145)/+.20905851445812023852e-01/
      data P(146)/+.40059606975775710702e+00/
      data P(147)/+.21956366305317824939e-01/
      data P(148)/+.44550486736806745112e+00/
      data P(149)/+.22940964229387748761e-01/
      data P(150)/+.49231224246628339785e+00/
      data P(151)/+.23854052106038540080e-01/
      data P(152)/+.54086998801016766712e+00/
      data P(153)/+.24690524744487676909e-01/
      data P(154)/+.59102017877011132759e+00/
      data P(155)/+.25445769965464765813e-01/
      data P(156)/+.64259616216846784762e+00/
      data P(157)/+.26115673376706097680e-01/
      data P(158)/+.69542355844328595666e+00/
      data P(159)/+.26696622927450359906e-01/
      data P(160)/+.74932126969651682339e+00/
      data P(161)/+.27185513229624791819e-01/
      data P(162)/+.80410249728889984607e+00/
      data P(163)/+.27579749566481873035e-01/
      data P(164)/+.85957576684743982540e+00/
      data P(165)/+.27877251476613701609e-01/
      data P(166)/+.91554595991628911629e+00/
      data P(167)/+.28076455793817246607e-01/
      data P(168)/+.97181535105025430566e+00/
      data P(169)/+.28176319033016602131e-01/
c
c     Corrections, nodes and weights for the 255-point formula.
c
c     Corrections for F(4,1), F(5,1), F(6,1-2), F(7,1-4).
      data P(170)/+.3326                e-18/
      data P(171)/+.114094770478        e-11/
      data P(172)/+.2952436056970351    e-08/
      data P(173)/+.51608328            e-15/
      data P(174)/-.110177219650597323  e-06/
      data P(175)/-.58656987416475      e-10/
      data P(176)/-.23340340645         e-13/
      data P(177)/-.1248950             e-16/
c     Nodes and weights for F(8,1-64).
      data P(178)/+.24036202515353807630e-05/
      data P(179)/+.69379364324108267170e-05/
      data P(180)/+.56003792945624240417e-04/
      data P(181)/+.53275293669780613125e-04/
      data P(182)/+.23950907556795267013e-03/
      data P(183)/+.13575491094922871973e-03/
      data P(184)/+.61966197497641806982e-03/
      data P(185)/+.24921240048299729402e-03/
      data P(186)/+.12543855319048853002e-02/
      data P(187)/+.38974528447328229322e-03/
      data P(188)/+.21946455040427254399e-02/
      data P(189)/+.55429531493037471492e-03/
      data P(190)/+.34858540851097261500e-02/
      data P(191)/+.74028280424450333046e-03/
      data P(192)/+.51684971993789994803e-02/
      data P(193)/+.94536151685852538246e-03/
      data P(194)/+.72786557172113846706e-02/
      data P(195)/+.11674841174299594077e-02/
      data P(196)/+.98486295992298408193e-02/
      data P(197)/+.14049079956551446427e-02/
      data P(198)/+.12907472045965932809e-01/
      data P(199)/+.16561127281544526052e-02/
      data P(200)/+.16481342421367271240e-01/
      data P(201)/+.19197129710138724125e-02/
      data P(202)/+.20593718329137316189e-01/
      data P(203)/+.21944069253638388388e-02/
      data P(204)/+.25265540247597332240e-01/
      data P(205)/+.24789582266575679307e-02/
      data P(206)/+.30515340497540768229e-01/
      data P(207)/+.27721957645934509940e-02/
      data P(208)/+.36359378430187867480e-01/
      data P(209)/+.30730184347025783234e-02/
      data P(210)/+.42811783890139037259e-01/
      data P(211)/+.33803979910869203823e-02/
      data P(212)/+.49884702478705123440e-01/
      data P(213)/+.36933779170256508183e-02/
      data P(214)/+.57588434808916940190e-01/
      data P(215)/+.40110687240750233989e-02/
      data P(216)/+.65931563842274211999e-01/
      data P(217)/+.43326409680929828545e-02/
      data P(218)/+.74921067092924347640e-01/
      data P(219)/+.46573172997568547773e-02/
      data P(220)/+.84562412844234959360e-01/
      data P(221)/+.49843645647655386012e-02/
      data P(222)/+.94859641186738404810e-01/
      data P(223)/+.53130866051870565663e-02/
      data P(224)/+.10581543166444097714e+00/
      data P(225)/+.56428181013844441585e-02/
      data P(226)/+.11743115975265809315e+00/
      data P(227)/+.59729195655081658049e-02/
      data P(228)/+.12970694445188609414e+00/
      data P(229)/+.63027734490857587172e-02/
      data P(230)/+.14264168911376784347e+00/
      data P(231)/+.66317812429018878941e-02/
      data P(232)/+.15623311732729139895e+00/
      data P(233)/+.69593614093904229394e-02/
      data P(234)/+.17047780536259859981e+00/
      data P(235)/+.72849479805538070639e-02/
      data P(236)/+.18537121234486258656e+00/
      data P(237)/+.76079896657190565832e-02/
      data P(238)/+.20090770903915859819e+00/
      data P(239)/+.79279493342948491103e-02/
      data P(240)/+.21708060588171698360e+00/
      data P(241)/+.82443037630328680306e-02/
      data P(242)/+.23388218069623990928e+00/
      data P(243)/+.85565435613076896192e-02/
      data P(244)/+.25130370638306339718e+00/
      data P(245)/+.88641732094824942641e-02/
      data P(246)/+.26933547875781873867e+00/
      data P(247)/+.91667111635607884067e-02/
      data P(248)/+.28796684463774796540e+00/
      data P(249)/+.94636899938300652943e-02/
      data P(250)/+.30718623022088529711e+00/
      data P(251)/+.97546565363174114611e-02/
      data P(252)/+.32698116976958152079e+00/
      data P(253)/+.10039172044056840798e-01/
      data P(254)/+.34733833458998250389e+00/
      data P(255)/+.10316812330947621682e-01/
      data P(256)/+.36824356228880576959e+00/
      data P(257)/+.10587167904885197931e-01/
      data P(258)/+.38968188628481359983e+00/
      data P(259)/+.10849844089337314099e-01/
      data P(260)/+.41163756555233745857e+00/
      data P(261)/+.11104461134006926537e-01/
      data P(262)/+.43409411457634557737e+00/
      data P(263)/+.11350654315980596602e-01/
      data P(264)/+.45703433350168850951e+00/
      data P(265)/+.11588074033043952568e-01/
      data P(266)/+.48044033846254297801e+00/
      data P(267)/+.11816385890830235763e-01/
      data P(268)/+.50429359208123853983e+00/
      data P(269)/+.12035270785279562630e-01/
      data P(270)/+.52857493412834112307e+00/
      data P(271)/+.12244424981611985899e-01/
      data P(272)/+.55326461233797152625e+00/
      data P(273)/+.12443560190714035263e-01/
      data P(274)/+.57834231337383669993e+00/
      data P(275)/+.12632403643542078765e-01/
      data P(276)/+.60378719394238406082e+00/
      data P(277)/+.12810698163877361967e-01/
      data P(278)/+.62957791204992176986e+00/
      data P(279)/+.12978202239537399286e-01/
      data P(280)/+.65569265840056197721e+00/
      data P(281)/+.13134690091960152836e-01/
      data P(282)/+.68210918793152331682e+00/
      data P(283)/+.13279951743930530650e-01/
      data P(284)/+.70880485148175331803e+00/
      data P(285)/+.13413793085110098513e-01/
      data P(286)/+.73575662758907323806e+00/
      data P(287)/+.13536035934956213614e-01/
      data P(288)/+.76294115441017027278e+00/
      data P(289)/+.13646518102571291428e-01/
      data P(290)/+.79033476175681880523e+00/
      data P(291)/+.13745093443001896632e-01/
      data P(292)/+.81791350324074780175e+00/
      data P(293)/+.13831631909506428676e-01/
      data P(294)/+.84565318851862189130e+00/
      data P(295)/+.13906019601325461264e-01/
      data P(296)/+.87352941562769803314e+00/
      data P(297)/+.13968158806516938516e-01/
      data P(298)/+.90151760340188079791e+00/
      data P(299)/+.14017968039456608810e-01/
      data P(300)/+.92959302395714482093e+00/
      data P(301)/+.14055382072649964277e-01/
      data P(302)/+.95773083523463639678e+00/
      data P(303)/+.14080351962553661325e-01/
      data P(304)/+.98590611358921753738e+00/
      data P(305)/+.14092845069160408355e-01/
c
c     *****     Executable Statements     ******************************
c
      icheck=0
      delta=b-a
      diff=0.5*delta
      ip=1
      jh=0
c
c     Apply 1-point Gauss formula (Midpoint rule).
c
      fncval=f(a+diff)
c     Don't write "0.5*(b+a)" above if the radix of arithmetic isn't 2.
      npts=1
      work(1)=fncval
      acum=fncval*delta
      result(1)=acum
c
      do 40 k = 2, 8
c
c       Go on to the next formula.
c
        pacum=acum
        acum=0.0
c
c       Compute contribution to current estimate due to function
c       values used in previous formulae.
c
        do 20 kk = kx(k-1)+1, kx(k)
          do 10 j = kl(kk), kh(kk)
            acum=acum+dble(p(ip)*work(j))
            ip=ip+1
10        continue
20      continue
c
c       Compute contribution from new function values.
c
        jl=jh+1
        jh=jl+jl-1
        j1=fl(k)
        j2=fh(k)
        do 30 j = jl, jh
          x=p(ip)*diff
          fncval=f(a+x)+f(b-x)
          npts=npts+2
          acum=acum+dble(p(ip+1)*fncval)
          if (j1.le.j2) then
            work(j1)=fncval
            j1=j1+1
          end if
          ip=ip+2
30      continue
        acum=dble(diff)*acum+0.5d0*pacum
        result(k)=acum
        if (abs(result(k)-result(k-1)).le.abs(epsil*result(k))) go to 50
40    continue
      icheck=1
      k=8
50    return
      end