Fortran Are these bugs for the GNU Fortran, g++, and clang compilers?

[Phyty]

TL;DR

I tried to calculate the smallest number in my machine with gfortran, g++ and clang. I did it by looping as:
a = a / 1.73
In gfortran after 1358 loops, I found that a = 4.9406564584124654E-324, and the value of a never changes from loop 1358 till the end 1500 loop.
It is questionable for me that a non zero real number here a = 4.9406564584124654E-324 divided by a non zero number here 1.73 will keep the same value.

I tried to calculate the smallest number in my machine with gfortran, g++ and clang. I did it by looping as:

a = a / 1.73

In gfortran after 1358 loops, I found that a = 4.9406564584124654E-324, and the value of a never changes from loop 1358 till the end 1500 loop.

It is questionable for me that a non zero real number here a = 4.9406564584124654E-324 divided by a non zero number here 1.73 will keep the same value. The phenomenon seems to the qualification of the ground state energy in a quantum system.

Is this bugs for the gfortran compiler?

I have tested similar codes with gfortran (GNU Fortran 8.2.0 2018, 5.40 201609), g++ (GNU 5.4.0 201609, clang-1000.10.44.4) on Linux (Ubuntu, KylinOS,) and Macintosh. I also tried a = a / 1.7.

All have similar results. A ground state qualification and unchanged value of a.

I guess that it is related to the binary system used by the computer system.

Here with my codes.

 

[Vanadium 50]

This is not a bug:

https://en.wikipedia.org/wiki/Floating-point_arithmetic

 

[DrClaude]

And let me also point out in particular:
https://en.wikipedia.org/wiki/Machine_epsilon

 

[sysprog]

From: https://www.math.utah.edu/~beebe/software/ieee/fpinfo.f $-$

      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