CODEWITHSUNDEEP
c++floating-point
mas
mas

Reputation: 1255

Given x>0, is it possible for sqrt(x) = 0 to give a floating point error?

Question in title.

I have a section of code:

double ccss = c * c + s * s;
double sqrtCCSS = sqrt(ccss);
if (sqrtCCSS != 0)
{
    n = n1 / sqrtCCSS;
}

and am just wondering if this is safe:

double ccss = c * c + s * s;
if (ccss != 0)
{
    n = n1 / sqrt(ccss);
}

My gut tells me yes, but floating point error is still somewhat mysterious.

Update:

At least for Python so far this seems impossible:

[ins] In [4]: x = np.nextafter(0, 1)
         ...: print(x)
         ...: print(np.sqrt(x))
5e-324
2.2227587494850775e-162

If the smallest possible float in numpy (given numpy is mostly in C, this seems relevant), cannot result in zero, the only possible number that would result in zero would have to be something that introduces grave floating point error but is not the smallest floating point number.

Update: I changed the variables to satisfy comments

Upvotes: 14

Views: 678

Answers (1)

phuclv
phuclv

Reputation: 41952

and am just wondering if this is safe:

double xxyy = x * x + y * y;
if (xxyy != 0)
{
    n = n1 / sqrt(xxyy);
}

It's always safe because floating-point math doesn't trap and you can freely divide by zero or even Inf/NaN, unless you tell the compiler to trap on those cases

Anyway if you just want to check whether the denominator is zero then the question on the title is actually different from what you're talking about in the question's content. The answer to the question

Given x>0, is it possible for sqrt(x) to be zero?

is yes, if denormals-are-zero (DAZ) and/or flush-to-zero (FTZ) are turned on. In fact most non-x86 architectures have DAZ/FTZ enabled by default, or don't even support denormals and always turn DAZ/FTZ on, because operation on denormal numbers are very slow (see Why does changing 0.1f to 0 slow down performance by 10x?). In x86 you can also turn those flags on and some compilers will turn them on by default

For example the below sample code

f = 0x0Cp-1022;
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
printf("%a %a\n", f, sqrt(f));

may print 0x0.000000000000cp-1022 0x0p+0 because f is treated as zero by sqrt due to DAZ

Some demo on Godbolt. And gcc demo on ideone

The above snippet can also print 0x0p+0 0x0p+0 because f is treated as zero by printf due to FTZ even though it contains a non-zero denormalized value

If you set DAZ/FTZ before assigning f like this

_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
f = 0x0Cp-1022;
printf("%a %a\n", f, sqrt(f));

then it'll likely also print 0x0p+0 0x0p+0, but in this case f can be assigned an actual zero value instead, and printing its bit pattern will result in all zero bits

However if you've asked this

Given x and y which aren't both zero at the same time, is it possible for sqrt(x*x + y*y) to be zero

then the answer to that question is yes even if DAZ/FTZ aren't enabled, because the numbers can be so small that their product became underflow. Try something like this and see

double x = 0x1p-1021;
double y = 0x2p-1020; // or y = 0
double xxyy = x * x + y * y;
printf("%a %a\n", x, y);
printf("%a %a\n", xxyy, sqrt(xxyy));

Upvotes: 10

Related Questions