*     (Qtest)
*     This is an accuracy benchmark proposed by W. Kahan (UC Berkeley)
*     in ``Lecture Notes on the Status of IEEE 754'', May 23, 1995.
*     [24-May-1995]
      integer j, N
      parameter (N = 15)
      double precision e, qtrial, r(N), t
      external qdrtc

*                              for 24 sig. bits
      r(1) = 2**12 + 2.0d0

*                              for 6 hex
      r(2) = 2**12 + 2.5d0

*                              6 hex IBM
      r(3) = 16**3 + 1 + 1.0d0/16**2

*                              48 bits Cray
      r(4) = 2**24 + 2.0d0

*                              48 bits Cray rounded
      r(5) = 2**24 + 2.25d0

*                              48 bits Cray chopped
      r(6) = 2**24 + 3.0d0

*                              53 sig bits
      r(7) = 94906267.0d0

*                              53 sig bits
      r(8) = 94906267 + 0.25d0

*                              PowerPC, i860
      r(9) = 2**28 - 5.5d0

*                              PowerPC, i860
      r(10) = 2**28 - 4.5d0

*                              56 sig bits
      r(11) = 2**28 + 2.0d0

*                              and 14 hex
      r(12) = 2**28 + 2.25d0

*                              14 hex IBM
      r(13) = 16**7 + 1 + 1.0d0/16**6

*                              64 sig bits
      r(14) = 2.0d0**32 + 2.0d0

*                              64 sig bits
      r(15) = 2.0d0**32 + 2.25d0

*                              e = +Infinity
      e = 0.0d0
      e = 1.0d0/e

      do 10 j = 1, N
*                              Could be NaN
           t = qtrial(qdrtc, r(j))

           if ((t .lt. e) .or. .not.(t .eq. t)) e = t
   10 continue

      print *,'Worst accuracy is ', e, ' sig. bits'
      end

      double precision function log2(x)
      double precision x
      log2 = dlog(dabs(x))/dlog(2.0d0)
      end

      double precision function qtrial(qdrtc, r)
      double precision e1, e2, log2, r, p, q, x1, x2
      external qdrtc

      p = r - 2
      q = r - 1
      qtrial = 0

      print *, 'qdrtc() for r = ', r
      if (p .le. 0) then
          print *,'qtrial(...,r) expects r > 2'
      else if (.not.(((r - q) .eq. 1) .and. ((q - p) .eq. 1))) then
          print *, 'r is too big for qtrial(...,r)'
      else
          call qdrtc(p, q, r, x1, x2)
*                                   Could be NaN
          e1 = -log2(x1 - 1)
*                                   Heed parentheses!
          e2 = -log2((x2 - 1) - 2/p)
*                                   Min(NaN,NaN) is NaN
          qtrial = min(e1,e2)
          print *,'    gets ', e1, ' and ', e2, ' sig. bits'
          if (.not.(x1 .ge. 1.0d0))
     &        print *,'    and root ', x1, ' isn''t at least 1.'
      endif

      end

      subroutine qdrtc(p, q, r, x1, x2)
      double precision p, q, r, s, x1, x2
*                                   NaN if sqrt( < 0)
      s = dsqrt(q*q - p*r)
      if (q .gt. 0) then
           s = q + s
      else
           s = q - s
      endif
*                                   NaNs if not real
      x1 = r/s
      x2 = s/p

      end