ORTHPOLq
   
  A QUDRUPLE-PRECISION VERSION OF THE FORTRAN PACKAGE ORTHPOL

     (See Algorithm 726: ORTHPOL -- A package of routines
     for generating orthogonal polynomials and Gauss-type
     quadrature rules, ACM Trans. Math. Software 20 (1994),
     21-62.)


c
c
      real*16 function q1mach(i)
c
c quadruple machine constants
c
      real*16 qmach(5)
      data qmach/.33621031431120935062626778173217526q-4931,
     *.11897314953572317650857593266280070q4933,
     *.96296497219361792652798897129246366q-34,
     *.19259299443872358530559779425849273q-33,
     *.30102999566398119521373889472449302q+00/
      q1mach=qmach(i)
      return
      end
c
c
      subroutine qrecur(n,ipoly,dal,dbe,da,db,iderr)
c
c This is a quadruple-precision version of the routine  recur.
c
      external qgamma
      real*16 dal,dbe,da,db,dlmach,q1mach,dkm1,dalpbe,dt,
     *qlga,dal2,dbe2,qgamma
      dimension da(n),db(n)
      if(n.lt.1) then
        iderr=3
        return
      end if
      dlmach=qlog(q1mach(2))
      iderr=0
      do 10 k=1,n
        da(k)=0.q0
   10 continue
      if(ipoly.eq.1) then
        db(1)=2.q0
        if (n.eq.1) return
        do 20 k=2,n
          dkm1=qfloat(k-1)
          db(k)=1.q0/(4.q0-1.q0/(dkm1*dkm1))
   20   continue
        return
      else if(ipoly.eq.2) then
        da(1)=.5q0
        db(1)=1.q0
        if(n.eq.1) return
        do 30 k=2,n
          da(k)=.5q0
          dkm1=qfloat(k-1)
          db(k)=.25q0/(4.q0-1.q0/(dkm1*dkm1))
   30   continue
        return
      else if(ipoly.eq.3) then
        db(1)=4.q0*qatan(1.q0)
        if(n.eq.1) return
        db(2)=.5q0
        if(n.eq.2) return
        do 40 k=3,n
          db(k)=.25q0
   40   continue
        return
      else if(ipoly.eq.4) then
        db(1)=2.q0*qatan(1.q0)
        if(n.eq.1) return
        do 50 k=2,n
          db(k)=.25q0
   50   continue
        return
      else if(ipoly.eq.5) then
        db(1)=4.q0*qatan(1.q0)
        da(1)=.5q0
        if(n.eq.1) return
        do 60 k=2,n
          db(k)=.25q0
   60   continue
        return
      else if(ipoly.eq.6) then
        if(dal.le.-1.q0 .or. dbe.le.-1.q0) then
          iderr=1
          return
        else
          dalpbe=dal+dbe
          da(1)=(dbe-dal)/(dalpbe+2.q0)
          dt=(dalpbe+1.q0)*qlog(2.q0)+qlga(dal+1.q0)+qlga(dbe+1.q0)-
     *      qlga(dalpbe+2.q0)
          if(dt.gt.dlmach) then
            iderr=2
            db(1)=q1mach(2)
          else
            db(1)=qexp(dt)
          end if
          if(n.eq.1) return
          dal2=dal*dal
          dbe2=dbe*dbe
          da(2)=(dbe2-dal2)/((dalpbe+2.q0)*(dalpbe+4.q0))
          db(2)=4.q0*(dal+1.q0)*(dbe+1.q0)/((dalpbe+3.q0)*(dalpbe+
     *      2.q0)**2)
          if(n.eq.2) return
          do 70 k=3,n
            dkm1=qfloat(k-1)
            da(k)=.25q0*(dbe2-dal2)/(dkm1*dkm1*(1.q0+.5q0*dalpbe/dkm1)
     *        *(1.q0+.5q0*(dalpbe+2.q0)/dkm1))
            db(k)=.25q0*(1.q0+dal/dkm1)*(1.q0+dbe/dkm1)*(1.q0+dalpbe/
     *        dkm1)/((1.q0+.5q0*(dalpbe+1.q0)/dkm1)*(1.q0+.5q0*(dalpbe
     *      -1.q0)/dkm1)*(1.q0+.5q0*dalpbe/dkm1)**2)
   70     continue
          return
        end if
      else if(ipoly.eq.7) then
        if(dal.le.-1.q0) then
          iderr=1
          return
        else
          da(1)=dal+1.q0
          db(1)=qgamma(dal+1.q0,iderr)
          if(iderr.eq.2) db(1)=q1mach(2)
          if(n.eq.1) return
          do 80 k=2,n
            dkm1=qfloat(k-1)
            da(k)=2.q0*dkm1+dal+1.q0
            db(k)=dkm1*(dkm1+dal)
   80     continue
          return
        end if
      else if(ipoly.eq.8) then
        db(1)=qsqrt(4.q0*qatan(1.q0))
        if(n.eq.1) return
        do 90 k=2,n
          db(k)=.5q0*qfloat(k-1)
   90   continue
        return
      else
        iderr=4
      end if
      end

      real*16 function qlga(dx)
      real*16 dbnum,dbden,dx,q1mach,dc,dp,dy,dt,ds
      dimension dbnum(8),dbden(8)
c 
c This routine evaluates the logarithm of the gamma function by a
c combination of recurrence and asymptotic approximation.
c
c The entries in the next data statement are the numerators and
c denominators, respectively, of the quantities B[16]/(16*15),
c B[14]/(14*13),..., B[2]/(2*1), where B[2n] are the Bernoulli
c numbers.
c
      data dbnum/-3.617q3,1.q0,-6.91q2,1.q0,-1.q0,1.q0,-1.q0,1.q0/,
     *     dbden/1.224q5,1.56q2,3.6036q5,1.188q3,1.68q3,1.26q3,3.6q2,
     *1.2q1/
c
c The quantity  dprec  in the next statement is the number of decimal
c digits carried in quadruple-precision floating-point arithmetic.
c
      dprec=-alog10(snglq(q1mach(3)))
      dc=.5q0*qlog(8.q0*qatan(1.q0))
      dp=1.q0
      dy=dx
      y=snglq(dy)
