On Sep 7, 2019, at 18:07, firstname.lastname@example.org wrote:
Since the gamma function (see https://en.wikipedia.org/wiki/Gamma_function) is defined for complex numbers too, wouldn't it be normal to have it implemented on the `cmath` module? The `math` module has this function, but it doesn't support complex numbers (as expected).
I believe this is just because math and cmath started off mostly thin wrappers around the C stdlib functions, and C had a gamma for double but not for complex.
(Actually, IIRC, it’s a bit more complicated, more like “what POSIX guaranteed that Windows also provided as of the mid 90s”, which goes beyond the C89 standard but doesn’t match the C99 standard.)
But that doesn’t mean it couldn’t be added. The `math` and `cmath` modules have added lots of things over the decades that are useful and aren’t thin wrappers, like `fsum` and `isclose`, and hasn’t included everything that was added in later versions of C or POSIX.
In fact, math.gamma itself is effectively one of those things; it’s no longer a wrapper around a libm function, but a custom Lancsoz implementation that gives the same values on all platforms.
But on the third hand, you’d presumably want cmath.gamma to give the same values as math.gamma for compatible inputs. And the algorithm in mathmodule.c doesn’t translate smoothly to complex. I believe libraries like GSL deal with this by using a different algorithm that does translate smoothly, for both real and complex values, but that means giving up either speed or accuracy or both.
I did some quick tests on scipy.special.gamma. It takes roughly 12x as long as math.gamma, and gives larger errors (within 3 ulp instead of 1 even near 1, and as much as 117ulp instead of no more than 10 for extreme values), and doesn’t give the same values for float and complex (e.g., for 1/3, they differ by 6ulp). Is that good enough for the stdlib?