50
$\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$
  • 9
    $\begingroup$ Voting to reopen. $\endgroup$ – user17762 Nov 17 '12 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 '16 at 0:19
40
$\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$
  • $\begingroup$ Isn't $\Delta(x) = \mu / M(x)$? $\endgroup$ – WimC Nov 18 '12 at 7:58
  • $\begingroup$ @WimC: Thanks, fixed. Fortunately, it works exactly the same. $\endgroup$ – maaartinus Nov 18 '12 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 '18 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 '18 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.