Given that PI/2 can never be represented precisely in floating point, is it safe to assume that cos(a) can never return exact zero?
If this is the case, then the following pseudo-code will never enter the block (and it could be safely removed):
...
y = h/cos(a);
if (!isfinite(a))
{
// handle infinite y
}
Other than zero, the double precision value that comes closest to an exact multiple of π/2 is 6381956970095103 * 2^797, which is equal to:
(an odd integer) * π/2 + 2.983942503748063...e−19
Thus, for all double-precision values x, we have the bound:
|cos(x)| >= cos(2.983942503748063...e−19)
Note that this is a bound on the mathematically exact value, not on the value returned by the library function cos
. On a platform with a good-quality math library, this bound is sufficiently good that we can say that cos(x)
is not zero for any double-precision x
. In fact, it turns out that this is not unique to double; this property holds for all IEEE-754 basic types, if cos
is faithfully rounded.
However, that's not to say that this could never occur on a platform that had a spectacularly poor implementation of trigonometric argument reduction.
Even more importantly, it's critical to note that in your example y
can be infinite without cos(a)
being zero:
#include <math.h>
#include <stdio.h>
int main(int argc, char *argv[]) {
double a = 0x1.6ac5b262ca1ffp+849;
double h = 0x1.0p1022;
printf("cos(a) = %g\n", cos(a));
printf("h/cos(a) = %g\n", h/cos(a));
return 0;
}
compile and run:
scanon$ clang example.c && ./a.out
cos(a) = -4.68717e-19
h/cos(a) = -inf
Zero is one of several values that can be represented exactly. Many systems have a lookup table for common values of sin and cos, so it's not inconceivable that exactly zero could be returned.
But you are safer using a delta compare, before performing the division:
if (Abs(cos(a)) < 0.0000001)
{
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With