Given x>0, is it possible for sqrt(x) = 0 given floating point error?-CodePudding

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

CodePudding user response：

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 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) = 0

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

Some demo on Godbolt. And gcc demo on ideone

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) = 0

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));