When does underflow occur?

问题:

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.

回答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

回答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.