c
c The quantity  y0  below is the threshold value beyond which asymptotic
c evaluation gives sufficient accuracy; see Eq. 6.1.42 in M. Abramowitz
c and I.A. Stegun,Handbook of Mathematical Functions''. The constants 
c are .12118868... = ln(10)/19 and .05390522... = ln(|B[20]|/190)/19.
c
      y0=exp(.121189*dprec+.053905)
   10 if(y.gt.y0) goto 20
      dp=dy*dp
      dy=dy+1.q0
      y=snglq(dy)
      goto 10
   20 dt=1.q0/(dy*dy)
c
c The right-hand side of the next assignment statement is B[18]/(18*17).
c
      ds=4.3867q4/2.44188q5
      do 30 i=1,8
        ds=dt*ds+dbnum(i)/dbden(i)
   30 continue
      qlga=(dy-.5q0)*qlog(dy)-dy+dc+ds/dy-qlog(dp)
      return
      end

      real*16 function qgamma(dx,iderr)
c
c This evaluates the gamma function for real positive  dx, using the
c function subroutine  qlga.
c
      real*16 dx,dlmach,q1mach,dt,qlga
      dlmach=qlog(q1mach(2))
      iderr=0
      dt=qlga(dx)
      if(dt.ge.dlmach) then
        iderr=2
        qgamma=q1mach(2)
        return
      else
        qgamma=qexp(dt)
        return
      end if
      end

c
c
      subroutine qmcdis(n,ncapm,mc,mp,dxp,dyp,qquad,deps,iq,idelta,
     *irout,finld,finrd,dendl,dendr,dxfer,dwfer,dalpha,dbeta,ncap,
     *kount,ierrd,ied,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2)
c
c This is a quadruple-precision version of the routine  mcdis.
c
      real*16 dxp,dyp,deps,dendl,dendr,dxfer,dwfer,dalpha,
     *dbeta,dbe,dx,dw,dxm,dwm,dp0,dp1,dp2
      dimension dxp(*),dyp(*),dendl(mc),dendr(mc),dxfer(ncapm),
     *dwfer(ncapm),dalpha(n),dbeta(n),dbe(n),dx(ncapm),dw(ncapm),
     *dxm(*),dwm(*),dp0(*),dp1(*),dp2(*)
      logical finld,finrd
c
c The arrays  dxp,dyp  are assumed to have dimension  mp  if mp > 0,
c the arrays  dxm,dwm,dp0,dp1,dp2  dimension  mc*ncapm+mp.
c
      if(idelta.le.0) idelta=1
      if(n.lt.1) then
        ierrd=-1
        return
      end if
      incap=1
      kount=-1
      ierr=0
      do 10 k=1,n
        dbeta(k)=0.q0
   10 continue
      ncap=(2*n-1)/idelta
   20 do 30 k=1,n
        dbe(k)=dbeta(k)
   30 continue
      kount=kount+1
      if(kount.gt.1) incap=2**(kount/5)*n
      ncap=ncap+incap
      if(ncap.gt.ncapm) then
        ierrd=ncapm
        return
      end if
      mtncap=mc*ncap
      do 50 i=1,mc
        im1tn=(i-1)*ncap
        if(iq.eq.1) then
          call qquad(ncap,dx,dw,i,ierr)
        else
          call qqgp(ncap,dx,dw,i,ierr,mc,finld,finrd,dendl,dendr,dxfer,
     *      dwfer)
        end if
        if(ierr.ne.0) then
          ierrd=i
          return
        end if
        do 40 k=1,ncap
          dxm(im1tn+k)=dx(k)
          dwm(im1tn+k)=dw(k)
   40   continue
   50 continue
      if(mp.ne.0) then
        do 60 k=1,mp
          dxm(mtncap+k)=dxp(k)
          dwm(mtncap+k)=dyp(k)
   60   continue
      end if
      if(irout.eq.1) then
        call qsti(n,mtncap+mp,dxm,dwm,dalpha,dbeta,ied,dp0,dp1,dp2)
      else
        call qlancz(n,mtncap+mp,dxm,dwm,dalpha,dbeta,ied,dp0,dp1)
      end if
      do 70 k=1,n
        if(qabs(dbeta(k)-dbe(k)).gt.deps*qabs(dbeta(k))) goto 20
   70 continue
      return
      end
c
c
      subroutine qqgp(n,dx,dw,i,ierr,mcd,finld,finrd,dendl,dendr,
     *  dxfer,dwfer)
c
c This is a quadruple-precision version of the routine  qgp. The user
c has to supply the routine
c
c              real*16 function qwf(dx,i),
c
c which evaluates the weight function in quadruple precision at the
c point  dx  on the i-th component interval.
c
      real*16 dx,dw,dendl,dendr,dxfer,dwfer,dphi,dphi1,qwf
      dimension dx(n),dw(n),dendl(mcd),dendr(mcd),dxfer(*),dwfer(*)
      logical finld,finrd
c
c The arrays  dxfer,dwfer  are dimensioned in the routine  qmcdis.
c
      ierr=0
      if(i.eq.1) call qfejer(n,dxfer,dwfer)
      if(i.gt.1 .and. i.lt.mcd) goto 60
      if(mcd.eq.1) then
        if(finld.and.finrd) goto 60
        if(finld) goto 20
        if(finrd) goto 40
        do 10 k=1,n
          call qsymtr(dxfer(k),dphi,dphi1)
          dx(k)=dphi
          dw(k)=dwfer(k)*qwf(dphi,i)*dphi1
   10   continue
        return
      else
        if((i.eq.1.and.finld).or.(i.eq.mcd.and.finrd)) goto 60
        if(i.eq.1) goto 40
      end if
   20 do 30 k=1,n
        call qtr(dxfer(k),dphi,dphi1)
        dx(k)=dendl(mcd)+dphi
        dw(k)=dwfer(k)*qwf(dx(k),mcd)*dphi1
   30 continue
      return
   40 do 50 k=1,n
        call qtr(-dxfer(k),dphi,dphi1)
        dx(k)=dendr(1)-dphi
        dw(k)=dwfer(k)*qwf(dx(k),1)*dphi1
   50 continue
      return
   60 do 70 k=1,n
        dx(k)=.5d0*((dendr(i)-dendl(i))*dxfer(k)+dendr(i)+dendl(i))
        dw(k)=.5d0*(dendr(i)-dendl(i))*dwfer(k)*qwf(dx(k),i)
   70 continue
      return
      end

      subroutine qsymtr(dt,dphi,dphi1)
