When does underflow occur?

Question

I get into a situation where calculating 1.77e-308/10 triggers an underflow exception, but calculating 1.777e-308/10 does not. This is strange because:

Underflow occurs when the true result of a floating point operation is smaller in magnitude (that is, closer to zero) than the smallest value representable as a normal floating point number in the target datatype (from Arithmetic Underflow, Wikipedia)

In other words, if we calculate x/y where both x and y are double, then underflow should occur if 0 < |x/y| < 2.2251e-308 (the smallest positive normalized double is 2.2251e-308). In theory, therefore, both 1.77e-308/10 and 1.777e-308/10 should trigger an underflow exception. The theory contradicts with what I have tested with the C program below.

#include <stdio.h>
#include <fenv.h>
#include <math.h>


int main(){
  double x,y;

  // x = 1.77e-308 => underflow
  // x = 1.777e-308 gives  ==> no underflow
  x=1.77e-308;

  feclearexcept(FE_ALL_EXCEPT);
  y=x/10.0;
  if (fetestexcept(FE_UNDERFLOW)) {
    puts("Underflow\n");
  }
  else puts("No underflow\n");
}

To compile the program, I used gcc program.c -lm; I also tried Clang, which gave me the same result. Any explanation?

[Edits] I have shared the code above via this online IDE.

Cacciatore

Can you show y value? – Masuria

How did you determine the smallest normalized double on your machine? – Davinadavine

On my platform is the opposite: 1.77e-308 trigger an underflow while 1.777e-308;` doesn't. g++ (Debian 4.9.2-10) 4.9.2 – Barker

@Davinadavine I determined the smallest normalized double via std::numeric_limits<double>::min() (with a separate c++ program). – Cacciatore

@Masuria I just tested 'y' value. It was as expected,namely y = 1.77e-309 or y = 1.777e-309 if x=1.77e-308 or x = 1.777e-308, respectively. I used "%g" of printf to print the 'y' value. – Cacciatore

@Barker Thanks. Can you try gcc or clang on your platform? – Cacciatore

Just for kicks, maybe check the min in C instead of in C++? en.cppreference.com/w/c/types/limits – Davinadavine

@Brick. Thanks. Just tried. It gives 2.2251e-308. I used printf ("%.5g", DBL_MIN). – Cacciatore

@Cacciatore I tested it with gcc also: still the same on my platform. And DBL_MIN is 2.2251e-308. – Barker

The text in the question and the comments in the code don't match as to which gives the underflow and which doesn't. I suspect that the code comments are correct, and that would match what @Barker is reporting. – Davinadavine

@Brick. You are right. I edited the text. I also shared the code via an online IDE so you can try. – Cacciatore

Answer 1

Underflow is not only a question of range, but also of precision/rounding.

7.12.1 Treatment of error conditions
The result underflows if the magnitude of the mathematical result is so small that the mathematical result cannot be represented, without extraordinary roundoff error, in an object of the specified type. C11 §7.12.1 6

1.777e-308, converted to the nearest binary64 0x1.98e566222bcfcp-1023, happens to have a significand (0x198E566222BCFC, 7193376082541820) that is a multiple of 10. So dividing by 10 is exact. No roundoff error.

I find this easier to demo with hex notation. Note that dividing by 2 is always exact, except for the smallest value.

#include <float.h>
#include <stdio.h>
#include <fenv.h>
#include <math.h>

int uf_test(double x, double denominator){
  printf("%.17e %24a ", x, x);
  feclearexcept(FE_ALL_EXCEPT);
  double y=x/denominator;
  int uf = !!fetestexcept(FE_UNDERFLOW);
  printf("%-24a %s\n", y, uf ? "Underflow" : "");
  return uf;
}

int main(void) {
  uf_test(DBL_MIN, 2.0);
  uf_test(1.777e-308, 2.0);
  uf_test(1.77e-308, 2.0);
  uf_test(DBL_TRUE_MIN, 2.0);

  uf_test(pow(2.0, -1000), 10.0);
  uf_test(DBL_MIN, 10.0);
  uf_test(1.777e-308, 10.0);
  uf_test(1.77e-308, 10.0);
  uf_test(DBL_TRUE_MIN, 10.0);
  return 0;
}

Output

2.22507385850720138e-308                0x1p-1022 0x1p-1023                
1.77700000000000015e-308  0x1.98e566222bcfcp-1023 0x1.98e566222bcfcp-1024  
1.77000000000000003e-308  0x1.97490d21e478cp-1023 0x1.97490d21e478cp-1024  
4.94065645841246544e-324                0x1p-1074 0x0p+0                   Underflow

// No underflow as inexact result is not too small
9.33263618503218879e-302                0x1p-1000 0x1.999999999999ap-1004  
// Underflow as result is too small and inexact
2.22507385850720138e-308                0x1p-1022 0x1.99999999999ap-1026   Underflow
// No underflow as result is exact
1.77700000000000015e-308  0x1.98e566222bcfcp-1023 0x1.471deb4e8973p-1026   
1.77000000000000003e-308  0x1.97490d21e478cp-1023 0x1.45d40a818394p-1026   Underflow
4.94065645841246544e-324                0x1p-1074 0x0p+0                   Underflow

Nieshanieto

Excellent explanation! Thank you. – Cacciatore

Maybe I am missing something obvious, but why the double negation in front of fetestexcept()? – Ania

@DavidBowling fetestexcept() returns a int: "the value of the bitwise OR". From uf_test() i only wanted to return 0 or 1, so used !!, but then did not use that result in this posted code. Better code would have used bool uf = fetestexcept(FE_UNDERFLOW); to achieve the same purpose. – Nieshanieto

Answer 2

Checking the documentation for the function that you called, leads to the definition:

FE_UNDERFLOW the result of an earlier floating-point operation was subnormal with a loss of precision

http://en.cppreference.com/w/c/numeric/fenv/FE_exceptions

You've verified that your number is subnormal, I think. The test also includes loss of precision. If you print more significant figures, you'll find that the one reporting the overflow does seem to lose precision at about 16 decimal places. I'm not clear how many significant figure to expect on a subnormal number, but I think this must be your answer.

Davinadavine