square root of square of x equals x

Question

Given a double precision floating-point (non-negative) number x, does square root of its square always equals to itself?

In other words, is there any loss of precision if doing the following:

x = <non-negative double>
y = x^2
z = sqrt(y)

so that:

x == z

I'm not interested in the case when the square becomes infinity or zero, just numbers that fit the double-precision.

score 2 · Answer 1 · edited Jun 20 '20 at 09:12

2

Squaring a number producing a value twice the number of bits in the original values. Hence if x is too large then some bits are lost in x^2 and x cannot be fully recovered from y [Edit: it's still possible to get x from y with proper rounding]. In case of IEEE-754 double precision then if x has more than 26 bits in the significand part then the result of y will be truncated. That's the simplest case.

If x has few significand bits but very large or very small exponent then x^2 might be too large for double precision and will become inf or denormal number, in which case there's no way to recover x.

If x is not too large or too small then sqrt(y) would be equal to x because IEEE-754 standard requires +, -, *, / and sqrt to be properly rounded.

Examples:

#include <iostream>
#include <ios>
#include <iomanip>
#include <cmath>

using std::fixed;
using std::hexfloat;
using std::cout;

int main() {
    double x = 1.25e155;
    double y = x*x;
    cout << hexfloat << "x = " << x << ", y = " << y << ", sqrt(y) = " << sqrt(y) << '\n';
    
    x = 1.25e-155;
    y = x*x;
    cout << hexfloat << "x = " << x << ", y = " << y << ", sqrt(y) = " << sqrt(y) << '\n';
    return 0;
}

edited Jun 20 '20 at 09:12

Community

1
1

answered Mar 13 '17 at 16:23

phuclv

37,963
15
156
475

Could you please provide an example? This is what I thought as well but I was unable to come up with such number. `812387342387487234.98132985798345 == sqrt(812387342387487234.98132985798345^2)` – Ecir Hana Mar 13 '17 at 16:30
@EcirHana First thought: a number so small that its square is not representable and produces zero. – Mar 13 '17 at 16:45
@hvd Yes, infinity and zero would not be recoverable but what about other values? Please, see the edit. – Ecir Hana Mar 13 '17 at 16:48
2

"x cannot be fully recovered from y" <- I don't think this is true: for IEEE 754, round-ties-to-even, etc, one can show that the error in `y` is small enough that taking square root *always* recovers `x`. Roughly: without loss of generality, assume `0.5 < x < 1.0`; show (assuming double-precision format) `|y - x^2| <= 2^-54`, deduce `|sqrt(y) - x| < 2^-54`, hence that under correct rounding `sqrt(y)` rounds to `x`. – Mark Dickinson Mar 13 '17 at 17:28
“if x has more than 21 bits in the significand part then the result of y will be truncated” is not true. Double-precision multiplication, when it does not overflow or underflow, is exact for factors of up to 26 binary digits each. – Pascal Cuoq Mar 13 '17 at 18:56
@PascalCuoq yes that was the result of writing too late at midnight – phuclv Mar 14 '17 at 01:22

score 1 · Accepted Answer · edited Nov 30 '19 at 07:53

1

#include <stdio.h>
#include <math.h>

int main(void) {
  double x = 1.0000000000000001E-160;
  double square = x*x;
  double root = sqrt(square);
  if (root != x) {
    printf("%.20g\n", x);
    printf("%.20g\n", root);
  }
}

outputs

1.0000000000000001466e-160
9.9999443357584897793e-161

What's happening here is that x is large enough that its square is non-zero, but small enough that its square is only representable as a denormalised number, which reduces the available precision.

I get the impression that @MarkDickinson's comment on @LưuVĩnhPhúc's answer is largely correct though. If both x and x*x are positive normalised numbers, then I'm not able to find examples where x != sqrt(x*x) even with a quick brute force (over a few small ranges), though this should not be considered proof.

edited Nov 30 '19 at 07:53

phuclv

37,963
15
156
475

answered Mar 13 '17 at 18:18

3

The property that if `x*x` does not underflow or overflow, then `sqrt(x*x) == x` is proved formally in https://hal.inria.fr/hal-01148409/document – Pascal Cuoq Mar 13 '17 at 18:58
@PascalCuoq That's interesting, looks like I have some reading to do. :) The OP didn't mention the restriction on denormalised numbers that the other question does mention, so it's not an exact duplicate, but that other question plus answer combined do answer this one and your answer there is far better than mine here. – Mar 13 '17 at 19:09

square root of square of x equals x

2 Answers2

Examples:

Linked