program fpinfo ************************************************************************ * (Floating-point information) * Fortran program to report the size of single-, double-, and * quadruple-precision datatypes, and the floating-point precision * in bits, the corresponding machine epsilon and smallest * floating-point numbers for each. * (29-Nov-2001) ************************************************************************ call ftest call dtest * call qtest end subroutine ftest real alog2, fstore real x * real*16 fppow2 double precision fppow2 character*36 note write (6,'(''single precision:'')') x = 1.0e+00 10 if (fstore(1.0e+00 + x/2.0e+00) .ne. 1.0e+00) then x = x / 2.0e+00 go to 10 end if write (6,'(8x,''sizeof(REAL) = '',i2)') 4 write (6,'(8x,''FLT_MANT_DIG = '',i2)') x (1 + nint(-alog2(x))) if (x .eq. real(fppow2(-23))) then note = '[IEEE 754 32-bit macheps]' else note = '[not IEEE 754 conformant]' end if write (6,'(8x,''machine epsilon = '', 1pe14.5e4, 1x, a)') x x, note x = 1.0e+00 20 if (fstore(x/2.0e+00) .ne. 0.0e+00) then x = x / 2.0e+00 go to 20 end if if (x .eq. real(fppow2(-149))) then note = '[IEEE 754 smallest 32-bit subnormal]' else if (x .eq. real(fppow2(-126))) then note = '[IEEE 754 smallest 32-bit normal]' else note = '[not IEEE 754 conformant]' end if write (6,'(8x,''smallest positive number = '', 1pe14.5e4, 1x, a)') x x, note write (6,'()') end subroutine dtest double precision dlog2, dstore double precision x * real*16 fppow2 double precision fppow2 character*36 note write (6,'(''double precision:'')') x = 1.0d+00 10 if (dstore(1.0d+00 + x/2.0d+00) .ne. 1.0d+00) then x = x / 2.0d+00 go to 10 end if write (6,'(8x,''sizeof(DOUBLE PRECISION) = '',i2)') 8 write (6,'(8x,''DBL_MANT_DIG = '',i2)') x (1 + nint(-dlog2(x))) if (x .eq. dble(fppow2(-52))) then note = '[IEEE 754 64-bit macheps]' else note = '[not IEEE 754 conformant]' end if write (6,'(8x,''machine epsilon = '', 1pe14.5e4, 1x, a)') x x, note x = 1.0d+00 20 if (dstore(x/2.0d+00) .ne. 0.0d+00) then x = x / 2.0d+00 go to 20 end if if (x .eq. dble(fppow2(-1074))) then note = '[IEEE 754 smallest 64-bit subnormal]' else if (x .eq. dble(fppow2(-1022))) then note = '[IEEE 754 smallest 64-bit normal]' else note = '[not IEEE 754 conformant]' end if write (6,'(8x,''smallest positive number = '', 1pe14.5e4, 1x, a)') x x, note write (6,'()') end * subroutine qtest * real*16 qlog2, qstore * real*16 x * real*16 fppow2 * character*42 note * * write (6,'(''quadruple precision:'')') * * x = 1.0q+00 * 10 if (qstore(1.0q+00 + x/2.0q+00) .ne. 1.0q+00) then * x = x / 2.0q+00 * go to 10 * end if * * write (6,'(8x,''sizeof(REAL*16) = '',i2)') 16 * * write (6,'(8x,''LDBL_MANT_DIG = '',i3)') * x (1 + nint(-qlog2(x))) * * if (x .eq. fppow2(-52)) then * note = '[IEEE 754 64-bit macheps]' * else if (x .eq. fppow2(-63)) then * note = '[IEEE 754 80-bit macheps]' * else if (x .eq. fppow2(-112)) then * note = '[IEEE 754 128-bit macheps]' * else * note = '[not IEEE 754 conformant]' * end if * write (6,'(8x,''machine epsilon = '', 1pe14.5e4, 1x, a)') * x x, note * * x = 1.0q+00 * 20 if (qstore(x/2.0q+00) .ne. 0.0q+00) then * x = x / 2.0q+00 * go to 20 * end if * if (x .eq. fppow2(-1074)) then * note = '[IEEE 754 smallest 64-bit subnormal]' * else if (x .eq. fppow2(-1022)) then * note = '[IEEE 754 smallest 64-bit normal]' * else if (x .eq. fppow2(-16446)) then * note = '[IEEE 754 smallest 80-bit subnormal]' * else if (x .eq. fppow2(-16382)) then * note = '[IEEE 754 smallest 80- and 128-bit normal]' * else if (x .eq. fppow2(-16494)) then * note = '[IEEE 754 smallest 128-bit subnormal]' * else if (x .eq. fppow2(-16382)) then * note = '[IEEE 754 smallest 128-bit normal]' * else * note = '[not IEEE 754 conformant]' * end if * write (6,'(8x,''smallest positive number = '', 1pe14.5e4, 1x, a)') * x x, note * * write (6,'()') * * end real function fstore(x) real x fstore = x end double precision function dstore(x) double precision x dstore = x end * real*16 function qstore(x) * real*16 x * qstore = x * end * real*16 function fppow2(n) * integer n * integer p * real*16 x * if (n .lt. 0) then * p = -n * x = 1.0q+00 / 2.0q+00 * else * p = n * x = 2.0q+00 * end if * fppow2 = 1.0q+00 * 10 if (p .gt. 0) then * p = p - 1 * fppow2 = fppow2 * x * go to 10 * end if * end double precision function fppow2(n) integer n integer p double precision x if (n .lt. 0) then p = -n x = 1.0d+00 / 2.0d+00 else p = n x = 2.0d+00 end if fppow2 = 1.0d+00 10 if (p .gt. 0) then p = p - 1 fppow2 = fppow2 * x go to 10 end if end real function alog2(x) real x alog2 = alog(x)/alog(2.0e+00) end double precision function dlog2(x) double precision x dlog2 = dlog(x)/dlog(2.0d+00) end * real*16 function qlog2(x) * real*16 x * qlog2 = qlog(x)/qlog(2.0q+00) * end