c
c This is a quadruple-precision version of  symtr.
c
      real*16 dt,dphi,dphi1,dt2
      dt2=dt*dt
      dphi=dt/(1.q0-dt2)
      dphi1=(dt2+1.q0)/(dt2-1.q0)**2
      return
      end

      subroutine qtr(dt,dphi,dphi1)
c
c This is a quadruple-precision version of  tr.
c
      real*16 dt,dphi,dphi1
      dphi=(1.q0+dt)/(1.q0-dt)
      dphi1=2.q0/(dt-1.q0)**2
      return
      end

      subroutine qfejer(n,dx,dw)
c
c This is a quadruple-precision version of  fejer.
c
      real*16 dx,dw,dpi,dn,dc1,dc0,dt,dsum,dc2
      dimension dx(n),dw(n)
      dpi=4.q0*qatan(1.q0)
      nh=n/2
      np1h=(n+1)/2
      dn=qfloat(n)
      do 10 k=1,nh
        dx(n+1-k)=qcos(.5q0*qfloat(2*k-1)*dpi/dn)
        dx(k)=-dx(n+1-k)
   10 continue
      if(2*nh.ne.n) dx(np1h)=0.q0
      do 30 k=1,np1h
        dc1=1.q0
        dc0=2.q0*dx(k)*dx(k)-1.q0
        dt=2.q0*dc0
        dsum=dc0/3.q0
        do 20 m=2,nh
          dc2=dc1
          dc1=dc0
          dc0=dt*dc1-dc2
          dsum=dsum+dc0/qfloat(4*m*m-1)
   20   continue
        dw(k)=2.q0*(1.q0-2.q0*dsum)/dn
        dw(n+1-k)=dw(k)
   30 continue
      return
      end

c
c
      subroutine qsti(n,ncap,dx,dw,dalpha,dbeta,ierr,dp0,dp1,dp2)
c
c This is a quadruple-precision version of the routine  sti.
c
      real*16 dx,dw,dalpha,dbeta,dp0,dp1,dp2,dtiny,q1mach,
     *dhuge,dsum0,dsum1,dsum2,dt
      dimension dx(ncap),dw(ncap),dalpha(n),dbeta(n),dp0(ncap),
     *dp1(ncap),dp2(ncap)
      dtiny=10.q0*q1mach(1)
      dhuge=.1q0*q1mach(2)
      ierr=0
      if(n.le.0 .or. n.gt.ncap) then
        ierr=1
        return
      end if
      nm1=n-1
      dsum0=0.q0
      dsum1=0.q0
      do 10 m=1,ncap
        dsum0=dsum0+dw(m)
        dsum1=dsum1+dw(m)*dx(m)
   10 continue
      dalpha(1)=dsum1/dsum0
      dbeta(1)=dsum0
      if(n.eq.1) return
      do 20 m=1,ncap
        dp1(m)=0.q0
        dp2(m)=1.q0
   20 continue
      do 40 k=1,nm1
        dsum1=0.q0
        dsum2=0.q0
        do 30 m=1,ncap
          if(dw(m).eq.0.q0) goto 30
          dp0(m)=dp1(m)
          dp1(m)=dp2(m)
          dp2(m)=(dx(m)-dalpha(k))*dp1(m)-dbeta(k)*dp0(m)
          if(qabs(dp2(m)).gt.dhuge .or. qabs(dsum2).gt.dhuge) then
            ierr=k
            return
          end if
          dt=dw(m)*dp2(m)*dp2(m)
          dsum1=dsum1+dt
          dsum2=dsum2+dt*dx(m)
   30   continue
        if(qabs(dsum1).lt.dtiny) then
          ierr=-k
          return
        end if
        dalpha(k+1)=dsum2/dsum1
        dbeta(k+1)=dsum1/dsum0
        dsum0=dsum1
   40 continue
      return
      end

c
c
      subroutine qlancz(n,ncap,dx,dw,dalpha,dbeta,ierr,dp0,dp1)
c
c This is a quadruple-precision version of the routine  lancz.
c
      real*16 dx(ncap),dw(ncap),dalpha(n),dbeta(n),
     *dp0(ncap),dp1(ncap),dpi,dgam,dsig,dt,dxlam,drho,dtmp,
     *dtsig,dtk
      if(n.le.0 .or. n.gt.ncap) then
        ierr=1
        return
      else
        ierr=0
      end if
      do 10 i=1,ncap
        dp0(i)=dx(i)
        dp1(i)=0.q0
   10 continue
      dp1(1)=dw(1)
      do 30 i=1,ncap-1
        dpi=dw(i+1)
        dgam=1.q0
        dsig=0.q0
        dt=0.q0
        dxlam=dx(i+1)
        do 20 k=1,i+1
          drho=dp1(k)+dpi
          dtmp=dgam*drho
          dtsig=dsig
          if(drho.le.0.q0) then
            dgam=1.q0
            dsig=0.q0
          else
            dgam=dp1(k)/drho
            dsig=dpi/drho
          end if
          dtk=dsig*(dp0(k)-dxlam)-dgam*dt
          dp0(k)=dp0(k)-(dtk-dt)
          dt=dtk
          if(dsig.le.0.q0) then
            dpi=dtsig*dp1(k)
          else
            dpi=(dt**2)/dsig
          end if
          dtsig=dsig
          dp1(k)=dtmp
   20   continue
   30 continue
      do 40 k=1,n
        dalpha(k)=dp0(k)
        dbeta(k)=dp1(k)
   40 continue
      return
      end 

c
c
      subroutine qcheb(n,da,db,dnu,dalpha,dbeta,ds,iderr,ds0,ds1,ds2)
