Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
5  
Voting to reopen. –  user17762 Nov 17 '12 at 18:32
add comment

1 Answer

up vote 13 down vote accepted

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.

share|improve this answer
    
Isn't $\Delta(x) = \mu / M(x)$? –  WimC Nov 18 '12 at 7:58
    
@WimC: Thanks, fixed. Fortunately, it works exactly the same. –  maaartinus Nov 18 '12 at 12:03
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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