      subroutine erd__rys_roots_weights(NT, NTGQP, NGQP, TVAL, RTS, WTS)
      implicit none
      integer nt, ntgqp, ngqp, i, n, nrts
      double precision tval(nt), rts(ntgqp), wts(ntgqp)
      double precision t, rtstemp(ngqp), wtstemp(ngqp)
      nrts = 0
      do n = 1, nt
        T = TVAL(N)
        call ryspol(NGQP, T, RTSTEMP, WTSTEMP)
        do i = 1, ngqp
            rts(nrts+i) = rtstemp(i)
            wts(nrts+i) = wtstemp(i)
        enddo
        nrts = nrts + ngqp
      enddo
      end

      subroutine ryspol(n, xx, rts, wts, iprint)
      implicit real*8(a-h,o-z)
      parameter (mxquad=25)
      dimension rts(mxquad),wts(mxquad),wrk(mxquad), &
               alpha(0:mxquad-1),beta(0:mxquad-1), &
               xasymp(mxquad)
      parameter (mxaux=100)
      dimension rtsaux(mxaux),wtsaux(mxaux), &
               alpaux(0:mxaux),betaux(0:mxaux), &
               p0(mxaux),p1(mxaux),p2(mxaux), &
               rgrid(mxaux),wgrid(mxaux),nauxs(mxaux)
      character*10 type,code,auxqd
      logical RPprt
!
!         Rys' RMS agreement between asymptotic code and Gautschi code
!         goes below 1d-13 for both roots and weights above these, as
!         it seems Ishida-san was a bit too conservative about this.
!
      data (xasymp(I),I=1,13)/ &
          29.0d+00,  37.0d+00,  43.0d+00,  49.0d+00,  55.0d+00, &
          60.0d+00,  65.0d+00,  71.0d+00,  76.0d+00,  81.0d+00, &
          86.0d+00,  91.0d+00,  96.0d+00/
!
!         Rys auxiliary grid size determined at the xasymp values,
!         eyeballing against using a naux=80 "converged" grid size.
!         It is likely that at lower X values, a smaller auxiliary
!         quadrature would be adequate.
!
      data (nauxs(I),I=1,13)/20,25,30,30,35,  40,40,40,45,50,  50,55,55/
!
!     Program to develop Gaussian quadrature roots and weights
!       assembled/written by Mike Schmidt in November 2005.
!
!     ================================================================
!     for the theory of Orthogonal Polynomials, see the following:
!
!             "Calculation of Gauss Quadrature Rules"
!                 Gene H. Golub, John H. Welsch
!           Mathematics of Computation, 23, 221-230(1969)
!
!     "Algorithm 726: ORTHPOL - A Package of Routines for Generating
!         Orthogonal Polynomials and Gauss-Type Quadrature Rules"
!                         Walter Gautschi
!       ACM Transactions on Mathematical Software  20, 21-62(1994)
!
!         Orthogonal Polynomials: Computation and Approximation
!            W.Gautschi, Oxford University Press, Oxford 2004
!              see Theorem 1.31 and Table 1.1 in this book.
!     
!     The orthogonal polynomial family is defined on an interval [a,b]
!     with respect to a weight function, w(t), 
!         integral p (t) * p (t) * w(t) = delta(i,j) == <p |p >
!           a->b    i       j                             i  j
!     The inner product notation <|> implies the use of the weight
!     function (metric) in the integration.  Any such family will
!     obey the same recursion formula, namely,
!        p   (t) = (t-a )p (t) - b p   (t)      k=0,1,2,...
!         k+1          k  k       k k-1
!     where a and b are called alpha,beta in this program, and are
!                <t*p |p >               <p |p >
!                    k  k                  k  k
!           a  = -----------      b  = -----------
!            k    <p |p >          k   <p   |p   >
!                   k  k                 k-1  k-1
!     We must define the first two members in each family as
!        p  (t) = 0    p (t) = 1
!         -1            0
!     and state that k >= 0 in the a-k formula, while k >=1 in the b-k
!     formula.  Then, stating that b-0 = <p-0|p-0>, e.g. the integral
!     of the weight function only, completes the definition of the
!     orthogonal polynomial.
!
!     Thus, forming alpha,beta recurrence coefficients is tantamount
!     to forming the orthogonal polynomial.  Furthermore, use of the
!     Golub-Walsh library routine 'gauss' will give the roots as the
!     eigenvalues of a tridiagonal matrix (alpha=diagonal, beta=next
!     to diag), and the weights are the first element of each vector.
!     ================================================================
!
      NROOTS = N
      PI=3.14159265358979323846264338327950288419716939937510D+00
      zero= 0.0D+00
      one = 1.0D+00
      two = 2.0D+00
      four= 4.0D+00