c
c This is a quadruple-precision version of the routine  cheb.
c
      real*16 da,db,dnu,dalpha,dbeta,ds,ds0,ds1,ds2,dtiny,
     *q1mach,dhuge
      dimension da(*),db(*),dnu(*),dalpha(n),dbeta(n),ds(n),
     *ds0(*),ds1(*),ds2(*)
c
c The arrays  da,db  are assumed to have dimension  2*n-1, the arrays
c dnu,ds0,ds1,ds2  dimension  2*n.
c
      nd=2*n
      dtiny=10.q0*q1mach(1)
      dhuge=.1q0*q1mach(2)
      iderr=0
      if(qabs(dnu(1)).lt.dtiny) then
        iderr=1
        return
      end if
      if(n.lt.1) then
        iderr=2
        return
      end if
      dalpha(1)=da(1)+dnu(2)/dnu(1)
      dbeta(1)=dnu(1)
      if(n.eq.1) return
      ds(1)=dnu(1)
      do 10 l=1,nd
        ds0(l)=0.q0
        ds1(l)=dnu(l)
   10 continue
      do 40 k=2,n
        lk=nd-k+1
        do 20 l=k,lk
          ds2(l)=ds1(l+1)-(dalpha(k-1)-da(l))*ds1(l)-dbeta(k-1)*ds0(l)
     *      +db(l)*ds1(l-1)
        if(l.eq.k) ds(k)=ds2(k)
   20   continue
        if(qabs(ds(k)).lt.dtiny) then
          iderr=-(k-1)
          return
        else if(qabs(ds(k)).gt.dhuge) then
          iderr=k-1
          return
        end if
        dalpha(k)=da(k)+(ds2(k+1)/ds2(k))-(ds1(k)/ds1(k-1))
        dbeta(k)=ds2(k)/ds1(k-1)
        do 30 l=k,lk
          ds0(l)=ds1(l)
          ds1(l)=ds2(l)
   30   continue
   40 continue
      return
      end

c
c
      program qtest1
c
c
      real*16 qnu,q,q0,qrr,qa,qb,qalpha,qbeta,qs,qs0,
     *qs1,qs2,qkk2,qeps,q1mach,qk2
      dimension qnu(160),q(80),q0(80),qrr(80),qa(159),
     *qb(159),qalpha(80),qbeta(80),qs(80),qs0(160),
     *qs1(160),qs2(160),qkk2(4)
      data qkk2/.1q0,.5q0,.9q0,.999q0/
c
c This test generates the first n=80 beta-coefficients in the 
c recurrence relation for the orthogonal polynomials relative to 
c the weight function
c
c         ((1-k2*x**2)*(1-x**2))**(-1/2)  on (-1,1)
c
c for k2=.1,.5,.9,.999, in quadruple precision, using
c modified Chebyshev moments. 
c
      qeps=q1mach(3)
      n=80
      ndm1=2*n-1
      do 30 ik=1,4
        qk2=qkk2(ik)
c
c Compute the modified Chebyshev moments using Eqs.(3.7)
c of the companion paper. 
c
        call qmm1(n,qeps,qk2,qnu,ierr,q,q0,qrr) 
        if(ierr.ne.0) then
          write(*,2) ierr,qk2
    2     format(/5x,'ierr in qmm1 = ',i1,
     *      '  for k2 = ',d12.4/)
          goto 30
        end if 
c
c Generate the recursion coefficients for the polynomials 
c defining the modified moments.
c
        call qrecur(ndm1,3,0.q0,0.q0,qa,qb,ierr)
c
c Compute the desired recursion coefficients by means of the modified
c Chebyshev algorithm; for the latter, see, e.g., Section 2.4 of
c W. Gautschi, ``On generating orthogonal polynomials'', SIAM J. Sci.
c Statist. Comput. 3, 1982, 289-317.
c
        call qcheb(n,qa,qb,qnu,qalpha,qbeta,qs,ierr,qs0,qs1,qs2)
        do 20 k=1,n
          write(*,8) qbeta(k)
    8     format(1x,d33.25)
   20   continue
   30 continue
      stop
      end

      subroutine qmm1(n,qeps,qk2,qnu,ierr,q,q0,qrr)
c
c This is a quadruple-precision version of the routine  fmm.
c
      real*16 qeps,qk2,qnu(*),q(n),q0(n),qrr(n),qpi,qq,
     *qq1,qr,qs,qn1,qc0,qc
c
c The array  qnu  is assumed to have dimension  2*n.
c
      ierr=0
      nd=2*n
      ndm1=nd-1
      qpi=4.q0*qatan(1.q0)
      qq=qk2/(2.q0-qk2+2.q0*qsqrt(1.q0-qk2))
      qq1=(1.q0+qq*qq)/qq
      do 10 k=1,n
        q(k)=0.q0
   10 continue
      nud=nd
   20 nud=nud+10
      do 30 k=1,n
        q0(k)=q(k)
   30 continue
      if(nud.gt.2000) then
        ierr=1
        return
      end if
      qr=0.q0
      qs=0.q0
      do 40 k=1,nud
        n1=nud-k+1
        qn1=qreal(n1)
        qr=-(qn1-.5q0)/(qn1*qq1+(qn1+.5q0)*qr)
        qs=qr*(2.q0+qs)
        if(n1.le.n) qrr(n1)=qr
   40 continue
      qc0=1.q0/(1.q0+qs)
      q(1)=qrr(1)*qc0
      if(n.gt.1) then
        do 50 k=2,n
          q(k)=qrr(k)*q(k-1)
   50   continue
      end if
      do 60 k=1,n
        if(qabs(q(k)-q0(k)).gt.qeps*qabs(q(k))) goto 20
   60 continue
      qnu(1)=qpi*qc0
      if(n.eq.1) return
      qnu(2)=0.q0
      if(n.eq.2) return
      qc=2.q0*qpi
      do 70 k=3,ndm1,2
        k1=(k-1)/2
        qc=-.25q0*qc
        qnu(k)=qc*q(k1)
        qnu(k+1)=0.q0
   70 continue
      end

c
c
      program qtest2
