49
$\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.

  • 9
    Voting to reopen.2012-11-17
  • 2
    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](http://www.sciencedirect.com/science/article/pii/S1571066115000456))2016-06-23

1 Answers 1

39

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.

  • 0
    Isn't $\Delta(x) = \mu / M(x)$?2012-11-18
  • 0
    @WimC: Thanks, fixed. Fortunately, it works exactly the same.2012-11-18
  • 0
    "$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 `double`s for equality, the latter is unlikely IMHO.2018-06-05
  • 0
    @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.2018-06-05