!          epsilon is a relative error tolerance in the QL routine
      eps = 1.0D-14
!
      do i=14,mxquad
         nauxs(i) = mxaux
      enddo
!
!        Gauss-Rys quadrature
!        interval [0,+1]         weight exp(-XX*t*t)
!        Introduced by John Rys, Michel Dupuis, Harry King in 1976,
!        to evaluate Cartesian Gaussian integrals in quantum chemistry:
!           M.Dupuis, J.Rys, H.F.King  J.Chem.Phys. 65, 111-116(1976)
!           J.Rys, M.Dupuis, H.F.King  J.Comput.Chem. 4, 154-157(1983)
!        The numerical parameter XX appearing in the weight function
!        depends on the distance between basis function centers, and
!        the exponents in the four Gaussians, 0 <= XX < infinite.
!
!        Rys polynomial computation was first discussed by
!          1. H.F.King, M.Dupuis  J.Comput.Phys. 21, 144-165(1976)
!        Other papers about the computation of Rys polynomials:
!          2. J.D.Augspurger, D.E.Bernholdt, C.E.Dykstra
!                J.Comput.Chem. 11, 972-977(1990)
!          3. K.Ishida  J.Chem.Phys. 95, 5198-5205(1991)
!                       J.Chem.Phys. 98, 2176-2181(1993)
!          4. R.B.Sagar, V.H.Smith Int.J.Quantum Chem. 42, 827-836(1992)
!          5. R.C.Y.Chin  J.Comput.Phys. 99, 321-336(1992)
!          6. P.Carsky, M.Polasek  J.Comput.Phys. 143, 266-277(1998)
!          7. B.I.Schneider, N.Nygaard J.Phys.Chem.A 106, 10773-6(2002)
!        of which, Sagar and Smith's is particularly illuminating.
!
!        It is known from paper 1 that in the large XX limit, the
!        Rys polynomial approaches a scaled Hermite polynomial.
!        See the 2nd paper by Ishida for more information on the
!        large XX limit.
!        The limit XX=0 is manifestly the shifted Legendre case.
!
!        The approach here is to provide a number of different
!        implementations, allowing convenient testing of their
!        speed and numerical accuracy.  
!
!        Sample numerical results can be found in references
!        1, 4, and 7, for NROOTS=5, 13, and 20.  Pictures in 1.
!
!        Using the Golub-Walsh procedure to obtain the roots and
!        weights, from recurrence coefficients -alpha- and -beta-
!        that are generated by the so-called "bootstrap method",
!        using a Stieltjes discretization of the Rys polynomial on
!        an auxiliary grid, is considered the "correct answer":
!                     code=Gautschi.
!
!
!        Choose code from: Dupuis, Sierra, Schmidt, or Gautschi.
!            The last two can print the Rys polynomial, see RPprt.
!            The last two codes need an auxiliary quadrature on
!               our interval [0,1], select from: Fejer or ShiftLeg.
!            The first two codes are dimension limited to NROOTS=13.
!
         RPprt=.false.
         if(XX.ge.xasymp(NROOTS)) then
            call rys_hermite(NROOTS,XX,rts,wts)
         end if
!          
!            For the standard test case, namely NROOTS=13 and XX=10.0,
!            ShiftLeg is much more accurate than Fejer around naux=25,
!            and converges to 14 significant digits about naux=30.
!            Fejer converges to 14 figures only near naux=60 or so.
!            The auxiliary grid must be denser for larger XX, and the
!            values stored in -nauxs- were determined at the XX for
!            which we are able to shift to the asymptotic code.
!            These are therefore overly accurate for smaller XX!
! 
         naux = nauxs(NROOTS)
         call fejer(naux,rtsaux,wtsaux)
