From Py3000 list: in C99 standard math library, but need code when notyet in particular implementation.  Dickinson: open and assign to meStutzbach: I suggest using the versions from newlib's libm.  They containextensive 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

here is an attempt to make erf, erfc, lgamma and tgammaaccessible from math module.erf and erfc are crude translation of the code pointedout by Daniel Stutbach.lgamma is also taken from sourceware.org, but doesn'tuse 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 inf1 0.0 1.02 0.0 1.03 0.69314718056 2.04 1.79175946923 6.05 3.17805383035 24.06 4.78749174278 120.07 6.57925121201 720.08 8.52516136107 5040.09 10.6046029027 40320.010 12.8018274801 362880.011 15.1044125731 3628800.012 17.5023078459 39916800.013 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

and 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 filePython/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 leastPython/pymath.c,Modules/mathmodule.c, configure.in,Lib/test/test_math.py, andDoc/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.

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 differentapproximation methods based on the input value.+1 on kicking it around and comparing the output with that from othermathematics 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().

Mark Dickinson :> So a full patch for this should touch at leastPython/pymath.c, >Modules/mathmodule.c, configure.in,Lib/test/test_math.py, and >Doc/Library/math.rst.here are patches forPython/pymath.c,Modules/mathmodule.c,Lib/test/test_math.py,Doc/Library/math.rst and alsoInclude/pymath.h.i haven't touched configure.in as i've erf, erfc, lgamma and tgammaon 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 thesecoefficients. this may be done by following the comments inthe original codes, but not that easy and not sure one will obtainthe 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 frexpit 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.

I 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 libraryyou use (sigh).I'm +1 on making gamma() be the true gamma function and not carryingover 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.

ouch!i was asleep, at 3AM, when i edited the issueand didn't really know what i was doing.i just see that i removed, instead of edited themathmodule.diff and didn't check after.so here it is again, with the url for original codeadded.if you want to test the erf, erfc, lgamma and gammaplease use this mathmodule.diff and the math_private.h headerthe other diffs are incomplete and don't work, unless you haveerfs and gammas functions implemented in your platform.i'll try to work on these files this weekend andwill rename tgamma to gamma.

Can you also implement blending of approximations:  (1-t)*f1(x) + t*f2(x)

On Fri, Jul 25, 2008 at 1:37 PM, Raymond Hettinger<report@bugs.python.org>wrote:> Raymond Hettinger <rhettinger@users.sourceforge.net> added the comment:>> Can you also implement blending of approximations:  (1-t)*f1(x) + t*f2> (x)>Is this necessary?  Are the approximations significantly different near thetransition points?(Haven't had a chance to download the patch and try it myself as a I'mgetting on a plane in the morning--sorry)

i wrote pure Python equivalent of the patch and willpost it to the ASPN cookbook as suggested. i'm wonderingif i have to upload the code here too, or just the linkwill suffice?Raymond Hettinger:> Can you also implement blending of approximations:> (1-t)*f1(x) + t*f2(x)i do this with the Python code for now, and my question is:how to choose the values for the borders?with erf, for example, one has typically 5 domainsof approximation and one can write something like:    def erf(x):        f = float(x)        if (f == inf):            return 1.0        elif (f == -inf):            return -1.0        elif (f == nan):            return nan        else:            if (abs(x)<0.84375):                return erf1(x)            elif (0.84375<=abs(x)<1.25):                return erf2(x)            elif (1.25<=abs(x)<2.857142):                return erf3(x)            elif (2.857142<=abs(x)<6):                return erf4(x)            elif (abs(x)>=6):                return erf5(x)implementing the blending of approximations,one may write:    def blend(x, a, b, f1, f2):        if x < a:            return f1(x)        if x > b:            return f2(x)        return (((b-x)*f1(x))-((a-x)*f2(x)))/(b-a)and the definition becomes:    def erf(x):f=float(x)if (abs(x)<0.83):return erf1(x)elif (0.83<=abs(x)<0.85):return blend(x,0.83,0.85,erf1,erf2)elif (0.85<=abs(x)<1.24):return erf2(x)elif (1.24<=abs(x)<1.26):return blend(x,1.24,1.26,erf2,erf3)elif (1.26<=abs(x)<2.84):return erf3(x)elif (2.84<=abs(x)<2.86):return blend(x,2.84,2.86,erf3,erf4)elif (2.86<=abs(x)<5.99):return erf4(x)elif (5.99<=abs(x)<6.01):return blend(x,5.99,6.01,erf4,erf5)elif (abs(x)>=6.01):return erf5(x)but the choice of these values is somewhat arbitrary.or i completely misunderstood what you mean by "blending of approximations"!for erf functions, blending of approximations, as i understand it,may be implemented easily as the functions are monotonics.for gammas, this will be a bit complicated, especially fornegative values, where there are extremums for half integers anddiscontinuities for integers.in the other hand, i'm wondering if these coefficients had been chosenwith optimization of discontinuities already taken into account.

