program qchalleps
c
c This evaluates the SIAM 100-digit challenge
c integral using the epsilon algorithm
c
      parameter(n=43,mcap1=41,mcap2=28)
      integer iarray(3)
      real*16 qeps,qpi,s(n),sa((n+1)/2),a0(mcap1),
     *b0(mcap1),a(mcap2),b(mcap2),u(mcap1),x(mcap1),
     *w(mcap1),e0(mcap1),quofx,sum,a1,sgn, qk,
     *E(n,n+1)
      write(*,1)
    1 format(/)
      call itime(iarray)
      write(*,2) iarray
    2 format(1x,'  system time ='3i3)
      qeps=.5q-33
      qpi=4.q0*qatan(1.q0)
      call qrecur(mcap1,2,0.q0,0.q0,a0,b0,ierr0)
      call qrecur(mcap2,1,0.q0,0.q0,a,b,ierr)
      call qgauss(mcap1,a0,b0,qeps,x,w,ierr,e0)
      do 10 i=1,mcap1
        x(i)=qpi*x(i)/2.q0
        w(i)=qpi*w(i)/2.q0
   10 continue
      do 20 i=1,mcap1
        u(i)=quofx(x(i))
        write(*,11) u(i)
   11   format(1x,d40.32)
   20 continue
      stop
      sum=0.q0
      do 30 i=1,mcap1
        sum=sum+w(i)*qcos(x(i))/(x(i)+u(i))
   30 continue
      s(1)=sum
c      write(*,12) sum
c   12 format(1x,d40.32)
c      stop
      sgn=1.q0
      do 70 k=1,n
        qk=qreal(k)
        if(k.ge.2) s(k)=s(k-1)+sgn*a1
        call qgauss(mcap2,a,b,qeps,x,w,ierr,e0)
        do 40 i=1,mcap2
          x(i)=qpi*x(i)/2.q0
          w(i)=qpi*w(i)/2.q0
   40   continue
        do 50 i=1,mcap2
          u(i)=quofx(x(i)+qk*qpi)
   50   continue
        sum=0.q0
        do 60 i=1,mcap2
          sum=sum+w(i)*qcos(x(i))/(x(i)+u(i)+
     *      qk*qpi)
   60   continue
        a1=sum
        sgn=-sgn
   70 continue
      call qepsalg(n,s,E)
      call itime(iarray)
      do 80 i=1,(n+1)/2
        sa(i)=E(n-2*i+2,2*i)
        write(*,3) sa(i)
    3   format(1x,d41.33)
   80 continue
      write(*,1)
      write(*,2) iarray
      stop
      end