Overflow issues: geomean amean

[Zeda]

The way geometric mean (sqrt(x*y)) and arithmetic mean (a+b)/2 is implemented in this library will erroneously cause overflow with some range of inputs. Neither geometric mean nor arithmetic mean (geomean and amean in this library) should overflow except when one or more arguments is +inf or -inf.

Because these are used to compute the Borchardt-Gauss mean, which is used for natural logarithm and inverse trig and inverse hyperbolics and other routines, all of these will erroneously result in overflow for some range of inputs.

As such, this is an important set of bugs to fix.

---

Ideas:

For amean, if the exponent of either input is close to the maximum, subtract 1 from both exponents, then add the inputs. Otherwise, add the inputs, then subtract 1 from the result exponent. Note: If we go with the first method for all inputs, then we will have erroneous underflow issues.

For geometric mean, it is a bit more complicated. You could compute it as sqrt(x)*sqrt(y), but that erroneously returns NaN if both inputs are negative. It also costs an extra square root.

I think the easiest method is to wrap the existing geomean code in a routine that calculates the output exponent, then sets input exponents to either 0 or 1 as needed. Basically, if e1 and e2 are the respective exponents, then first do (e1+e2)>>1 → etemp. Set e1 to the carry from computing etemp (so 0 or 1) and e2 to 0. Now compute the product, then the square root, then add etemp to the exponent.

Another method that I contemplated a few years ago here might be useful for the lower-precision floats, but that'll require some investigation.