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.