Repository navigation
Add gamma function, error functions and other C99 math.h functions to math module #47616
Description
Activity
From Py3000 list: in C99 standard math library, but need code when not
yet in particular implementation.
Dickinson: open and assign to me
Stutzbach: I suggest using the versions from newlib's libm. They contain
extensive comments explaining the math and have a generous license,
e.g.,:
http://sourceware.org/cgi-bin/cvsweb.cgi/~checkout~/src/newlib/libm/math/s_erf.c?rev=1.1.1.1&cvsroot=srcnirinA commented
on Jul 21, 2008 nirinAmannequinMannequinhere is an attempt to make erf, erfc, lgamma and tgamma
accessible from math module.
erf and erfc are crude translation of the code pointed
out by Daniel Stutbach.
lgamma is also taken from sourceware.org, but doesn't
use the reentrant call for the sign.i tested only on gcc-4.3.1/linux-2.6.26.
i'll write a testsuite soon.some results:
Python 3.0b2 (r30b2:65080, Jul 21 2008, 13:13:13)
[GCC 4.3.1] on linux2>>> print(math.tgamma(0.5)) 1.77245385091 >>> print(math.sqrt(math.pi)) 1.77245385091 >>> print(math.tgamma(1.5)) 0.886226925453 >>> print(math.sqrt(math.pi)/2) 0.886226925453 >>> print(math.sqrt(math.pi)*3/4) 1.32934038818 >>> print(math.tgamma(2.5)) 1.32934038818 >>> for i in range(14): print(i, math.lgamma(i), math.tgamma(i)) 0 inf inf 1 0.0 1.0 2 0.0 1.0 3 0.69314718056 2.0 4 1.79175946923 6.0 5 3.17805383035 24.0 6 4.78749174278 120.0 7 6.57925121201 720.0 8 8.52516136107 5040.0 9 10.6046029027 40320.0 10 12.8018274801 362880.0 11 15.1044125731 3628800.0 12 17.5023078459 39916800.0 13 19.9872144957 479001600.0 >>> for i in range(-14,0): print(i/5, math.lgamma(i/5), math.tgamma(i/5)) -2.8 0.129801291662 -1.13860211111 -2.6 -0.118011632805 -0.888685714647 -2.4 0.102583615968 -1.10802994703 -2.2 0.790718673676 -2.20498051842 -2.0 inf inf -1.8 1.15942070884 3.18808591111 -1.6 0.837499812222 2.31058285808 -1.4 0.978052353322 2.65927187288 -1.2 1.57917603404 4.85095714052 -1.0 inf -inf -0.8 1.74720737374 -5.73855464 -0.6 1.30750344147 -3.69693257293 -0.4 1.31452458994 -3.72298062203 -0.2 1.76149759083 -5.82114856863 >>> math.erf(INF) 1.0 >>> math.erf(-INF) -1.0 >>> math.erfc(INF) 0.0 >>> math.erfc(-INF) 2.0 >>> math.erf(0.6174) 0.61741074841139409
nirinA commented
on Jul 21, 2008 nirinAmannequinMannequinand this is the header.
it is stolen from glibc.Thanks for the patch! I probably won't get to this properly until after
2.6 final, but it won't get lost. It seems like there's pretty good
support for adding these functions.By the way, there's already a setup in Python 2.6 (ad 3.0) for substitutes
for missing math functions: take a look at the file Python/pymath.c,
which provides inverse hyperbolic functions for those platforms that need
them, as well as the bits in configure.in that check for these functions.
This is the really the right place for the tgamma/lgamma/erf/erfc code.So a full patch for this should touch at least Python/pymath.c,
Modules/mathmodule.c, configure.in, Lib/test/test_math.py, and
Doc/Library/math.rst. If you have time to fill in any of these pieces, it
would be greatly appreciated.Since we're not in a hurry for Py2.7 and 3.1, I would like to this
kicked around a bit on the newsgroup and in numpy forums (personally, I
would also post a pure python equivalent to the ASPN cookbook for
further commentary).There are several different approximations to choose from. Each of
them has their own implications for speed and accuracy. IIRC, the one
I used in test.test_random.gamma() was accompanied by a proof that its
maximum error never exceeded a certain amount. I think there were some
formulas that made guarantees only over a certain interval and others
that had nice properties in the first and second derivatives (one that
don't have those properties can throw newtonian solvers wildly off the
mark).Let's let the community use its collective wisdom to select the best
approach and not immediately commit ourselves to the one in this patch.At one point, Tim was reluctant to add any of these functions because
it is non-trivial to do well and because it would open a can of worms
about why python gets a different answer (possibly not as close to
true) and some other favorite tool (matlab, excel, numpy, and
engineering calculator, etc).FWIW, I changed the 2.6 version of test_random.gamma() to take
advantage of msum() to increase its accuracy when some of the summands
have nearly cancelling values.stutzbach commented
on Jul 21, 2008 stutzbachmannequinMannequinOn Mon, Jul 21, 2008 at 5:34 PM, Raymond Hettinger
<report@bugs.python.org> wrote:There are several different approximations to choose from. Each of
them has their own implications for speed and accuracy.FWIW, the current patch dynamically selects one of a few different
approximation methods based on the input value.+1 on kicking it around and comparing the output with that from other
mathematics tools.It would be nice if we knew the error bounds for each of the
approximation methods. Do we know how the coefficients were generated?When switching from one method to another, it might be nice to have a
range where the results slowly blend from one method to another:
if x < a: return f1(x)
if x > b: return f2(x)
t = (x - a) / (b - a)
return (1-t)*f1(x) + t*f2(x)
This minimizes discontinuities in the first and second derivatives.The lowword/highword macros look to be tightly tied to the internal
processor representation for floats. It would be more portable and
maintainable to replace that bit twiddling with something based on frexp
().nirinA commented
on Jul 25, 2008 nirinAmannequinMannequinMark Dickinson :
So a full patch for this should touch at least Python/pymath.c,
Modules/mathmodule.c, configure.in, Lib/test/test_math.py, and
Doc/Library/math.rst.
here are patches for Python/pymath.c, Modules/mathmodule.c,
Lib/test/test_math.py, Doc/Library/math.rst and also Include/pymath.h.
i haven't touched configure.in as i've erf, erfc, lgamma and tgamma
on my platform, so the patches can be tested.
the doc may need better wording and test_math.py some additionnal tests.Raymond Hettinger:
It would be nice if we knew the error bounds for each of the
approximation methods. Do we know how the coefficients were generated?
it will be really great if we're able to regenerate all these
coefficients. this may be done by following the comments in
the original codes, but not that easy and not sure one will obtain
the same things.The lowword/highword macros look to be tightly tied to the internal
processor representation for floats. It would be more portable and
maintainable to replace that bit twiddling with something based on frexp
it seems possible to avoid the use of these macros.
i'll try to find time to dig up these possibilities.Why isn't tgamma() simply named gamma()? The t prefix does nothing for
me except raise questions.stutzbach commented
on Jul 25, 2008 stutzbachmannequinMannequinI think he just carried the names over from C, where:
tgamma() is the "true gamma function"
lgamma() is the log of the gamma functionand gamma() might be tgamma() or lgamma() depending on which C library
you use (sigh).I'm +1 on making gamma() be the true gamma function and not carrying
over this brain-damage to Python.I'm +1 on making gamma() be the true gamma function and not carrying
over this brain-damage to Python.+1 from me, too.
nirinA commented
on Jul 25, 2008 nirinAmannequinMannequinouch!
i was asleep, at 3AM, when i edited the issue
and didn't really know what i was doing.
i just see that i removed, instead of edited the
mathmodule.diff and didn't check after.
so here it is again, with the url for original code
added.if you want to test the erf, erfc, lgamma and gamma
please use this mathmodule.diff and the math_private.h headerthe other diffs are incomplete and don't work, unless you have
erfs and gammas functions implemented in your platform.
i'll try to work on these files this weekend and
will rename tgamma to gamma.30 remaining items
nirinA: thanks for prodding this issue. Yes, it is still alive (just).
:)About adding the Euler constant: I agree this should be added. Do you
have time to expand your patch to include docs and a test or two?For expm1, I was planning to use the libm function when it exists (which
I expect it will on most platforms), and just use a quick hack based on
the identityexpm1(x) = sinh(x/2)*exp(x/2)
for platforms where it's missing.
I'll look at the erf and erfc implementations.
Daniel: re incomplete beta and gamma functions, I'd prefer to limit
this particular issue to the C99 math functions. It's difficult to know
where to draw the line, but it seems to me that these are really the
domain of numpy/scipy. And they involve multiple input arguments, which
makes them *darned hard* to write and test! Anyway, if you'd like to
see those added to Python's math module, maybe this could be brought up
on python-dev to see how much support there is.expm1(x) = sinh(x/2)*exp(x/2)
Should be 2*sinh(x/2)*exp(x/2).
lgamma patch committed to trunk (with looser tests) in r76755.
nirinA: thanks for the erf patch. It's fine as far as it goes; the main
technical issue is that the series development only converges at a
sensible rate for smallish x; for larger x something more is needed.Here's a prototype patch for the erf and erfc functions that uses a power
series (not quite the same as the series nirinA used; by pulling out a
factor of exp(-x*2) you can make all coefficients in the remaining series
positive, which is good for the numerics) for small x, and a continued
fraction expansion for larger x. The erf and erfc functions are
currently written in Python. I'll do some more testing (and possibly
tweaking) of this patch and then translate the Python to C.Here's the C version.
Here's a patch to add expm1. Rather than putting the code in pymath.c,
which is the wrong place (see issue bpo-7518), I've added a new file
Modules/_math.c for this; log1p, atanh, etc. should also eventually be
moved from pymath.c into this file.I'm fairly confident that the maths and numerics are okay, but less
confident about my changes to the build structure; if anyone has a chance
to test this patch, especially on Windows, I'd welcome feedback. I've
tested it on OS X, and will test on Linux.I committed the patch for expm1 in r76861 (trunk) and r76863 (py3k), with
one small change: instead of using 2 * exp(x/2) * sinh(x/2), expm1 is now
computed using a method due to Kahan that involves only exp and log. This
seems a little safer, and guards against problems on platforms with a
naive implementation of sinh. (I *hope* no such platforms exist, but
suspect that's wishful thinking. :-)Thanks to Eric Smith for testing the build system changes on Windows.
Added erf and erfc in r76879 (trunk), r76880 (py3k).
Oops. That's r76878 (trunk) and r76879 (py3k).
The last two functions to consider adding are exp2 and log2. Does anyone
care about these?stutzbach commented
on Jan 2, 2010 stutzbachmannequinMannequinAny time I've ever needed log2(x), log(x)/log(2) was sufficient.
In Python, exp2(x) can be spelled 2.0**x. What would exp2(x) gain us?
In Python, exp2(x) can be spelled 2.0**x. What would exp2(x) gain us?
Not much, I suspect. :)
I'd expect (but am mostly guessing) exp2(x) to have better accuracy than
pow(2.0, x) for some math libraries; I'd further guess that it's somewhat
more likely to give exact results for (small) integral x.Similarly for log2: log2(n) should be a touch more accurate than
log(n)/log(2), and the time you're most likely to notice the difference is
when n is an exact power of 2.But we've already got the 'bit_length' method for integers, which fills
some of the potential uses for log2. So unless there's a feeling that
these functions are needed, I'd rather leave them out.Will close this unless there's an outcry of support for exp2 and log2.
nirinA, if you're still interested in adding the euler constant, please open another issue.
Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.
Show more details
GitHub fields:
bugs.python.org fields: