c - 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
");
  }
  else puts("No underflow
");
}

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.

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