6
$\begingroup$

It is conventional wisdom that programmers should use std::hypot whenever one implements an expression of the form $r = \sqrt{x^2 + y^2}$

In my example, my expression is $r = \sqrt{1 + c^2}$, where $0 \le c \le 1, c \in \Re$. My first instinct is to reach for hypot() to implement - I am interested in both double and float variants. In my problem, for most instances $c \ll 1$.

However, a quick read of Wikipedia suggests I might be wasting my (CPU) time. Since the suggested implementation where $x > y$ is:

$$ r = \sqrt { x^2 + y^2 } \\ = \sqrt { x^2 ( 1 + (y/x)^2) } \\ = |x| \sqrt {1 + (y/x)^2 } $$

However, in my instance $x = 1$ and therefore it is algebraically the same as the naive implementation $r = \sqrt{1 + y^2}$.

What would be the impact of not using hypot() in this problem, and are there any other numerical pitfalls to look for?

$\endgroup$
10
$\begingroup$

Summary: it is cheapest and most accurate to use sqrt(fma(c, c, 1)) if you have FMA, and sqrt(1+c*c) otherwise. In my testing, though, the difference is extremely marginal: of the 1065353216 32-bit floating point numbers $0\leq c\leq 1$, the first formula is better 532509 times (0.05%), and worse 382159 times (0.035%), and neither formula has error worse than 1 ulp, so the obvious sqrt(1+c*c) is good enough.

It is a little bit of a fine point, and Wikipedia can sometimes be a little unreliable on this. One good way to settle these types of questions is to go to a mature library that already implements hypot, such as openlibm (https://github.com/JuliaLang/openlibm).

Quoting the source code comment that explains it (https://github.com/JuliaLang/openlibm/blob/master/src/e_hypot.c):

/* __ieee754_hypot(x,y)
 *
 * Method :                  
 *  If (assume round-to-nearest) z=x*x+y*y 
 *  has error less than sqrt(2)/2 ulp, than 
 *  sqrt(z) has error less than 1 ulp (exercise).
 *
 *  So, compute sqrt(x*x+y*y) with some care as 
 *  follows to get the error below 1 ulp:
 *
 *  Assume x>y>0;
 *  (if possible, set rounding to round-to-nearest)
 *  1. if x > 2y  use
 *      x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
 *  where x1 = x with lower 32 bits cleared, x2 = x-x1; else
 *  2. if x <= 2y use
 *      t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
 *  where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 
 *  y1= y with lower 32 bits chopped, y2 = y-y1.
 *      
 *  NOTE: scaling may be necessary if some argument is too 
 *        large or too tiny
 *
 * Special cases:
 *  hypot(x,y) is INF if x or y is +INF or -INF; else
 *  hypot(x,y) is NAN if x or y is NAN.
 *
 * Accuracy:
 *  hypot(x,y) returns sqrt(x^2+y^2) with error less 
 *  than 1 ulps (units in the last place) 
 */

So if you can compute $1+c^2$ sufficiently accurately (and in your case with $|c|\leq 1$ you will not have over-/underflow), the final result will be accurate.

If you have a modern CPU with a fused multiply-add (FMA) instruction, this is trivial using fma, which most languages' standard math libraries have, so you can use sqrt(fma(c, c, 1)) (cost: just 1 flop plus the cost of a square root), this is even marginally cheaper than what hypot does. The error in fma(c, c, 1) is at most $\frac12$ an ulp with round-to-nearest, so you'll get the same accuracy as hypot, with error $<1$ ulp, so this is the best.

Regarding the formulas that hypot actually uses, I don't really understand what they're doing. It chooses between $1+c^2$ and $2c+(c-1)^2$ depending on whether $c\leq\frac12$. It almost looks like a special case of double-double arithmetic, where you compute a product accurately by writing it as a sum of two numbers, but I'm not sure. I imagine if they could use fma, the formulas would be a lot simpler. In my testing they're at most as accurate as the plain $1+c^2$.

What would happen if you tried $1+c^2$ directly? The relative error of evaluating $\hat w = \mathrm{fl}(1+\mathrm{fl}(c^2))$ is at most $\epsilon_1\leq \frac34\mathrm{ulp}$, which gives the error in $\mathrm{fl}(\sqrt{\hat w})$ as $\sqrt{1+\epsilon_1}-1+\epsilon_2 \leq \tfrac78\mathrm{ulp}$, which is at most a single ulp, so it's accurate too, and as good as hypot.

$\endgroup$
  • $\begingroup$ @Damien You're welcome. Incidentally, it's not obvious from your question if you want this, and it might be overkill, but you could probably increase accuracy all the way to half an ulp by using compensated arithmetic. That's probably close to something like what openlibm's hypot is trying to do. $\endgroup$ – Kirill Aug 31 '17 at 21:48
  • $\begingroup$ @Krill, I have clarified the final question a little - if nothing else but to make it fit your excellent answer. My ultimate motivation is that I need 5 decimal digits of precision when using 32-bit float. However, the $\sqrt{1+c^2}$ is in the inner loop of an $O^3$ operation running on an embedded system, so cycles matter, and there will be other contributors to the error budget. $\endgroup$ – Damien Aug 31 '17 at 22:47
  • $\begingroup$ So sqrt(fma(c,c,1)) or sqrt(1+c*c) (depending on having FMA, and sometimes your compiler can rewrite 1+c*c into fma(c,c,1) on its own, but it can depend on compiler flags) is about as good as it can get. By comparison, hypot does a lot of work to handle the harder cases that you don't have. Also, it usually doesn't really get much better than 1ulp error, most common arithmetic expression you'd normally have are not even as accurate as that. $\endgroup$ – Kirill Aug 31 '17 at 22:55
  • $\begingroup$ Better yet, C++11 now defines std::fma() which should delegate the instruction to hardware if available (which I just checked it is for 32-bit float on my platform) $\endgroup$ – Damien Aug 31 '17 at 22:57
  • $\begingroup$ @Damien Right, but bear in mind that if FMA is not available, then std::fma may default to a slower software implementation. (Although for Float32, it might look like Float32(Float64(a)*Float64(b)+Float64(c)), but it depends.) Here's openlibm's: github.com/JuliaLang/openlibm/blob/master/src/s_fmaf.c (compare with Float64 version: github.com/JuliaLang/openlibm/blob/master/src/s_fma.c), so normally you have to check for CPU support before relying on it too much. $\endgroup$ – Kirill Aug 31 '17 at 22:58

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.