c
c
      real*16 qa,qb,qnu,qalpha,qbeta,qs,qs0,qs1,qs2,qsigma
      dimension qa(199),qb(199),qnu(200),qalpha(100),qbeta(100),
     *qs(100),qs0(200),qs1(200),qs2(200)
      logical intexp
c
c This test generates the first n=100 recursion coefficients for the
c orthogonal polynomials relative to the weight function
c
c        (x**sigma)*ln(1/x)  on (0,1],  sigma = -.5, 0, .5,
c
c using modified Legendre moments.  It prints the quadruple-
c precision values of the coefficients.
c
c
c Generate the recursion coefficients for the polynomials 
c defining the modified moments.
c
      n=100
      ndm1=2*n-1
      call qrecur(ndm1,2,0.q0,0.q0,qa,qb,ierr)
      do 30 is=1,3
        qsigma=-.5q0+.5q0*qreal(is-1)
        if(is.eq.2) then
          intexp=.true.
        else
          intexp=.false.
        end if
c
c Compute the modified moments using Eqs. (3.12) of the 
c companion paper. 
c
        call qmm2(n,intexp,qsigma,qnu)
c
c Compute the desired recursion coefficients by means of the modified
c Chebyshev algorithm; for the latter, see, e.g., Section 2.4 of
c W. Gautschi, ``On generating orthogonal polynomials'', SIAM J. Sci.
c Statist. Comput. 3, 1982, 289-317.
c
        call qcheb(n,qa,qb,qnu,qalpha,qbeta,qs,ierr,qs0,qs1,qs2)
        do 25 k=1,n
          write(*,2) qalpha(k),qbeta(k)
    2     format(1x,2d33.25)
   25   continue
   30 continue
      stop
      end

      subroutine qmm2(n,intexp,qsigma,qnu)
      real*16 qsigma,qnu(*),qsigp1,qc,qk,qp,qs,qi,qq,qiq,qkm1
      logical intexp
c
c The array  qnu  is assumed to have dimension  2*n.
c
      nd=2*n
      qsigp1=qsigma+1.q0
      isigma=idint(qsigma)
      isigp1=isigma+1
      isigp2=isigma+2
      isigp3=isigma+3
      if(intexp .and. isigp1.lt.nd) then
        kmax=isigp1
      else
        kmax=nd
      end if
      qc=1.q0
      do 20 k=1,kmax
        km1=k-1
        qk=qreal(k)
        qp=1.q0
        qs=1.q0/qsigp1
        if(kmax.gt.1) then
          do 10 i=1,km1
            qi=qreal(i)
            qp=(qsigp1-qi)*qp/(qsigp1+qi)
            qs=qs+1.q0/(qsigp1+qi)-1.q0/(qsigp1-qi)
   10     continue
        end if
        qnu(k)=qc*qs*qp/qsigp1
        qc=qk*qc/(4.q0*qk-2.q0)
   20 continue
      if(.not.intexp .or. isigp1.ge.nd) return
      qq=-.5q0
      if(isigma.gt.0) then
        do 30 iq=1,isigma
          qiq=qreal(iq)
          qq=qiq*qiq*qq/((2.q0*qiq+1.q0)*(2.q0*qiq+2.q0))
   30   continue
      end if
      qnu(isigp2)=qc*qq
      if(isigp2.eq.nd) return
      do 40 k=isigp3,nd
        km1=k-1
        qkm1=qreal(km1)
        qnu(k)=-qkm1*(qkm1-qsigp1)*qnu(km1)/((4.q0*qkm1-2.q0)*
     *    (qkm1+qsigp1))
   40 continue
      return
      end

c
c
      subroutine qmccheb(n,ncapm,mcd,mp,dxp,dyp,qquad,deps,iq,
     *idelta,finld,finrd,dendl,dendr,dxfer,dwfer,da,db,dnu,dalpha,
     *dbeta,ncap,kount,ierrd,dbe,dx,dw,dxm,dwm,ds,ds0,ds1,ds2)
c
c This is a quadruple-precision version of the routine  mccheb.
c
      real*16 dxp,dyp,deps,dendl,dendr,dxfer,dwfer,da,db,
     *dnu,dalpha,dbeta,dbe,dx,dw,dxm,dwm,ds,ds0,ds1,ds2,dsum,dp1,
     *dp,dpm1
      dimension dxp(*),dyp(*),dendl(mcd),dendr(mcd),dxfer(ncapm),
     *dwfer(ncapm),da(*),db(*),dnu(*),dalpha(n),dbeta(n),dbe(n),
     *dx(ncapm),dw(ncapm),dxm(*),dwm(*),ds(n),ds0(*),ds1(*),ds2(*)
      logical finld,finrd
c
c The arrays  dxp,dyp  are assumed to have dimension  mp  if mp > 0,
c the arrays  da,db  dimension 2*n-1, the arrays  dnu,ds0,ds1,ds2
c dimension  2*n, and the arrays  dxm,dwm  dimension  mc*ncapm+mp.
c
      nd=2*n
      if(idelta.le.0) idelta=1
      if(n.lt.1) then
        ierrd=-1
        return
      end if
      incap=1
      kount=-1
      ierrd=0
      do 10 k=1,n
        dbeta(k)=0.q0
   10 continue
      ncap=(nd-1)/idelta
   20 do 30 k=1,n
        dbe(k)=dbeta(k)
   30 continue
      kount=kount+1
      if(kount.gt.1) incap=2**(kount/5)*n
      ncap=ncap+incap
      if(ncap.gt.ncapm) then
        ierrd=ncapm
        return
      end if
      mtncap=mcd*ncap
      do 50 i=1,mcd
        im1tn=(i-1)*ncap
        if(iq.eq.1) then
          call qquad(ncap,dx,dw,i,ierr)
        else
          call qqgp(ncap,dx,dw,i,ierr,mcd,finld,finrd,dendl,dendr,
     *      dxfer,dwfer)
        end if
        if(ierr.ne.0) then
          ierrd=i
          return
        end if
        do 40 k=1,ncap
          dxm(im1tn+k)=dx(k)
          dwm(im1tn+k)=dw(k)
   40   continue
   50 continue
      if(mp.ne.0) then
        do 60 k=1,mp
          dxm(mtncap+k)=dxp(k)
          dwm(mtncap+k)=dyp(k)
   60   continue
      end if
      mtnpmp=mtncap+mp
      do 90 k=1,nd
        km1=k-1
        dsum=0.q0
        do 80 i=1,mtnpmp
          dp1=0.q0
          dp=1.q0
          if(k.gt.1) then
            do 70 l=1,km1
              dpm1=dp1
              dp1=dp
              dp=(dxm(i)-da(l))*dp1-db(l)*dpm1
   70       continue
          end if
          dsum=dsum+dwm(i)*dp
   80   continue
        dnu(k)=dsum
   90 continue
      call qcheb(n,da,db,dnu,dalpha,dbeta,ds,ierr,ds0,ds1,ds2)
      do 100 k=1,n
        if(qabs(dbeta(k)-dbe(k)).gt.deps*qabs(dbeta(k))) goto 20
  100 continue
      return
      end

