73 Numbers — <math.h>, <fenv.h>, <tgmath.h>
What to know first
Looking back
Chapter 50 said not to compare reals with ==, and chapter 8 said 0.1 is not exactly representable. Then how does a mathematical function tell you when it receives “input it cannot calculate”?
A. There are two paths. It gives NaN or an infinity as the return value, and at the same time leaves the reason in errno — EDOM for outside the domain, ERANGE when the result exceeds the representable range. But an implementation may choose to report through the floating-point exception flags (<fenv.h>) instead of errno, so to check portably you must be ready to look at both. In the field it is usually simpler to check the return value with isnan and isinf.
The need for this chapter, and its context
By the end of this chapter
The questions this chapter answers
- How do the functions of
math.hreport failure — the return value alone cannot say?
73.1 The properties of NaN and infinity
examples-en/ch73/math.c
#include <stdio.h>
#include <math.h>
#include <errno.h>
#include <string.h>
int main(void)
{
/* (1) NaN is not even equal to itself */
double nan_v = nan("");
printf("NaN == NaN : %s\n", nan_v == nan_v ? "true" : "false");
printf("isnan(NaN) : %s\n", isnan(nan_v) ? "true" : "false");
/* (2) infinity and division by zero (for reals this is not UB) */
double inf_v = 1.0 / 0.0;
printf("1.0/0.0 : %f (isinf: %d)\n", inf_v, isinf(inf_v));
printf("0.0/0.0 : %s\n", isnan(0.0 / 0.0) ? "NaN" : "a number");
/* (3) a call outside the domain reports through errno and NaN */
errno = 0;
double r = sqrt(-1.0);
printf("sqrt(-1) : %s, errno=%s\n",
isnan(r) ? "NaN" : "a number", errno == EDOM ? "EDOM" : "0");
/* (4) beyond the range: ERANGE and infinity */
errno = 0;
double big = exp(1000.0);
printf("exp(1000) : %s, errno=%s\n",
isinf(big) ? "inf" : "a number", errno == ERANGE ? "ERANGE" : "0");
/* (5) the classification function sorts values by kind */
double vals[] = { 1.0, 0.0, -0.0, inf_v, nan_v, 1e-320 };
const char *names[] = { "1.0", "0.0", "-0.0", "inf", "NaN", "1e-320" };
for (size_t i = 0; i < sizeof vals / sizeof vals[0]; i++) {
const char *kind = "?";
switch (fpclassify(vals[i])) {
case FP_NORMAL: kind = "normal"; break;
case FP_SUBNORMAL: kind = "subnormal"; break;
case FP_ZERO: kind = "0"; break;
case FP_INFINITE: kind = "infinite"; break;
case FP_NAN: kind = "NaN"; break;
}
printf(" %-7s -> %-9s (sign bit %d)\n", names[i], kind, signbit(vals[i]) != 0);
}
/* (6) 0.0 and -0.0 compare equal, but their signs differ */
printf("0.0 == -0.0 : %s\n", 0.0 == -0.0 ? "true" : "false");
return 0;
}
Output
NaN == NaN : false
isnan(NaN) : true
1.0/0.0 : inf (isinf: 1)
0.0/0.0 : NaN
sqrt(-1) : NaN, errno=EDOM
exp(1000) : inf, errno=ERANGE
1.0 -> normal (sign bit 0)
0.0 -> 0 (sign bit 0)
-0.0 -> 0 (sign bit 1)
inf -> infinite (sign bit 0)
NaN -> NaN (sign bit 0)
1e-320 -> subnormal (sign bit 0)
0.0 == -0.0 : true
Four things to point out in the output.
① NaN is not equal to itself. IEEE 754 settled it so. Hence the old idiom that if x != x is true then x is NaN, while the standard function is isnan(x). Because of this property, sorting an array containing NaN with qsort breaks the comparator’s total order and the result collapses (chapter 66).
② Dividing a real by zero is not outside the contract. Unlike integer division (chapter 28), in an IEEE 754 environment it yields an infinity or a NaN. But the same holds that the very fact of dividing by zero is usually a bug.
③ sqrt(-1) is EDOM, exp(1000) is ERANGE. The former is outside the domain, the latter a case where the result exceeded the representable range. If you mean to look at errno, set it to 0 just before the call (chapter 75).
④ 0.0 and −0.0 are equal under ==. But the sign bit differs, and 1/0.0 and 1/-0.0 are +∞ and −∞ respectively. If the sign must be distinguished, use signbit.
A common misconception. “Comparing reals is safe if you use an epsilon”
The epsilon comparison learned in chapter 50 is not omnipotent. Absolute error (fabs(a-b) < eps) becomes meaningless when the values are large — near 1e9, 1e-9 is not even representable — and relative error collapses near zero. The prescription in the field is settling a tolerance that fits the situation, not using a universal constant. And it is better to ask first whether it can be handled with integers or fixed point so that the comparison is not needed at all (chapter 8′s story of calculating money).Q. How do the functions of math.h report failure — the return value alone cannot say?
A. In three ways. Outside the domain (say sqrt(-1)) they return NaN and set errno to EDOM. Beyond the range (say exp(1000)) they return infinity and set ERANGE. And the floating-point exception flags of <fenv.h> are raised.
The trouble is that how far each of the three is honoured varies between implementations. So the practical idiom is to clear errno = 0 before the call and check immediately after (chapter 75). To inspect the value itself use isnan and isinf — they say what they mean, unlike tricks such as x != x.
73.2 Functions often got wrong
| function | what it does | trap |
|---|---|---|
pow(x, y) | raising to a power | used for an integer power it can be slow and inexact |
round, nearbyint | rounding | round goes away from zero, nearbyint follows the current mode |
floor, ceil, trunc | cutting to an integer | the direction differs for negatives |
fmod, remainder | the remainder | their sign rules differ from each other |
abs, fabs | absolute value | ★ abs is for integers. used on a real it truncates |
atan2(y, x) | angle | the argument order is y, x |
isnan, isinf | classification | they are macros — they cannot be used as function pointers |
Table 74.1
pow(x, 2) is widely used, but for an integer square x * x is faster and exact. The compiler often optimises it, but not always.
The mistake of using abs on a real is especially quiet. <stdlib.h>‘s abs takes an int, so abs(-1.5) turns −1.5 into 1. Today’s compilers warn, but it is easy to miss in a file that does not include <math.h>.
73.3 Rounding modes and floating-point exceptions — <fenv.h>
Floating-point operations have two pieces of hidden global state.
The rounding mode — the default is “to the nearest value, ties to even”. It can be changed with fesetround, and once changed every subsequent real operation is affected.
The exception flags — flags are raised when division by zero, overflow, inexactness and so on occur. They are read with fetestexcept and cleared with feclearexcept. They are finer than errno, but to use this facility #pragma STDC FENV_ACCESS ON must be turned on — and then the compiler refrains from reordering real operations, so optimisation is reduced.
Counter-example. Turning on -ffast-math and checking for NaN
cc -O2 -ffast-math app.c # tells the compiler "take it that NaN and infinity do not exist"if (isnan(x)) { /* this branch can vanish entirely */ }Options of the -ffast-math family tell the compiler it may assume associativity and ignore the existence of NaN and −0.0. Speed is gained, but the checking code can vanish under optimisation — the real-number edition of the “bug that appears only in release” seen in chapter 17. In a program where numerical accuracy matters, not turning it on is the default.
73.4 Type-generic — <tgmath.h>
sqrt is for double, sqrtf for float, sqrtl for long double. Include <tgmath.h> and the edition fitting the argument’s type is chosen by sqrt(x) alone — the representative case of the _Generic seen in chapter 58 being used in the standard library.
It is convenient but has a price. Being macros, they cannot be passed as function pointers, and there may be implementations that evaluate the argument twice, so putting in an expression with side effects is dangerous.
In practice. The same calculation, a different answer — the history of excess precision
x86′s old floating-point unit (x87) calculated internally in 80 bits. So it happened that the same double operation differed depending on whether it was still in a register or had been stored to memory — change the optimisation level and the result changed minutely, and x == y that had been true could become false.
C99 made this circumstance explicit with FLT_EVAL_METHOD, and today’s 64-bit x86 uses SSE so the problem has greatly diminished. But the possibility of “the same code, a different answer” still remains in compilation options and the target machine — the reason chapter 50 said “real-number calculation needs reproducibility looked after separately.”
Recap
Numbers in summary.
| situation | what to use | what to beware of |
|---|---|---|
| checking for NaN | isnan | x == NaN is always false |
| checking for infinity | isinf | dividing a real by zero is not UB |
| the kind of a value | fpclassify | the existence of subnormal numbers |
| domain and range errors | the return value + errno | errno = 0 just before the call |
| integer squares | x * x | pow(x, 2) |
| absolute value of a real | fabs | abs (for integers) |
| per-type functions | <tgmath.h> | macros — no arguments with side effects |
| fast-math options | off by default | checking code vanishes |
Table 74.2
We have passed numbers. The next chapter is time — a place with unusually much that the standard does not settle for you.