Exponent(inf) trouble

This program gives one or other of two different outputs when compiled with 5 compilers in a Linux x86_64 system:

program testexpinf ! file testexpinf.f90
  implicit none
  character:: cinf*3 = "INF"
  real     :: infvalue
  read(cinf,*)infvalue
  print "(/,A,F3.0)" ,'Reading "INF" gives ',infvalue
  print "(A,B32.32)" ,'transfer(inf,[1]) = ',transfer(infvalue,[1])
  print "(A,B32.32)" ,'exponent(inf)     = ',exponent(infvalue)
  print "(A,I0)"     ,'exponent(inf)     = ',exponent(infvalue)
end program testexpinf

Output from gfortran 15.1.0, lfortran 0.62.0, or ifort 2021.9.0:

Reading "INF" gives Inf
transfer(inf,[1]) = 01111111100000000000000000000000
exponent(inf)     = 01111111111111111111111111111111
exponent(inf)     = 2147483647

Output from g95 0.94 or AOCC flang 5.1.0:

Reading "INF" gives Inf
transfer(inf,[1]) = 01111111100000000000000000000000
exponent(inf)     = 00000000000000000000000010000001
exponent(inf)     = 129

It seems that either

  1. the program is not standard-conforming, or
  2. either 2 or 3 of those compilers have bugs, or
  3. the program is standard-conforming, but those diffent outputs are OK.

Which case applies?

The description of EXPONENT states

EXPONENT(3) returns the value of the exponent part of X

  • If X is zero the value returned is zero.
  • If X is an IEEE infinity or NaN, the result has the value HUGE(0).

The method of getting an INF by reading the string “INF” is perhaps
the most portable and reliable method of getting an Infinite value
but technically it might not be an IEEE infinity. I added a method
using the IEEE_ARiTHMETIC intrinsic module as well although it is
very likely you will get the same value. The B format descriptor can
be used directly on the value, a TRANSFER to an integer is not
required that I know of.

So by my reading you should be getting HUGE(0).