c
c
      subroutine qgauss(n,dalpha,dbeta,deps,dzero,dweigh,ierr,de)
c
c This is a quadruple-precision version of the routine  gauss.
c
      real*16 dalpha,dbeta,deps,dzero,dweigh,de,dp,dg,dr,
     *ds,dc,df,db
      dimension dalpha(n),dbeta(n),dzero(n),dweigh(n),de(n)
      if(n.lt.1) then
        ierr=-1
        return
      end if
      ierr=0
      dzero(1)=dalpha(1)
      if(dbeta(1).lt.0.q0) then
        ierr=-2
        return
      end if
      dweigh(1)=dbeta(1)
      if (n.eq.1) return
      dweigh(1)=1.q0
      de(n)=0.q0
      do 100 k=2,n
        dzero(k)=dalpha(k)
        if(dbeta(k).lt.0.q0) then
          ierr=-2
          return
        end if
        de(k-1)=qsqrt(dbeta(k))
        dweigh(k)=0.q0
  100 continue
      do 240 l=1,n
        j=0
  105   do 110 m=l,n
          if(m.eq.n) goto 120
          if(qabs(de(m)).le.deps*(qabs(dzero(m))+qabs(dzero(m+1)))) 
     *      goto 120
  110   continue
  120   dp=dzero(l)
        if(m.eq.l) goto 240
        if(j.eq.30) goto 400
        j=j+1
        dg=(dzero(l+1)-dp)/(2.q0*de(l))
        dr=qsqrt(dg*dg+1.q0)
        dg=dzero(m)-dp+de(l)/(dg+qsign(dr,dg))
        ds=1.q0
        dc=1.q0
        dp=0.q0
        mml=m-l
        do 200 ii=1,mml
          i=m-ii
          df=ds*de(i)
          db=dc*de(i)
          if(qabs(df).lt.qabs(dg)) goto 150
          dc=dg/df
          dr=qsqrt(dc*dc+1.q0)
          de(i+1)=df*dr
          ds=1.q0/dr
          dc=dc*ds
          goto 160
  150     ds=df/dg
          dr=qsqrt(ds*ds+1.q0)
          de(i+1)=dg*dr
          dc=1.q0/dr
          ds=ds*dc
  160     dg=dzero(i+1)-dp
          dr=(dzero(i)-dg)*ds+2.q0*dc*db
          dp=ds*dr
          dzero(i+1)=dg+dp
          dg=dc*dr-db
          df=dweigh(i+1)
          dweigh(i+1)=ds*dweigh(i)+dc*df
          dweigh(i)=dc*dweigh(i)-ds*df
  200   continue
        dzero(l)=dzero(l)-dp
        de(l)=dg
        de(m)=0.q0
        goto 105
  240 continue
      do 300 ii=2,n
        i=ii-1
        k=i
        dp=dzero(i)
        do 260 j=ii,n
          if(dzero(j).ge.dp) goto 260
          k=j
          dp=dzero(j)
  260   continue
        if(k.eq.i) goto 300
        dzero(k)=dzero(i)
        dzero(i)=dp
        dp=dweigh(i)
        dweigh(i)=dweigh(k)
        dweigh(k)=dp
  300 continue
      do 310 k=1,n
        dweigh(k)=dbeta(1)*dweigh(k)*dweigh(k)
  310 continue
      return
  400 ierr=l
      return
      end

c
c
      subroutine qradau(n,dalpha,dbeta,dend,dzero,dweigh,ierr,de,
     *da,db)
c
c This is a quadruple-precision version of the routine  radau.
c
      real*16 dend,depsma,dp0,dp1,dpm1,dalpha(*),dbeta(*),
     *dzero(*),dweigh(*),de(*),da(*),db(*),q1mach
c
c The arrays  dalpha,dbeta,dzero,dweigh,de,da,db  are assumed to have
c dimension  n+1.
c
      depsma=q1mach(3)
c
c depsma is the machine quadruple precision.
c
      np1=n+1
      do 10 k=1,np1
        da(k)=dalpha(k)
        db(k)=dbeta(k)
   10 continue
      dp0=0.q0
      dp1=1.q0
      do 20 k=1,n
        dpm1=dp0
        dp0=dp1
        dp1=(dend-da(k))*dp0-db(k)*dpm1
   20 continue
      da(np1)=dend-db(np1)*dp0/dp1
      call qgauss(np1,da,db,depsma,dzero,dweigh,ierr,de)
      return
      end

c
c

      subroutine qlob(n,dalpha,dbeta,dleft,dright,dzero,dweigh,
     *ierr,de,da,db)
c
c This is a quadruple-precision version of the routine  lob.
c
      real*16 dleft,dright,depsma,dp0l,dp0r,dp1l,dp1r,dpm1l,
     *dpm1r,ddet,dalpha(*),dbeta(*),dzero(*),dweigh(*),de(*),da(*),
     *db(*),q1mach
c
c The arrays  dalpha,dbeta,dzero,dweigh,de,da,db  are assumed to have
c dimension  n+2.
c
      depsma=q1mach(3)
