ENH: np.hypot(x1: complex, x2: complex) -> complex

[keithbriggs]

hypot currently doesn't work for complex inputs (and this is not documented). I can't see any good reason for this to not be implemented.

>>> import numpy as np
>>> np.hypot(1,1j)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: ufunc 'hypot' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

[mhvk]

What would be wanted? sqrt(|z1|^2+|z2|^2) or without the absolute values? I assume with (i.e., output is real), but np.square does z^2 not |z|^2 (sadly).

[mattip]

I think this not implemented on purpose. Python complex math cmath module does not implement hypot, and as far as I can tell it is not part of the C99 complex spec. Is there an accepted (i.e. standards-based) way to compute this?

[rkern]

And what is the use case for this operation on complexes?

[keithbriggs]

We would simply want sqrt(x1**2+x2**2), but computed in a numerically stable way without risk of overflow. These are the same reasons the ordinary hypot is included in most math libraries. The expression sqrt(x1**2+x2**2) comes up frequently in numerical linear algebra, where the inputs can be complex. My use case is less common, in some power series calculations involving Bessel functions.

[mhvk]

So, you actually want x1**2, not |x1|^2? Could you give an actual example for the complex case?

[keithbriggs]

If you want - but it's very specialized and won't mean much. If you want to integrate J0(b*t)*cos(c*x)/t, where t=sqrt(a**2-x**2), over x=0...a, and where b and c are complex, then there is a closed formula in terms of hypot(b,c).

[rkern]

I'm curious about the usage in linear algebra. I'd expect to see a complex version in LAPACK, but I only see one for reals.

I guess I wouldn't reject a contribution of one that adapts the simple real-valued scheme that pulls the larger argument out of the square root, but I'd want to see some tests that confirm that it does help, which I think might be original research that I'm not sure anyone's going to volunteer to do. There will also need to be tests to ensure that the rearrangement is doing the right thing with respect to branch cuts and signed zeros and all of that annoying complex stuff.

[rkern]

Julia and MATLAB implement hypot(z1, z2) as sqrt(abs(z1)**2 + abs(z2)**2), FWIW.

[seberg]

Given the prior art, it seems to me that the only definition that we could consider adding is the one returning a real value, and that is not the one requested. (Whether it was thought through or not, but the Julia discussion seems to have a tendency towards that definition at least in hindsight.)

The way to go for you is probably a custom ufunc outside of NumPy. If you use such tools, you could and don't want to dive into how that works, you could also write it using a tool like @numba.vectorize.

[keithbriggs]

Ok, I can manage with my own implementation. But could at least the docs be updated to mention that hypot only works for floats?

[seberg]

But could at least the docs be updated to mention that hypot only works for floats?

Happy if someone has a look if it makes sense, we have probably quite a few other similar cases...

TBH, I think it is more useful to update the error message to include the complex information, chances are you see the error before you read the docs anyway. (But I am pretty sure we have an issue elsewhere for that.)