!
!           Gautschi's discretized Stieltjes "bootstrapping", which
!           bootstraps the numerical values of the Rys polynomials on
!           the grid, yielding numerically stable alpha,beta values.
!           Of course, the Rys polynomial can be generated using the
!           final alpha,beta values in the recurrence.  If you
!           want to be sure about having the polynomial exactly,
!           run in quadruple precision, at least in RysGen.
!
!DEC$ NOVECTOR 
         do i=1,naux
            t2 = rtsaux(i)*rtsaux(i)
            rgrid(i) = t2
            wgrid(i) = wtsaux(i)*exp(-XX*t2)
         enddo
         call sti(NROOTS,naux,rgrid,wgrid,alpha,beta,ierr,p0,p1,p2)
         if(ierr.ne.0) call cwrite('discretized Stieltjes failed')
         if(RPprt) call RysGen(alpha,beta,NROOTS)
         call gauss(NROOTS,alpha,beta,eps,rts,wts,ierr,wrk)
!
!     the following converts to everyone's convention for Tables.
!
!      do i=1,NROOTS
!         rts(i) = 0.1d0*sqrt(rts(i))
!      enddo
!
!      write(6,920) XX
!      do i=1,NROOTS
!         write(6,940) i,rts(i),wts(i)
!      enddo
!
!  900 format('Enter desired Gauss-Rys NROOTS and X: ',$)
  920 format('Roots and Weights for Gauss-Rys quadrature, at X=',f8.3)
  930 format('Roots and Weights for Gauss-',a,' quadrature.')
  940 format(1x,'i=',i4,1P,' root=',F20.8,' weight=',F20.8)
      end
!
      subroutine gauss(n,alpha,beta,eps,zero,weight,ierr,wrk)
      implicit real*8(a-h,o-z)
      dimension alpha(n),beta(n),zero(n),weight(n),wrk(n)
!
!         this routine is a supplement to the journal article,
!      Algorithm 726: ORTHPOL - A Package of Routines for Generating
!        Orthogonal Polynomials and Gauss-Type Quadrature Rules.
!                         Walter Gautschi, 
!      ACM Transactions on Mathematical Software, 20, 21-62(1994)
!
!            A footnote in Gautschi's article reveals that
!            this particular routine came from Gene Golub.
!
! Given  n  and a measure  dlambda, this routine generates the n-point
! Gaussian quadrature formula
! 
!     integral over supp(dlambda) of f(x)dlambda(x)
!
!        = sum from k=1 to k=n of w(k)f(x(k)) + R(n;f).
! 
! The nodes are returned as  zero(k)=x(k) and the weights as 
! weight(k)=w(k), k=1,2,...,n. The user has to supply the recursion 
! coefficients  alpha(k), beta(k), k=0,1,2,...,n-1, for the measure 
! dlambda. The routine computes the nodes as eigenvalues, and the 
! weights in term of the first component of the respective normalized 
! eigenvectors of the n-th order Jacobi matrix associated with  dlambda.
! It uses a translation and adaptation of the algol procedure  imtql2,
! Numer. Math. 12, 1968, 377-383, by Martin and Wilkinson, as modified 
! by Dubrulle, Numer. Math. 15, 1970, 450. See also Handbook for 
! Autom. Comput., vol. 2 - Linear Algebra, pp.241-248, and the eispack
! routine  imtql2.
!
!        Input:  n - - the number of points in the Gaussian quadrature	
!                      formula; type integer
!                alpha,beta - - arrays of dimension  n  to be filled 
!                      with the values of  alpha(k-1), beta(k-1), k=1,2,
!                      ...,n
!                eps - the relative accuracy desired in the nodes
!                      and weights
!
!        Output: zero- array of dimension  n  containing the Gaussian 
!                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
!                      ...,n
!                weight - array of dimension  n  containing the 
!                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
!                ierr- an error flag equal to  0  on normal return,
!                      equal to  i  if the QR algorithm does not
!                      converge within 30 iterations on evaluating the 
!                      i-th eigenvalue, equal to  -1  if  n  is not in
!                      range, and equal to  -2  if one of the beta's is 
!                      negative.
!
! The array  wrk  is needed for working space.
!
      if(n.lt.1) then
        ierr=-1
        return
      end if