pure Python codes are posted to ASPN:http://code.activestate.com/recipes/576391/http://code.activestate.com/recipes/576393/i also rewrote the patch so that they don'tuse the high_word/low_word macros anymore.

There have been requests to add other C99 libm functions to Python:expm1 is requested inissue 3501, and log2 came up in theissue 3724 discussion.  I'm absorbing those issues into this one.With these two, and the gamma and error functions, we're pretty close to implementing all the functions in C99s math.h.  Seems to me that we might as well aim for completeness and implement the whole lot.

Finally got around to looking at this properly.Here's a first attempt at a patch to add gamma, lgamma, erf and erfc to the math module.  It uses the pymath.c.diff code from nirinA, and adds docs, tests, and the extra infrastructure needed for the build system.The code has gone in a new file calledModules/math_cextra.c.  The old approach of sticking things in pymath.c really isn't right, since that ends up putting the extra code into the core Python executable instead of the math module.There are currently many test failures;  some of these are undoubtedly due to the tests being too strict, but some of them show problems with the code for these functions.  It should be possible to have erf and erfc produce results accurate to within a few ulps (<50, say) over the entire extended real line;  getting gamma and lgamma accurate for positive inputs should also be perfectly possible.  Negative arguments seem to be more difficult to get right.To test the code on non-Windows, build with something like:CC='gcc -DTEST_GAMMA' ./configure && makethe '-DTEST_GAMMA' just forces the build to use nirinA's code;  without this, the math module will use the libm functions if they're available.This is just work-in-progress;  apart from the test failures, there's still a lot of polishing to do before this can go in.The patch is against the trunk.

Here's a more careful patch for just the gamma function.  It's a fairly direct implementation of Lanczos' formula, with the choice of parameters (N=13, g=6.024680040776729583740234375) taken from the Boost library.  In testing of random values and known difficult values, it's consistently giving errors of < 5ulps across its domain.  This is considerably better than the error from the exp(lgamma(x)) method, since the exp call alone can easily introduce errors of hundreds of ulps for larger values of x.I'd appreciate any feedback, especially if anyone has the chance to test on Windows:  I'd like to know whether I got the build modifications right.This patch *doesn't* use the system gamma function, even when available.  At least on OS X 10.5.8, the supplied gamma function has pretty terrible accuracy, and some severe numerical problems near negative integers.  I don't imagine Linux is much better.Once this is working, I'll concentrate on lgamma.  Then erf and erfc.

New patch for gamma , with some tweaks: - return exact values for integral arguments: gamma(1) through gamma(23) - apply a cheap correction to improve accuracy of exp and pow   computations - use a different form of the reflection formula:     gamma(x) = -pi/sinpi(x)/x/gamma(x)   (the usual reflection formula has accuracy problems for    x close to a power of 2;  e.g., x in (-64,-63) or x    in (-128, -127)) - avoid duplication formula for large negative arguments - add a few extra testsOn my machine, testing with approx. 10**7 random samples, this version achieves an accuracy of <= 10 ulps across the domain (comparing with correctly-rounded results generated by MPFR).  Limiting the test to arguments in the range (-256.0, 1/256.0] + [1/256.0, 256.0) (with each float in that range equally likely), the error in ulps from 10**6 samples has mean -0.104 and standard deviation 1.230.I plan to check this in in a week or two.  Feedback welcome!  It would be especially useful to know whether this patch compiles correctly on Windows.

I'll test your patch on Windows.  Are you working against the trunk orthe py3k branch?

Thanks!  The patch is against the trunk.  (It doesn't quite apply cleanly to py3k, but the changes needed to make it do so should be minimal.)Hmm. Rereading my previous comment, I seem to have a blindness for negative signs:    gamma(x) = -pi/sinpi(x)/x/gamma(x)should have been    gamma(-x) = -pi/sinpi(x)/x/gamma(x)and    (-256.0, 1/256.0] + [1/256.0, 256.0)should have been    (-256.0, -1/256.0] + [1/256.0, 256.0)

I'm only setup to test the official Windows build setup under PCbuild. I'm not testing the legacy build setups underPC/VC6,PC/VS7.1, orPC/VS8.0.The patch against the trunk failed forPC/VC6/pythoncore.dsp.  I don'tneed that to build, though.Your patch built fine, with one warning about a loss of precision on thefollowing line in sinpi():    n = round(2.0*y);I recommend explicitly casting the result of round() to int, to get ridof the warning and because explicit is better than implicit. :)I ran the following tests, all of which passed:PCbuild/python.exeLib/test/regrtest.py -v test_mathI appreciate your work on this patch.  Let me know if there's anythingelse I can do.

Many thanks, Daniel.> The patch against the trunk failed forPC/VC6/pythoncore.dsp.  I don't> need that to build, though.I've no idea why that would happen.  A line-ending issue, perhaps?  If it doesn't stop me committing the change, then perhaps it's not a problem.> Your patch built fine, with one warning about a loss of precision on> the following line in sinpi():>>    n = round(2.0*y);>>I recommend explicitly casting the result of round() to int, to get>rid of the warning and because explicit is better than implicit. :)Done.  Thanks!