So if MINEXPONENT(0.0) is probably -125 then exponent(tiny(0.0) should be -125 so what is exponent(tiny(0.0)/2) supported to produce?

program testexpinf ! file testexpinf.f90
   implicit none
   character:: cinf*3 = "INF"
   real     :: infvalue
   read(cinf,*)infvalue
   print "(/,A,F3.0)" ,'Reading "INF" gives ',infvalue
   print "(A,b32.32)" ,'binary            = ',infvalue
   print "(A,B32.32)" ,'transfer(inf,[1]) = ',transfer(infvalue,[1])
   print "(A,B32.32)" ,'exponent(inf)     = ',exponent(infvalue)
   print "(A,I0)"     ,'exponent(inf)     = ',exponent(infvalue)
   call ieee()
contains
   subroutine ieee()
      use, intrinsic :: ieee_arithmetic
      real :: my_inf, my_neg_inf
      ! Generate positive infinity
      print *, 'IEEE'
      my_inf = ieee_value(my_inf, ieee_positive_inf)
      ! Generate negative infinity
      my_neg_inf = ieee_value(my_neg_inf, ieee_negative_inf)

      print "(A,b32.32)" ,'ieee infinity         = ',my_inf
      print "(A,B32.32)" ,'exponent(inf)         = ',exponent(my_inf)
      print "(A,I0)"     ,'exponent(inf)         = ',exponent(my_inf)

      print "(A,b32.32)" ,'negative ieee infinity= ',my_neg_inf
      print "(A,B32.32)" ,'exponent(neg_inf)         = ',exponent(my_neg_inf)
      print "(A,I0)"     ,'exponent(neg_inf)         = ',exponent(my_neg_inf)
   end subroutine ieee
end program testexpinf
NAME
  EXPONENT(3) ‐ [MODEL:COMPONENTS] Exponent of floating‐point number

SYNOPSIS
  result = exponent(x)

           elemental integer function exponent(x)

            real(kind=**),intent(in) :: x

CHARACTERISTICS
  •  X shall be of type real of any valid kind

  •  the result is a default integer type

DESCRIPTION
  EXPONENT(3) returns the value of the exponent part of X, provided the
  exponent is within the range of default integers.

OPTIONS
  •  X : the value to query the exponent of

RESULT
  EXPONENT(3) returns the value of the exponent part of X

  If X is zero the value returned is zero.

  If X is an IEEE infinity or NaN, the result has the value HUGE(0).

EXAMPLES
  Sample program:

      program demo_exponent
      implicit none
      real :: x = 1.0
      integer :: i
         i = exponent(x)
         print *, i
         print *, exponent(0.0)
         print *, exponent([10.0,100.0,1000.0,‐10000.0])
         ! beware of overflow, it may occur silently
         !print *, 2**[10.0,100.0,1000.0,‐10000.0]
         print *, exponent(huge(0.0))
         print *, exponent(tiny(0.0))
      end program demo_exponent

  Results:

       >            4           7          10          14
       >          128
       >         ‐125

STANDARD
  Fortran 95

SEE ALSO
  DIGITS(3),  EPSILON(3),  FRACTION(3),  HUGE(3),  MAXEXPONENT(3),
  MINEXPONENT(3),  NEAREST(3),  PRECISION(3), RADIX(3), RANGE(3),
  RRSPACING(3), SCALE(3), SET_EXPONENT(3), SPACING(3), TINY(3)

  Fortran intrinsic descriptions

Thank you @urbanjost. With the 5 compilers I tried, reading “INF” gave the same result every time as the IEEE positive infinity. The g95 bug is unlikely to be fixed but I hope the flang one will be. (Can’t the various flang compilers be given different names?)

program testexpinf ! file testexpinf.f90
   implicit none
   character(len=*),parameter :: gg='(*(g0))', ab='(a,*(b32.32,1x))'

   print gg, 'EXPONENT RANGE'
   print gg, 'minexponent(0.0)      = ',minexponent(0.0) 
   print gg, 'maxexponent(0.0)      = ',maxexponent(0.0) 
   print gg, 'TINY/HUGE'
   print gg, 'exponent(tiny(0.0))   = ',exponent(tiny(0.0))
   print gg, 'exponent(huge(0.0))   = ',exponent(huge(0.0))
   print gg, 'A VALUE OUTSIDE OF ABOVE RANGE IS SUSPECT'
   print gg, 'exponent(tiny(0.0)/2) = ',exponent(tiny(0.0)/2)

end program testexpinf

With gfortran the last line prints -126, but the previous lines indicate the range is -125 to 128.
It seems suspect if the returned value is not in the exponent range. So I think a few dusty corners remain.

I doubt this is covered by the fortran standard. In practice, it probably depends on whether gradual underflow is active and whether or not IEEE arithmetic is enabled. Your result of -126 seems reasonable with gradual underflow, but a return value of 0 would also seem likely with abrupt underflow.

I agree it is not as simple as i might appear at first glance, but given than minexporent(x) indicates the lowest exponent value is -125 on my platform I think it should say either only a pure zero is a special case and return -125 or say that this value is so close to zero it should be treated as such and return 0; but returning a value outside of the range seems like a bug; else minexponent(x) should return -126. So it probably is not explicitly covered and it did not produce an ICE but to me that answer is still a bug based on the aforesaid logic but it is interesting to consider other viewpoints, as I also cannot say there is anything wrong with your statements. So for now I think I will add all this to my Dusty_Corners/ collection! I keep a man-page of the core intrinsics and this seems like it should get at least a footnote in the associated example program. Might be useful in a Fortran QA suite too. I will at least add testing INF and NAN values as they are covered in the standard but not in my test demo program now that all this made me take a second look.

These discussions of subnormal numbers are always odd in some way or another within fortran because the fortran standard does not define them. Only normal floating point numbers are described in section 16.4, and for normal floating point numbers what you are saying holds. Namely, exponent() always returns an integer in the range (inclusive) emin to emax.

But people actually using floating point numbers have found that subnormal numbers are useful. The IEEE standard describes them in detail.

For IEEE 32-bit floating point numbers, the 8-bit exponent field stores a value between 0 and 255 (inclusive). To get the binary exponent of the floating point number you subtract 127. That is called a biased exponent. But the IEEE format uses the exponent field to flag special cases, so not all of those exponents are allowed for floating point numbers. In particular, the exponent field value of 0 is reserved for +0.0, -0.0, and for subnormal floating point values, and the field value of 255 is reserved for +INF, -INF, and for NaN values. So normal floating point numbers can only have exponent field values between 1 and 254 (inclusive). Subtracting the bias of 127 means that the floating point exponents range from -126 to +127. But, there is a hidden bit in the floating point representation, the leading bit which is always assumed to be 1, so that shifts the exponent range to -125 to +128, which is what your program prints. For each of those exponent values, there are 2^23 possible fractions, all corresponding to evenly spaced floating point values.

Then the question is what should happen with subnormal values? For these, the exponent field is 0 and there are 2^23 possible fractions. Each of these correspond to an exact real number value (actually a rational number value), and if subnormals are enabled, then arithmetic (addition, subtraction, multiplication, division) can be performed along with normal floating point values and with other subnormal values. However, they are limited in the sense that both underflows (to zero) and overflows can occur, so when used the programmer must be aware of their limitations.

That is why exponent() returns a value that is smaller than that for tiny(), because their exponent really is smaller. They are nonzero, so when compared to 0.0 they are not equal.

Thanks for that wonderful exposition. It cleared up other questions I had not even voiced.

I wrote a little program to demonstrate how subnormals behave.

program ieee
   use, intrinsic :: iso_fortran_env, wp=>real32
   use, intrinsic :: ieee_arithmetic
   implicit none
   integer :: i
   real(wp) :: x
   logical :: gradual_mode, gradual_tmp
   character(*), parameter :: cfmt = '(*(g0))'

   call ieee_get_underflow_mode( gradual_mode )
   call ieee_set_underflow_mode( gradual = .true. )
   call ieee_get_underflow_mode( gradual_tmp )

   print cfmt, 'original mode=', gradual_mode, ', new mode=', gradual_tmp

   x = tiny(x)  ! smallest normal number.
   do i = 1, 25
      call printx( x )
      x = x * 0.5
   enddo

contains

   subroutine printx( x )
      real(wp), intent(in) :: x
      integer :: ix, s, e, f
      character(*), parameter :: xfmt = '(e14.7,1x, b1.1,1x,b8.8,1x,b23.23,1x,l1,1x,i0)'
      ix = transfer(x,ix)
      s = ibits(ix,31,1)
      e = ibits(ix,23,8)
      f = ibits(ix,0,23)
      print xfmt, x, s, e, f, ieee_is_normal(x), exponent(x)
      return
   end subroutine printx

end program ieee

$ gfortran ieee.f90 && a.out
original mode=T, new mode=T
 0.1175494E-37 0 00000001 00000000000000000000000 T -125
 0.5877472E-38 0 00000000 10000000000000000000000 F -126
 0.2938736E-38 0 00000000 01000000000000000000000 F -127
 0.1469368E-38 0 00000000 00100000000000000000000 F -128
 0.7346840E-39 0 00000000 00010000000000000000000 F -129
 0.3673420E-39 0 00000000 00001000000000000000000 F -130
 0.1836710E-39 0 00000000 00000100000000000000000 F -131
 0.9183550E-40 0 00000000 00000010000000000000000 F -132
 0.4591775E-40 0 00000000 00000001000000000000000 F -133
 0.2295887E-40 0 00000000 00000000100000000000000 F -134
 0.1147944E-40 0 00000000 00000000010000000000000 F -135
 0.5739719E-41 0 00000000 00000000001000000000000 F -136
 0.2869859E-41 0 00000000 00000000000100000000000 F -137
 0.1434930E-41 0 00000000 00000000000010000000000 F -138
 0.7174648E-42 0 00000000 00000000000001000000000 F -139
 0.3587324E-42 0 00000000 00000000000000100000000 F -140
 0.1793662E-42 0 00000000 00000000000000010000000 F -141
 0.8968310E-43 0 00000000 00000000000000001000000 F -142
 0.4484155E-43 0 00000000 00000000000000000100000 F -143
 0.2242078E-43 0 00000000 00000000000000000010000 F -144
 0.1121039E-43 0 00000000 00000000000000000001000 F -145
 0.5605194E-44 0 00000000 00000000000000000000100 F -146
 0.2802597E-44 0 00000000 00000000000000000000010 F -147
 0.1401298E-44 0 00000000 00000000000000000000001 F -148
 0.0000000E+00 0 00000000 00000000000000000000000 T 0

However, the output surprised me a little. I was expecting the same exponent, -126, to be printed for all of the subnormal values and 0 exponent for the 0.0 value after underflow occurs. But that’s not what gets printed. Instead, the exponents that are printed scale with the exponent in the floating point number combined with the leading zeros in the subnormal fraction.

NaN produced a greater variety of results than Inf did. Finding NaN in 5 ways (reading ‘NAN’, 0/0, 0**0, IEEE quiet NaN and signaling NaN) gave me 4 different values with format Z8.8. IEEE quiet NaN was the same as 0/0 for ifx, but it was the same as the read-in NaN for gfortran and lfortran. Test program:

program nanieee
  use ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_signaling_nan
  implicit none
  real      :: x0, readnan, calcnan(2), ieeenan(2)
  character :: charnan*3 = 'NAN', fmt*8 = "(A,Z8.8)"
  x0 = 0 
  read(charnan,*) readnan
  calcnan = [x0/x0, x0**x0]
  ieeenan = ieee_value(x0,[ieee_quiet_nan, ieee_signaling_nan])
  print fmt,'read-in        NaN = ',readnan
  print fmt,'x0/x0          NaN = ',calcnan(1)
  print fmt,'x0**x0         NaN = ',calcnan(2)
  print fmt,'IEEE quiet     NaN = ',ieeenan(1)
  print fmt,'IEEE signaling NaN = ',ieeenan(2)
end program nanieee

Nice. If the platform is IEEE-conformant, x0**x0 would not be a NAN but 1; but the x0/x0 one is surprising . in particular. That ends up in gfortran that sign(1.0,x0/x0) is -1.0 but that ieee_is_
negative(x0/x0) is .FALSE., as the x0/x0 value has the sign bit set

There might be something but as far as I remember nothing would guarantee whether you got a signaling or quiet NaN except for the IEEE functions so it is probably OK if the NaN done in any other manner is someines a quiet NaN and sometimes a signaling NaN depending on the platform and compiler options used I think. Rereading a few parts of the standard with that in mind to see if that is left up to the platform or not. Interesting.

There are 2**23 possible fractions with the exponent field of 255. f==0 corresponds to +Inf for the sign bit 0 and to -Inf for the sign bit 1. All of the other ones will be identified as NaN values. I guess I’m unsure exactly how much of this is defined by IEEE and how much are just conventions used by the compilers (fortran, C, C++, etc.), but generally any floating point value with an exponent field of 255 that is not +Inf or -Inf will be identified as a NaN. You can test that with this little program.

program NaN
   use, intrinsic :: iso_fortran_env, wp=>real32
   use, intrinsic :: ieee_arithmetic
   implicit none
   integer :: i, s, e, f, knt
   character(*), parameter :: cfmt = '(*(g0))'

   s = 1   ! sign bit.
   e = 255 ! exponent field.
   knt = 0
   do i = 0, 2**23-1
      f = i
      call printx( s, e, f )
   enddo
   print cfmt, 'knt=', knt, ' 2**23=', 2**23
contains

   subroutine printx( s, e, f )
      ! assemble x from the sef fields and test for NaN.
      integer, intent(in) :: s, e, f
      real(wp) :: x
      integer :: ix
      character(*), parameter :: xfmt = '(e14.7,1x, b1.1,1x,b8.8,1x,b23.23,1x,l1,1x,i0)'
      ix = 0
      call mvbits( s, 0,  1, ix, 31 )
      call mvbits( e, 0,  8, ix, 23 )
      call mvbits( f, 0, 23, ix,  0 )
      x = transfer(ix,x)
      if ( .not. ieee_is_nan(x) ) then
         print xfmt, x, s, e, f, ieee_is_normal(x), exponent(x)
      else
         knt = knt + 1
      endif
      return
   end subroutine printx

end program NaN

$ flang NaN.f90 && a.out
          -Inf 1 11111111 00000000000000000000000 F 2147483647
knt=8388607 2**23=8388608

As you can see, all but one of the possible f values are identified as NaN, and that one exception is -Inf in this case. The same happens with s=0. The exponent() intrinsic prints huge() instead of 128 (=255-127) for the two Inf cases, but that is a special case required by the fortran standard. If x=0.0 or x=-0.0, then the exponent() function returns 0 instead of -126, so those are also special cases defined by the standard.

There are also set_exponent() and fraction() intrinsic functions that can be used to manipulate the fields within a floating point number. I did it with mvbits() in the above code instead just to be explicit. These other intrinsic functions are probably better to use in actual application codes.

Perhaps incidental but having initially been confused by special treatment of some values by intrinisics some caution is in order when using SIGN(), FRACTION(), and EXPONENT() to parse
values They work well enough to query normal numbers, but because of prescribed behaviors like FRACTION() being required to return a Nan if the input is a Nan or Inf means they do not work well when used to investigate the bits of Nan and Inf values, so I think the bit-level routines and transfer and BOZ -formatted internal I/O as used in everyone’s examples so far are more useful for the bit-level investigations being discussed.

“If the platform is IEEE-conformant, x0**x0 would not be a NAN but 1” raises a problem. If the platform is a Fortran compiler it has to calculate x0**x0 as exp(x0*log(x0)) = exp(0.0*(-Inf)) = exp(NaN) = NaN. That makes good sense mathematically because the limit of 0.0**a as a → 0 is 0.0 if a>0, infinity if a<0. There is a beautiful picture in Jahnke & Emde’s “Tables of Functions” (Dover 1945) on p1 (after p.306!) of the 3-dimensional surface z = x**y with x>0. You can find the 5 straight lines on it. I have not seen a better version published more recently.

In my test cases the signaling NaN was never one of the other 4 NaNs, but the quiet NaN was sometimes the same as the read-in NaN and sometimes the same as x0/x0. x0**x0 was never the same as another NaN.

Sorry I can’t try NAG Fortran. Would anyone like to?

0**0 is problematic both mathematically and as a floating-point operation. There are times even in mathematical series where it is treated as one. But if you just enter:

print *, x0**x0

I will be very surprised if you do not print “1” instead of “NaN”.
What result do you get?

Borrowing a few things from the previous examples …

program nanieee
   use, intrinsic :: iso_fortran_env, wp=>real32
   use ieee_arithmetic, only: ieee_value,   & 
    & ieee_quiet_nan,  ieee_signaling_nan,  &
    & ieee_is_nan,     ieee_is_normal,      &
    & ieee_is_finite,  ieee_is_negative     
   implicit none
   real(kind=wp,parameter     :: x0=0.0
   real(kind=wp)              :: readnan
   character                  :: charnan*3 = 'NAN'

   read(charnan,*) readnan

   call display('read-in         ',readnan )
   call display('x0/x0           ',x0/x0 )
   call display('x0**x0          ',x0**x0 )
   call display('one             ',1.0 )
   call display('IEEE quiet      ',ieee_value(0.0,ieee_quiet_nan) )
   call display('IEEE signaling  ',ieee_value(0.0,ieee_signaling_nan) )

contains

   subroutine display(title,value)
      character(len=*),intent(in) :: title
      real(kind=wp),intent(in)    :: value
      integer                     :: ix, s, e, f
      ix = transfer(value,ix)
      s = ibits(ix,31,1)
      e = ibits(ix,23,8)
      f = ibits(ix,0,23)
      write(*,'(a)',advance='no')                         title
      write(*,'(1x,f14.7)',advance='no')                  value
      write(*,'(1x,z8.8)',advance='no')                   value
      write(*,'(1x,b1.1,1x,b8.8,1x,b23.23)',advance='no') s, e, f 
      write(*,'(1x,a)',advance='no')                      merge('normal    ','not normal',ieee_is_normal(value))
      write(*,'(1x,a)',advance='no')                      merge('nan       ','not nan   ',ieee_is_nan(value))
      write(*,'(1x,a)',advance='no')                      merge('neg       ','not neg   ',ieee_is_negative(value))
      write(*,'(1x,a)',advance='no')                      merge('finite    ','not finite',ieee_is_finite(value))
      write(*,'(1x,g0)',advance='yes')                    exponent(value)
      !print *, int(sign(1.0,value)), merge('negative','        ',ieee_is_negative(value))
   end subroutine display

With default compiler switches gfortran produces “1” for 0.0**0.0:

end program nanieee
read-in                     NaN 7FC00000 0 11111111 10000000000000000000000 not normal nan        not neg    not finite 2147483647
x0/x0                       NaN FFC00000 1 11111111 10000000000000000000000 not normal nan        not neg    not finite 2147483647
x0**x0                1.0000000 3F800000 0 01111111 00000000000000000000000 normal     not nan    not neg    finite     1
one                   1.0000000 3F800000 0 01111111 00000000000000000000000 normal     not nan    not neg    finite     1
IEEE quiet                  NaN 7FC00000 0 11111111 10000000000000000000000 not normal nan        not neg    not finite 2147483647
IEEE signaling              NaN 7FA00000 0 11111111 01000000000000000000000 not normal nan        not neg    not finite 2147483647

I did get 1000000 or 1.00000000 depending on the compiler. I had failed to notice that in Z8.8 format 1.0 is 3F800000, which of course is not a NaN. Mea culpa. I have now changed my program to read ‘INF’ and get the IEEE +Inf, check whether they are the same, and that what is printed as Inf or NaN really is, and that Inf-Inf is NaN (replacing 0**0). My revised program is below. It works with lfortran, gfortran, ifx and AMD flang.

program infnan
  use ieee_arithmetic, only: ieee_value, ieee_positive_inf, &
       ieee_quiet_nan, ieee_signaling_nan
  implicit none
  real     :: x0, readinf, readnan, calcnan(2), ieeenan(2), ieeeinf
  character:: charnan*3 = 'NAN', charinf*3 = 'INF', fmt*15 = "(A,F3.0,A,Z8.8)"
  x0 = 0
  read(charnan,*) readnan
  read(charinf,*) readinf
  calcnan = [x0/x0, readinf-readinf]
  ieeenan = ieee_value(x0,[ieee_quiet_nan, ieee_signaling_nan])
  ieeeinf = ieee_value(x0, ieee_positive_inf)
  print fmt,'IEEE positive  ',ieeeinf   ,' = ',ieeeinf
  print fmt,'read-in        ',readinf   ,' = ',readinf,' '
  print fmt,'read-in        ',readnan   ,' = ',readnan
  print fmt,'0.0/0.0        ',calcnan(1),' = ',calcnan(1)
  print fmt,'inf-inf        ',calcnan(2),' = ',calcnan(2)
  print fmt,'IEEE quiet     ',ieeenan(1),' = ',ieeenan(1)
  print fmt,'IEEE signaling ',ieeenan(2),' = ',ieeenan(2)
end program infnan

Considering how problematic NaN and Inf have been to use with past compilers the fact you can do that with four readily available compilers shows in a small way how far things have progressed in the last few years.

Early on I used silent Nan for data masking and initialized everything to signaling Nan to provide an “Early Warning System” for bugs but discontinued that practice not because it was not useful but because I encountered portability problems over and over. Looks like it might be safe to resume that! Might have been for a while(?) but “once bitten twice shy.”!

Of course, trying gfortran --fast-math and the equivalent can cause problems, as that allows even a IEEE-compliant compiler to violate IEEE rules. In the past at least some compilers would eliminate tests like val==val, often used to test for a NaN; and so on as well.

I guess there are a lot of inconsistencies with whatever the return value is chosen to be, but there is also a good reason for a return vaue of 1.0. Namely

lim(x->0) x^x = lim(x->0) exp(x*log(x)) = lim(x->0) exp(log(x)/(1/x))

This last form has an Inf/Inf behavior, so l’hopital’s rule can be applied. Differentiating both the numerator and denominator gives

lim(x->0) exp( (1/x) / (-1/(x*x))) = lim(x->0) exp(-x) = exp(0) = 1

If a program needs different behavior from this formal limit, then it looks like the programmer must explicitly test for the special cases.

I get the same output as you for gfortran on MacOS, but I had to comment out the line with the x0/x0 computation. I would need to find a compiler option that allows that line to compile. But it makes me wonder if the results you see (namely the sign bit on the NaN) would change if that expression were evaluated somehow at run time rather than at compile time. As I said before, any number with that 255 exponent field and with a nonzero fraction appears to be flagged as a NaN, regardless of the sign bit, but I do not know which of those bit patterns are purposely generated and which ones are just incidental. The answer to that might be hardware dependent. This is not discussed, for example, at the wikipedia page. IEEE 754 - Wikipedia Does anyone have a llnk to the IEEE specification?

I had to buy a copy. The document is available in a print on demand format ($141 as a non member) and pdf format ($112 as a non member).