56
$\begingroup$

I verified experimentally that in Java the equality

Math.sqrt(x*x) = x

holds for all long x such that x*x doesn't overflow. Here, Java long is a $64$ bit signed type and double is a IEEE binary floating point type with at least $53$ bits mantissa and sufficiently long exponent.

Mathematically, there are two imprecise functions involved:

  • Conversion from long to double which loses precision due to the mantissa being only $53$ bits where $63$ bits would be needed. This operation is guaranteed to return the closest representable result.
  • Computing square root, which is also guaranteed to return the closest representable result.

Mathematically, this can be expressed like

$$ \mathop\forall_{x \in {\mathbb N} \atop x \le 3037000499} \mathrm{round}\left(\sqrt{\mathrm{round}(x^2)}\right) = x $$

where $\mathrm{round}$ is the rounding function from $\mathbb R$ into the set of all numbers representable as double.

I'm looking for a proof since no experiment can assure it works across all machines.

$\endgroup$
2
  • 9
    $\begingroup$ Voting to reopen. $\endgroup$
    – user17762
    Nov 17, 2012 at 18:32
  • 2
    $\begingroup$ See also: Sylvie Boldo, "Stupid is as Stupid Does: Taking the Square Root of the Square of a Floating-Point Number". Electronic Notes in Theoretical Computer Science Volume 317, 18 November 2015, Pages 27-32 (online) $\endgroup$
    – njuffa
    Jun 23, 2016 at 0:19

1 Answer 1

45
$\begingroup$

The idea is simple: Find upper and lower bounds for

$$X := \sqrt{\mathrm{round}(x^2)}$$

and show that $\mathrm{round}(X) = x$.


Let $\mathrm{ulp}(x)$ denote the unit of least precision at $x$ and let $E(x)$ and $M(x)$ denote the exponent and mantissa of $x$, i.e.,

$$x = M(x) \cdot 2^{E(x)}$$

with $1 \le M(x) < 2$ and $E(x) \in \mathbb Z$. Define

$$\Delta(x) = \frac{\mathrm{ulp}(x)}x = \frac{\mu \cdot 2^{E(x)}}x = \frac\mu{M(x)}$$

where $\mu=2^{-52}$ is the machine epsilon.

Expressing the rounding function by its relative error leads to

$$X = \sqrt{(1+\epsilon) \cdot x^2} = \sqrt{(1+\epsilon)} \cdot x < \big( 1+\frac\epsilon2 \big) \cdot x$$

We know that $|\epsilon| \le \frac12\Delta(x^2)$ and get (ignoring the trivial case $x=0$)

$$\frac Xx < 1 + \frac{\Delta(x^2)}4 = 1 + \frac\mu{4 M(x^2)}$$


By observing $M(x)$ and $M(x^2)$ e.g. over the interval $[1, 4]$, it can be easily be shown that $\frac{M(x)}{M(x^2)} \le \sqrt2$ which gives us

$$\frac Xx < 1 + \frac{\mu\sqrt2}{4 M(x)}$$

and therefore

$$X < x + \frac{\sqrt2}4 \frac{\mu}{M(x)} \cdot x < x + \frac12 \mathrm{ulp}(x)$$


Analogously we get the corresponding lower bound. Just instead of

$$\sqrt{(1+\epsilon)} < \big( 1+\frac\epsilon2 \big)$$

we use something like

$$\sqrt{(1-\epsilon)} > \big( 1 - (1+\epsilon) \cdot \frac\epsilon2 \big)$$

which suffices, since we used a very generous estimate ($\sqrt2/4<\frac12$) in the last step.


Because of $|X-x|$ being smaller than $\frac12 \mathrm{ulp}(x)$, $x$ is the double closest to $X$, therefore $\mathrm{round}(X)$ must equal to $x$, q.e.d.

$\endgroup$
4
  • $\begingroup$ Isn't $\Delta(x) = \mu / M(x)$? $\endgroup$
    – WimC
    Nov 18, 2012 at 7:58
  • $\begingroup$ @WimC: Thanks, fixed. Fortunately, it works exactly the same. $\endgroup$
    – maaartinus
    Nov 18, 2012 at 12:03
  • $\begingroup$ "$x$ is the double closest to $X$" - $x$ is a long, not a double, so I'm having difficulty following this line of the argument. Did you mean "$X$ is the double closest to $x$"? If so, that does not imply that $X = x$ unless you either establish that $\mathrm{ulp}(x) < 2$ or assume that we're converting $x$ into double (rather than converting $X$ into long). Since most programmers have a ritualistic fear of comparing doubles for equality, the latter is unlikely IMHO. $\endgroup$
    – Kevin
    Jun 5, 2018 at 6:49
  • $\begingroup$ @Kevin I guess, I indeed meant that "$X$ is the double closest to $x$". But yes, $\mathrm{ulp}(x)<2$ as $x$ is small. I wrote above: long x such that x*x doesn't overflow. $\endgroup$
    – maaartinus
    Jun 5, 2018 at 11:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.