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