C++ Sqrt Function Precision for Full Squares

C++ sqrt function precision for full squares

No, you cannot be guaranteed. For integers and their squares that fit in the dynamic range of the floating point type's mantissa (2^53 for a typical C/C++ double), you're likely to be OK, but not necessarily guaranteed.

You should avoid equals comparisons between floating point values and exact values, especially exact integer values. Floating point rounding modes and other such things can really get in your way.

You either want to use a "comparison range" to accept an "approximately equal" result, or recast your algorithm in terms of integers. There are multiple StackOverflow questions covering floating point equality comparisons. I suggest you search for them and read up.

For a certain class of problem, I wrote up an alternate solution here:
Find n-th root of all numbers within an interval

That solution took a different approach than relying on tricky floating point arithmetic.

Guaranteed precision of sqrt function in C/C++

For C99, there are no specific requirements. But most implementations try to support Annex F: IEC 60559 floating-point arithmetic as good as possible. It says:

An implementation that defines __STDC_IEC_559__ shall conform to the specifications in this annex.

And:

The sqrt functions in <math.h> provide the IEC 60559 square root operation.

IEC 60559 (equivalent to IEEE 754) says about basic operations like sqrt:

Except for binary <-> decimal conversion, each of the operations shall be performed as if it first produced an intermediate result correct to infinite precision and with unbounded range, and then coerced this intermediate result to fit in the destination's format.

The final step consists of rounding according to several rounding modes but the result must always be the closest representable value in the target precision.

sqrt, perfect squares and floating point errors

In IEEE 754 floating-point, if the double-precision value x is the square of a nonnegative representable number y (i.e. y*y == x and the computation of y*y does not involve any rounding, overflow, or underflow), then sqrt(x) will return y.

This is all because sqrt is required to be correctly-rounded by the IEEE 754 standard. That is, sqrt(x), for any x, will be the closest double to the actual square root of x. That sqrt works for perfect squares is a simple corollary of this fact.

If you want to check whether a floating-point number is a perfect square, here's the simplest code I can think of:

int issquare(double d) {
  if (signbit(d)) return false;
  feclearexcept(FE_INEXACT);
  double dd = sqrt(d);
  asm volatile("" : "+x"(dd));
  return !fetestexcept(FE_INEXACT);
}

I need the empty asm volatile block that depends on dd because otherwise your compiler might be clever and "optimise" away the calculation of dd.

I used a couple of weird functions from fenv.h, namely feclearexcept and fetestexcept. It's probably a good idea to look at their man pages.

Another strategy that you might be able to make work is to compute the square root, check whether it has set bits in the low 26 bits of the mantissa, and complain if it does. I try this approach below.

And I needed to check whether d is zero because otherwise it can return true for -0.0.

EDIT: Eric Postpischil suggested that hacking around with the mantissa might be better. Given that the above issquare doesn't work in another popular compiler, clang, I tend to agree. I think the following code works:

int _issquare2(double d) {
  if (signbit(d)) return 0;
  int foo;
  double s = sqrt(d);
  double a = frexp(s, &foo);
  frexp(d, &foo);
  if (foo & 1) {
    return (a + 33554432.0) - 33554432.0 == a && s*s == d;
  } else {
    return (a + 67108864.0) - 67108864.0 == a;
  }
}

Adding and subtracting 67108864.0 from a has the effect of wiping the low 26 bits of the mantissa. We will get a back exactly when those bits were clear in the first place.

sqrt function in C not returning exact value

The function distance returns an integer, therefore it truncates the result of the sqrt. You should change the declaration of:

 int distance(int x1, int y1, int x2, int y2)

to

double distance(int x1, int y1, int x2, int y2)

IEEE double such that sqrt(x*x) ≠ x

Sylvie Boldo has formally proved that a floating-point number satisfying the conditions in your question does not exist.

Quoting the abstract of the article:

Floating-point experts know that mathematical formulas may fail or
give imprecise results when implemented in floating-point arithmetic.
This article describes an example where, surprisingly, it is
absolutely not the case. Indeed, using radix 2 and an unbounded
exponent range, the computation of the square root of the square of a
floating-point number a is exactly |a|. A consequence is the fact that
the floating-point computation of a/ sqrt (a² + b²) is always in the
interval [−1, 1]. This removes the need for a test when calling an
arccos or an arcsin on this value. For more guarantees, this property
was formally checked using the Coq proof assistant and the Flocq
library. The conclusion will give hints on what happens without
assumptions and in other radices, where the behavior is very
different.

“using radix 2” was likely implicit in your question (although the IEEE has also standardized decimal floating-point number formats and operations), and “an unbounded exponent range” is equivalent to your “no overflow or underflow” restriction.

A reason making the property possible at all is that x*x “expands” (the interval [1,2] is mapped to [1,4], for instance) in a way such that, when there is no overflow or underflow, the rounding that can happen for * is benign and x is still the closest representable floating-point number to the real square root of the floating-point product x*x. This hand-wavy argument does not constitute a proof, so it's a good thing that the article linked above contains one.

Is there a value of type `double`, `K`, such that `K * K == 3.0`?

Testing with Python is valid I think, since both use the IEEE-754 representation for doubles along with the rules for operations on same.

The closest possible double to the square root of 3 is slightly low.

>>> Sqrt3 = 3**0.5
>>> Sqrt3*Sqrt3
2.9999999999999996

The next available value is too high.

>>> import numpy as np
>>> Sqrt3p = np.nextafter(Sqrt3,999)
>>> Sqrt3p*Sqrt3p
3.0000000000000004

If you could split the difference, you'd have it.

>>> Sqrt3*Sqrt3p
3.0

---

Related Topics