* (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