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