Add gamma function, error functions and other C99 math.h functions to math module

[terryjreedy]

BPO	3366
Nosy	@tim-one, @rhettinger, @terryjreedy, @mdickinson, @ned-deily, @stevendaprano
Files	math_private.h: header for erf, erfc, lgamma and tgamma pymath.c.diff mathmodule.c.diff pymath.h.diff test_math.py.diff math.rst.diff mathmodule.diff: patch for erf, erfc, lgamma and gamma pymath.c.diff: erf and gamma patches without math_private.h gamma.patch gamma3.patch: Implementation of math.gamma gamma4.patch gamma5.patch lgamma1.patch gamma-failures.txt pymath.h.diff: Euler constant mathmodule.c.diff: expm1 erf.py: error function mathmodule.c.diff erf.patch erf_c.patch: Add erf and erfc functions to math module. expm1.patch

^{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:

assignee = 'https://github.com/mdickinson'
closed_at = <Date 2010-01-09.18:57:05.131>
created_at = <Date 2008-07-15.17:37:53.559>
labels = ['type-feature', 'library']
title = 'Add gamma function, error functions and other C99 math.h functions to math module'
updated_at = <Date 2010-01-09.18:57:05.131>
user = 'https://github.com/terryjreedy'

bugs.python.org fields:

activity = <Date 2010-01-09.18:57:05.131>
actor = 'mark.dickinson'
assignee = 'mark.dickinson'
closed = True
closed_date = <Date 2010-01-09.18:57:05.131>
closer = 'mark.dickinson'
components = ['Library (Lib)']
creation = <Date 2008-07-15.17:37:53.559>
creator = 'terry.reedy'
dependencies = []
files = ['10955', '10973', '10974', '10975', '10976', '10977', '10982', '11005', '14844', '14927', '14940', '14943', '15008', '15152', '15488', '15489', '15490', '15491', '15539', '15540', '15574']
hgrepos = []
issue_num = 3366
keywords = ['patch', 'needs review']
message_count = 54.0
messages = ['69698', '70106', '70108', '70120', '70125', '70128', '70131', '70232', '70233', '70241', '70262', '70266', '70267', '70270', '70361', '70412', '77075', '92294', '92844', '92928', '92935', '92936', '92939', '92944', '92948', '93228', '93370', '93373', '94165', '94166', '94168', '94169', '94170', '94171', '94172', '94175', '94224', '96108', '96109', '96110', '96111', '96266', '96267', '96268', '96332', '96333', '96485', '96497', '96599', '96600', '97127', '97131', '97133', '97160']
nosy_count = 8.0
nosy_names = ['tim.peters', 'rhettinger', 'terry.reedy', 'mark.dickinson', 'ned.deily', 'stutzbach', 'nirinA', 'steven.daprano']
pr_nums = []
priority = 'normal'
resolution = None
stage = 'commit review'
status = 'closed'
superseder = None
type = 'enhancement'
url = 'https://bugs.python.org/issue3366'
versions = ['Python 2.7', 'Python 3.2']

[terryjreedy]

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=src

[nirinA]

here 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]

and this is the header.
it is stolen from glibc.

[mdickinson]

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.

[rhettinger]

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]

On 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.

[rhettinger]

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]

Mark 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.

[rhettinger]

Why isn't tgamma() simply named gamma()? The t prefix does nothing for
me except raise questions.

[stutzbach]

I think he just carried the names over from C, where:

tgamma() is the "true gamma function"
lgamma() is the log of the gamma function

and 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.

[mdickinson]

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]

ouch!
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 header

the 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.