c
c depsma is the machine quadruple precision.
c
      np1=n+1
      np2=n+2
      do 10 k=1,np2
        da(k)=dalpha(k)
        db(k)=dbeta(k)
   10 continue
      dp0l=0.q0
      dp0r=0.q0
      dp1l=1.q0
      dp1r=1.q0
      do 20 k=1,np1
        dpm1l=dp0l
        dp0l=dp1l
        dpm1r=dp0r
        dp0r=dp1r
        dp1l=(dleft-da(k))*dp0l-db(k)*dpm1l
        dp1r=(dright-da(k))*dp0r-db(k)*dpm1r
   20 continue
      ddet=dp1l*dp0r-dp1r*dp0l
      da(np2)=(dleft*dp1l*dp0r-dright*dp1r*dp0l)/ddet
      db(np2)=(dright-dleft)*dp1l*dp1r/ddet
      call qgauss(np2,da,db,depsma,dzero,dweigh,ierr,de)
      return
      end

c
c
      subroutine qchri(n,iopt,da,db,dx,dy,dhr,dhi,dalpha,dbeta,ierr)
c
c This is a quadruple-precision version of the routine  chri.
c
      real*16 da,db,dx,dy,dhr,dhi,dalpha,dbeta,deps,q1mach,
     *de,dq,ds,dt,deio,dd,der,dei,deroo,deioo,dso,dero,deoo,deo,du,dc,
     *dc0,dgam,dcm1,dp2,deh,dqh,dalh
      dimension da(*),db(*),dalpha(n),dbeta(n)
c
c The arrays  da,db  are assumed to have dimension  n+1.
c
      deps=5.q0*q1mach(3)
      ierr=0
      if(n.lt.2) then
        ierr=1
        return
      end if
      if(iopt.eq.1) then
        de=0.q0
        do 10 k=1,n
          dq=da(k)-de-dx
          dbeta(k)=dq*de
          de=db(k+1)/dq
          dalpha(k)=dx+dq+de
   10   continue
        dbeta(1)=db(1)*(da(1)-dx)
        return
      else if(iopt.eq.2) then
        ds=dx-da(1)
        dt=dy
        deio=0.q0
        do 20 k=1,n
          dd=ds*ds+dt*dt
          der=-db(k+1)*ds/dd
          dei=db(k+1)*dt/dd
          ds=dx+der-da(k+1)
          dt=dy+dei
          dalpha(k)=dx+dt*der/dei-ds*dei/dt
          dbeta(k)=dt*deio*(1.q0+(der/dei)**2)
          deio=dei
   20   continue
        dbeta(1)=db(1)*(db(2)+(da(1)-dx)**2+dy*dy)
        return
      else if(iopt.eq.3) then
        dt=dy
        deio=0.q0
        do 30 k=1,n
          dei=db(k+1)/dt
          dt=dy+dei
          dalpha(k)=0.q0
          dbeta(k)=dt*deio
          deio=dei
   30   continue
        dbeta(1)=db(1)*(db(2)+dy*dy)
        return
      else if(iopt.eq.4) then
        dalpha(1)=dx-db(1)/dhr
        dbeta(1)=-dhr
        dq=-db(1)/dhr
        do 40 k=2,n
          de=da(k-1)-dx-dq
          dbeta(k)=dq*de
          dq=db(k)/de
          dalpha(k)=dq+de+dx
   40   continue
        return
      else if(iopt.eq.5) then
        nm1=n-1
        dd=dhr*dhr+dhi*dhi
        deroo=da(1)-dx+db(1)*dhr/dd
        deioo=-db(1)*dhi/dd-dy
        dalpha(1)=dx+dhr*dy/dhi
        dbeta(1)=-dhi/dy
        dalpha(2)=dx-db(1)*dhi*deroo/(dd*deioo)+dhr*deioo/dhi
        dbeta(2)=dy*deioo*(1.q0+(dhr/dhi)**2)
        if(n.eq.2) return
        dso=db(2)/(deroo**2+deioo**2)
        dero=da(2)-dx-dso*deroo
        deio=dso*deioo-dy
        dalpha(3)=dx+deroo*deio/deioo+dso*deioo*dero/deio
        dbeta(3)=-db(1)*dhi*deio*(1.q0+(deroo/deioo)**2)/dd
        if(n.eq.3) return
        do 50 k=3,nm1
          ds=db(k)/(dero**2+deio**2)
          der=da(k)-dx-ds*dero
          dei=ds*deio-dy
          dalpha(k+1)=dx+dero*dei/deio+ds*deio*der/dei
          dbeta(k+1)=dso*deioo*dei*(1.q0+(dero/deio)**2)
          deroo=dero
          deioo=deio
          dero=der
          deio=dei
          dso=ds
   50   continue
        return
      else if(iopt.eq.6) then
        nm1=n-1
        deoo=-db(1)/dhi-dy
        deo=db(2)/deoo-dy
        dalpha(1)=0.q0
        dbeta(1)=-dhi/dy
        dalpha(2)=0.q0
        dbeta(2)=dy*deoo
        if(n.eq.2) return
        dalpha(3)=0.q0
        dbeta(3)=-db(1)*deo/dhi
        if(n.eq.3) return
        do 60 k=3,nm1
          de=db(k)/deo-dy
          dbeta(k+1)=db(k-1)*de/deoo
          dalpha(k+1)=0.q0
          deoo=deo
          deo=de
   60   continue
        return
      else if(iopt.eq.7) then
        du=0.q0
        dc=1.q0
        dc0=0.q0
        do 70 k=1,n
          dgam=da(k)-dx-du
          dcm1=dc0
          dc0=dc
          if(qabs(dc0).gt.deps) then
            dp2=(dgam**2)/dc0
          else
            dp2=dcm1*db(k)
          end if
          if(k.gt.1) dbeta(k)=ds*(dp2+db(k+1))
          ds=db(k+1)/(dp2+db(k+1))
          dc=dp2/(dp2+db(k+1))
          du=ds*(dgam+da(k+1)-dx)
          dalpha(k)=dgam+du+dx
   70   continue
        dbeta(1)=db(1)*(db(2)+(dx-da(1))**2)
        return
      else if(iopt.eq.8) then
        deh=0.q0
        dq=da(1)-dx
        do 80 k=1,n
          de=db(k+1)/dq
          dqh=2.q0*dx+dq+de-deh
          dbeta(k)=dqh*deh
          dq=da(k+1)-de-dx
          deh=dq*de/dqh
          dalpha(k)=-dx+dqh+deh
   80   continue
        dbeta(1)=db(1)*(db(2)+da(1)**2-dx**2)
        return
      else if(iopt.eq.9) then
        dalpha(1)=0.q0
        dalh=-dx+db(1)/dhr
        dqh=dalh+dx
        dq=-dx
        deh=0.q0
        do 90 k=2,n
          de=dqh+deh-2.q0*dx-dq
          deh=da(k-1)+dx-dqh
          dbeta(k)=dq*de
          dq=dqh*deh/de
          dalpha(k)=dq+de+dx
          dqh=db(k)/deh
   90   continue
        dbeta(1)=-dhr/dx 
        return
      else
        ierr=2
        return
      end if
      end