!
      ierr=0
      zero(1)=alpha(1)
      if(beta(1).lt.0.0D+00) then
        ierr=-2
        return
      end if
      weight(1)=beta(1)
      if (n.eq.1) return
!
      weight(1)=1.0D+00
      wrk(n)=0.0D+00
      do 100 k=2,n
        zero(k)=alpha(k)
        if(beta(k).lt.0.0D+00) then
          ierr=-2
          return
        end if
        wrk(k-1)=sqrt(beta(k))
        weight(k)=0.0D+00
  100 continue
!
      do 240 l=1,n
        j=0
!
! Look for a small subdiagonal element.
!
  105   do 110 m=l,n
          if(m.eq.n) go to 120
          if(abs(wrk(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) go to 120
  110   continue
  120   dp=zero(l)
        if(m.eq.l) go to 240
        if(j.eq.30) go to 400
        j=j+1
!
! Form shift.
!
        dg=(zero(l+1)-dp)/(2.0D+00*wrk(l))
        dr=sqrt(dg*dg+1.0D+00)
        dg=zero(m)-dp+wrk(l)/(dg+sign(dr,dg))
        ds=1.0D+00
        dc=1.0D+00
        dp=0.0D+00
        mml=m-l
!
! For i=m-1 step -1 until l do ...
!
        do 200 ii=1,mml
          i=m-ii
          df=ds*wrk(i)
          db=dc*wrk(i)
          if(abs(df).lt.abs(dg)) go to 150
          dc=dg/df
          dr=sqrt(dc*dc+1.0D+00)
          wrk(i+1)=df*dr
          ds=1.0D+00/dr
          dc=dc*ds
          go to 160
  150     ds=df/dg
          dr=sqrt(ds*ds+1.0D+00)
          wrk(i+1)=dg*dr
          dc=1.0D+00/dr
          ds=ds*dc
  160     dg=zero(i+1)-dp
          dr=(zero(i)-dg)*ds+2.0D+00*dc*db
          dp=ds*dr
          zero(i+1)=dg+dp
          dg=dc*dr-db
!
! Form first component of vector.
!
          df=weight(i+1)
          weight(i+1)=ds*weight(i)+dc*df
          weight(i)=dc*weight(i)-ds*df
  200   continue
        zero(l)=zero(l)-dp
        wrk(l)=dg
        wrk(m)=0.0D+00
        go to 105
  240 continue
!         
! Order eigenvalues and eigenvectors.
!
      do 300 ii=2,n
        i=ii-1
        k=i
        dp=zero(i)
        do 260 j=ii,n
          if(zero(j).ge.dp) go to 260
          k=j
          dp=zero(j)
  260   continue
        if(k.eq.i) go to 300
        zero(k)=zero(i)
        zero(i)=dp
        dp=weight(i)
        weight(i)=weight(k)
        weight(k)=dp
  300 continue
      do 310 k=1,n
        weight(k)=beta(1)*weight(k)*weight(k)
  310 continue
      return
!
! Set error - no convergence to an eigenvalue after 30 iterations.
!
  400 ierr=l
      return
      end
!
      subroutine sti(n,ncap,x,w,alpha,beta,ierr,p0,p1,p2)
      implicit real*8(a-h,o-z)
      dimension x(ncap),w(ncap),alpha(n),beta(n), &
                p0(ncap),p1(ncap),p2(ncap)
!
!         this routine is a supplement to the journal article,
!      Algorithm 726: ORTHPOL - A Package of Routines for Generating
!        Orthogonal Polynomials and Gauss-Type Quadrature Rules.
!                         Walter Gautschi, 
!      ACM Transactions on Mathematical Software, 20, 21-62(1994)
!
! This routine applies Stieltjes's procedure (cf. Section 2.1 of
! W. Gautschi,On generating orthogonal polynomials'', SIAM J. Sci.
! Statist. Comput. 3, 1982, 289-317) to generate the recursion
! coefficients  alpha(k), beta(k) , k=0,1,...,n-1, for the discrete
! (monic) orthogonal polynomials associated with the inner product
!
!     (f,g)=sum over k from 1 to ncap of w(k)*f(x(k))*g(x(k)).
!
! The integer  n  must be between  1  and  ncap, inclusive; otherwise,
! there is an error exit with  ierr=1. The results are stored in the
! arrays  alpha, beta; the arrays  p0, p1, p2  are working arrays.
!
! This routine should be used with caution if  n  is relatively close 
! to  ncap, since there is a distinct possibility of numerical 
! instability developing. (See W. Gautschi,Is the recurrence relation 
! for orthogonal polynomials always stable?'', BIT, 1993, to appear.) 
! In that case, the routine  lancz  should be used.
!
      tiny = 1.0d-40
      huge = 1.0d+40
!
      ierr=0
      if(n.le.0 .or. n.gt.ncap) then
        ierr=1
        return
      end if
      nm1=n-1
!
! Compute the first alpha- and beta-coefficient.
!
      sum0=0.0d+00
      sum1=0.0d+00
      do 10 m=1,ncap
        sum0=sum0+w(m)
        sum1=sum1+w(m)*x(m)
   10 continue
      alpha(1)=sum1/sum0
      beta(1)=sum0
      if(n.eq.1) return
!
! Compute the remaining alpha- and beta-coefficients.
!
      do 20 m=1,ncap
        p1(m)=0.0d+00
        p2(m)=1.0d+00
   20 continue
      do 40 k=1,nm1
        sum1=0.0d+00
        sum2=0.0d+00
        do 30 m=1,ncap
!
! The following statement is designed to avoid an overflow condition
! in the computation of  p2(m)  when the weights  w(m)  go to zero 
! faster (and underflow) than the  p2(m)  grow. 
!
          if(w(m).eq.0.0d+00) goto 30
          p0(m)=p1(m)
          p1(m)=p2(m)
          p2(m)=(x(m)-alpha(k))*p1(m)-beta(k)*p0(m)
          t=w(m)*p2(m)*p2(m)
          sum1=sum1+t
          sum2=sum2+t*x(m)
   30   continue
        alpha(k+1)=sum2/sum1
        beta(k+1)=sum1/sum0
        sum0=sum1
   40 continue
!
      return
      end
!
      subroutine fejer(nh,x,w)
      implicit real*8(a-h,o-z)
      dimension x(nh),w(nh)
!
! This routine generates the n-point Fejer quadrature rule.
! Again, this came from Walter Gautschi's package, but was
! modified by MWS to store only the quadrature points on
! the positive half of the usual interval, which is [-1,+1].
!
!         input:   nh  - the number of quadrature nodes at positive x.
!         output:  x,w - arrays of dimension -nh- holding the quadrature
!                        nodes and weights, respectively; the nodes
!                        are ordered increasingly
!
      n=nh+nh
      pi=4.0d+00*atan(1.0d+00)
      np1h=(n+1)/2
      fn=n
!DEC$ NOVECTOR
      do 10 k=1,nh
        x(nh+1-k)=cos(0.5d+00*(2*k-1)*pi/fn)
   10 continue
!     if(2*nh.ne.n) x(np1h)=0.0d+00
      do 30 k=1,np1h
        c1=1.0d+00
        c0=2.0d+00*x(k)*x(k)-1.0d+00
        t=2.0d+00*c0
        sum=c0/3.0d+00
        do 20 m=2,nh
          c2=c1
          c1=c0
          c0=t*c1-c2
          sum=sum+c0/(4*m*m-1)
   20   continue
        w(k)=2.0d+00*(1.0d+00-2.0d+00*sum)/fn
   30 continue
      return
      end
!
      SUBROUTINE RYSSET()
      RETURN
      END
!
      subroutine RysGen(alpha,beta,NROOTS)
      implicit real*8(a-h,o-z)
      dimension alpha(0:*),beta(0:*)
      parameter (mxquad=25)
      dimension RP(0:mxquad,0:mxquad)
      parameter (zero=0.0d+00, one=1.0d+00)
!
      do k=0,NROOTS
         do j=0,NROOTS
            RP(j,k) = zero
         enddo
      enddo
!
!        by definition, the zero-th Rys Polynomial is unity
!
      RP(0,0) = one
!
!        then set the first Rys Polynomial
!
      RP(0,1) = -alpha(0)
      RP(1,1) = one
      if(NROOTS.eq.1) go to 200
!
!        upward recurrence generates the next Rys Polynomial, at k+1
!
      do k=1,NROOTS-1
         do j=0,k
            RP(j,k+1) = -alpha(k)*RP(j,k) - beta(k)*RP(j,k-1)
         enddo
         do j=0,k
            RP(j+1,k+1) = RP(j+1,k+1) + RP(j,k)
         enddo
      enddo
!
  200 continue
!      write(6,900) RP(0,NROOTS),(RP(j,NROOTS),2*j,j=1,NROOTS)
      return
  900 format(1P,e15.7,7x,2(e15.7,'*t**',i2,1x)/(1P,3(e15.7,'*t**',i2,1x)))
      end

     subroutine Rys_Hermite(N,X,TP2,WP2)
     implicit none
     Integer*8 I,J,K,IAD
     Integer N
     Real*8 X, XI, XROOTI, TP2(N), WP2(N)
     Real*8 One, HRoot2(36), HWeigh(36)
     Data HRoot2/0.707106781186547D+00,0.524647623275290D+00, &
     0.165068012388578D+01,0.436077411927616D+00, &
     0.133584907401369D+01,0.235060497367449D+01, &
     0.381186990207322D+00,0.115719371244678D+01, &
     0.198165675669584D+01,0.293063742025724D+01, &
     0.342901327223704D+00,0.103661082978951D+01, &
     0.175668364929988D+01,0.253273167423278D+01, &
     0.343615911883773D+01,0.314240376254359D+00, &
     0.947788391240163D+00,0.159768263515260D+01, &
     0.227950708050105D+01,0.302063702512088D+01, &
     0.388972489786978D+01,0.291745510672562D+00, &
     0.878713787329399D+00,0.147668273114114D+01, &
     0.209518325850771D+01,0.274847072498540D+01, &
     0.346265693360227D+01,0.430444857047363D+01, &
     0.273481046138152d+00,0.822951449144655d+00, &
     0.138025853919888d+01,0.195178799091625d+01, &
     0.254620215784748d+01,0.317699916197995d+01, &
     0.386944790486012d+01,0.468873893930581d+01/
     Data HWeigh/0.886226925452758D+00,0.804914090005512D+00, &
     0.813128354472451D-01,0.724629595224392D+00, &
     0.157067320322856D+00,0.453000990550884D-02, &
     0.661147012558241D+00,0.207802325814891D+00, &
     0.170779830074134D-01,0.199604072211367D-03, &
     0.610862633735325D+00,0.240138611082314D+00, &
     0.338743944554810D-01,0.134364574678123D-02, &
     0.764043285523262D-05,0.570135236262479D+00, &
     0.260492310264161D+00,0.516079856158839D-01, &
     0.390539058462906D-02,0.857368704358785D-04, &
     0.265855168435630D-06,0.536405909712090D+00, &
     0.273105609064246D+00,0.685055342234652D-01, &
     0.785005472645794D-02,0.355092613551923D-03, &
     0.471648435501891D-05,0.862859116812515D-08, &
     0.507929479016613d+00,0.280647458528533d+00, &
     0.838100413989858d-01,0.128803115355099d-01, &
     0.932284008624180d-03,0.271186009253788d-04, &
     0.232098084486521d-06,0.265480747401118d-09/
     Parameter (One=1.0d0)

     XI=ONE/X
     XROOTI= SQRT(XI)
     IAD=(N*(N-1))/2
     DO I=1,N
       IAD=IAD+1
       TP2(I)=HROOT2(IAD)**2*XI
       WP2(I)=HWEIGH(IAD)*XROOTI
     ENDDO
     RETURN
     END