gamma5.patch.  very minor changes from the last patch: - add (int) cast suggested by Daniel Stutzbach -Misc/NEWS entry - ..versionadded in docs - minor corrections to some comments

Gamma function added inr75117 (trunk),r75119 (py3k).I'm working on lgamma.

A first attempt for lgamma, using Lanczos' formula.In general, this gives very good accuracy.  As the tests show, there's one big problem area, namely when gamma(x) is close to +-1, so that lgamma(x) is close to 0.  This happens for x ~ 1.0,  x ~ 2.0, and various negative values of x (mathematically, lgamma(x) hits zero twice in each interval (n-1, n) for n <= -2).  In that case the accuracy is still good in absolute terms (around 1e-15 absolute error), but terrible in ulp terms.I think it's reasonable to allow an absolute error of a few times the machine epsilon in lgamma:  for many uses, the lgamma result will beexponentiated at some point (often after adding or subtracting other lgamma results), and at that point this absolute error just becomes a small relative error.

I agree with your reasoning.Any chance of also getting the incomplete beta and gamma functions,needed to compute the CDF of many commonly used probabilitydistributions? (Beta, Logarithmic, Poisson, Chi-Square, Gamma, etc.)

FYI, approximately 20 of the gamma test cases fail on PPC Macs.  Attached are snippets from regrtest runs with trunk and with py3k, both on a G4 ppc (32-bit) running OS X 10.5.  Identical failures for trunk (did not try py3k) were observed on a G3 ppc with OS X 10.4.

FYI, mysterious numeric differences on PPC are often due to the Ccompiler generated code to use the "fused multiply-add" HW instruction. In which case, find a way to turn that off :-)

For the record, these OS X installer build scripts (running on 10.5) results in these gcc options:gcc-4.0 -arch ppc -arch i386 -isysroot /Developer/SDKs/MacOSX10.4u.sdk -fno-strict-aliasing -fno-common -dynamic -DNDEBUG -g -O3 -I/tmp/_py/libraries/usr/local/include -I. -IInclude -I/usr/local/include -I/private/tmp/_t/Include -I/private/tmp/_py/_bld/python -c /private/tmp/_t/Modules/mathmodule.c -o build/temp.macosx-10.3-fat-3.2/private/tmp/_t/Modules/mathmodule.oand:$ gcc-4.0 --versionpowerpc-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5493)

Hmm.  It looks as though the gamma computation itself is fine;  it's the accuracy check that's going wrong.I'm suspecting an endianness issue and/or a struct module bug.Ned, do you still get these failures with a direct single-architecture build on PPC?  Or just for the universal build?

Adding-mno-fused-maddwould be worth trying.  It usually fixes PPC bugs ;-)

It was a stupid Mark bug.  I'd missed out a '<' from the struct.unpack call that was translating floats into integers, for the purposes of computing ulps difference.  Fixed inr75454 (trunk),r75455 (py3k).Thanks for reporting this, Ned.

BTW, the gamma computation shouldn't be at all fragile, unless I've messed up somewhere:  the Lanczos sum is evaluated as a rational function in which both numerator and denominator have only positive coefficients, so there's little danger of nasty cancellations boosting the relative error.I'd be surprised if use of fma instructions made more than 1 or 2 ulps difference to the result, but I'll take a look.

Mark, you needn't bother:  you found the smoking gun already!  From yourdescription, I agree it would be very surprising if FMA made asignificant difference in the absence of catastrophic cancellation.

Verified thatr75454 makes the gamma test errors go away on a PPC Mac.

having the gamma function, one certainly needs the other Euler constant:C = 0.5772....the C constant is taken from GSL

for expm1, we use the Taylor development near zero, butwe have to precise what means small value for x. here this meansabs(x) < math.log(2.0).one can also use abs(x) < C

finally, for the error function, with the same kind of idea as with expm1,here is a pure python definitionthese patches are against r27a1:76674if all these make sense, i'll check for NAN, infinity, test suite...

may be some error when i send this the first time,i cannot see the file, so i resend mathmodule.c.diffsorry if you get it twice

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 identity  expm1(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) inr76755.

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#7518), I've added a new fileModules/_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 inr76861 (trunk) andr76863 (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 inr76879 (trunk),r76880 (py3k).

Oops.  That'sr76878 (trunk) andr76879 (py3k).

The last two functions to consider adding are exp2 and log2.  Does anyone care about these?

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