c
c
      subroutine qgchri(n,iopt,nu0,numax,deps,da,db,dx,dy,dalpha,
     *dbeta,nu,ierr,ierrc,dnu,drhor,drhoi,droldr,droldi,ds,ds0,
     *ds1,ds2) 
c
c This is a quadruple-precision version of the routine  gchri.
c
      real*16 deps,da(numax),db(numax),dx,dy,dalpha(n),
     *dbeta(n),dnu(*),drhor(*),drhoi(*),droldr(*),droldi(*),
     *ds(n),ds0(*),ds1(*),ds2(*)
c
c The arrays  dnu,drhor,drhoi,droldr,droldi,ds0,ds1,ds2  are assumed
c to have dimension  2*n.
c
      if(n.lt.1) then
        ierr=-1
        return
      end if
      ierr=0
      nd=2*n
      ndm1=nd-1
      if(iopt.eq.1) then
        call qknum(ndm1,nu0,numax,dx,dy,deps,da,db,drhor,drhoi,nu,
     *    ierr,droldr,droldi)
        do 10 k=1,nd
          dnu(k)=-drhor(k)
   10   continue
        call qcheb(n,da,db,dnu,dalpha,dbeta,ds,ierrc,ds0,ds1,ds2)
        return
      else if(iopt.eq.2) then
        dy=qabs(dy)
        call qknum(ndm1,nu0,numax,dx,dy,deps,da,db,drhor,drhoi,nu,
     *    ierr,droldr,droldi)
        do 20 k=1,nd
          dnu(k)=-drhoi(k)/dy
   20   continue
        call qcheb(n,da,db,dnu,dalpha,dbeta,ds,ierrc,ds0,ds1,ds2)
        return
      else
        ierr=1
        return
      end if
      end

c
c
      subroutine qkern(n,nu0,numax,dx,dy,deps,da,db,dkerr,dkeri,
     *  nu,ierr,droldr,droldi)
c
c This is a quadruple-precision version of the routine  kern.
c
      real*16 dx,dy,deps,da(numax),db(numax),dkerr(*),
     *  dkeri(*),droldr(*),droldi(*),dp0r,dp0i,dpr,dpi,dpm1r,
     *  dpm1i,dden,dt
c
c The arrays  dkerr,dkeri,droldr,droldi  are assumed to have
c dimension  n+1.
c
      call qknum(n,nu0,numax,dx,dy,deps,da,db,dkerr,dkeri,nu,ierr,
     *  droldr,droldi)
      if(ierr.ne.0) return
      dp0r=0.q0
      dp0i=0.q0
      dpr=1.q0
      dpi=0.q0
      do 10 k=1,n
        dpm1r=dp0r
        dpm1i=dp0i
        dp0r=dpr
        dp0i=dpi
        dpr=(dx-da(k))*dp0r-dy*dp0i-db(k)*dpm1r
        dpi=(dx-da(k))*dp0i+dy*dp0r-db(k)*dpm1i
        dden=dpr**2+dpi**2
        dt=(dkerr(k+1)*dpr+dkeri(k+1)*dpi)/dden
        dkeri(k+1)=(dkeri(k+1)*dpr-dkerr(k+1)*dpi)/dden
        dkerr(k+1)=dt
   10 continue
      return
      end

c
c
      subroutine qknum(n,nu0,numax,dx,dy,deps,da,db,drhor,drhoi,nu,
     *ierr,droldr,droldi)
c
c This is a quadruple-precision version of the routine  knum.
c
      real*16 dx,dy,deps,da(numax),db(numax),drhor(*),
     *drhoi(*),droldr(*),droldi(*),drr,dri,dden,dt
c
c The arrays  drhor,drhoi,droldr,droldi  are assumed to have
c dimension  n+1.
c
      ierr=0
      np1=n+1
      if(nu0.gt.numax) then
        ierr=nu0
        return
      end if
      if(nu0.lt.np1) nu0=np1
      nu=nu0-5
      do 10 k=1,np1
        drhor(k)=0.q0
        drhoi(k)=0.q0
   10 continue
   20 nu=nu+5
      if(nu.gt.numax) then
        ierr=numax
        goto 60
      end if
      do 30 k=1,np1
        droldr(k)=drhor(k)
        droldi(k)=drhoi(k)
   30 continue
      drr=0.q0
      dri=0.q0
      do 40 j=1,nu
        j1=nu-j+1
        dden=(dx-da(j1)-drr)**2+(dy-dri)**2
        drr=db(j1)*(dx-da(j1)-drr)/dden
        dri=-db(j1)*(dy-dri)/dden
        if(j1.le.np1) then
          drhor(j1)=drr
          drhoi(j1)=dri
        end if
   40 continue
      do 50 k=1,np1
        if((drhor(k)-droldr(k))**2+(drhoi(k)-droldi(k))**2.gt.
     *    (deps**2)*(drhor(k)**2+drhoi(k)**2)) goto 20
   50 continue
   60 if(n.eq.0) return
      do 70 k=2,np1
        dt=drhor(k)*drhor(k-1)-drhoi(k)*drhoi(k-1)
        drhoi(k)=drhor(k)*drhoi(k-1)+drhoi(k)*drhor(k-1)
        drhor(k)=dt
   70 continue
      return
      end