Is there any point to using hypot() for $\sqrt{1+c^2}$, $0 \le c \le 1$ for real numbers

Question

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?

Answer 1

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.