-
-
Notifications
You must be signed in to change notification settings - Fork 31.1k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add gamma function, error functions and other C99 math.h functions to math module #47616
Comments
From Py3000 list: in C99 standard math library, but need code when not |
here is an attempt to make erf, erfc, lgamma and tgamma i tested only on gcc-4.3.1/linux-2.6.26. some results: Python 3.0b2 (r30b2:65080, Jul 21 2008, 13:13:13) >>> 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 |
and this is the header. |
Thanks for the patch! I probably won't get to this properly until after By the way, there's already a setup in Python 2.6 (ad 3.0) for substitutes So a full patch for this should touch at least Python/pymath.c, |
Since we're not in a hurry for Py2.7 and 3.1, I would like to this There are several different approximations to choose from. Each of Let's let the community use its collective wisdom to select the best At one point, Tim was reluctant to add any of these functions because FWIW, I changed the 2.6 version of test_random.gamma() to take |
On Mon, Jul 21, 2008 at 5:34 PM, Raymond Hettinger
FWIW, the current patch dynamically selects one of a few different +1 on kicking it around and comparing the output with that from other |
It would be nice if we knew the error bounds for each of the When switching from one method to another, it might be nice to have a The lowword/highword macros look to be tightly tied to the internal |
Mark Dickinson :
Raymond Hettinger:
|
Why isn't tgamma() simply named gamma()? The t prefix does nothing for |
I think he just carried the names over from C, where: tgamma() is the "true gamma function" and gamma() might be tgamma() or lgamma() depending on which C library I'm +1 on making gamma() be the true gamma function and not carrying |
+1 from me, too. |
ouch! if you want to test the erf, erfc, lgamma and gamma the other diffs are incomplete and don't work, unless you have |
Can you also implement blending of approximations: (1-t)*f1(x) + t*f2 |
On Fri, Jul 25, 2008 at 1:37 PM, Raymond Hettinger
Is this necessary? Are the approximations significantly different near the (Haven't had a chance to download the patch and try it myself as a I'm |
i wrote pure Python equivalent of the patch and will Raymond Hettinger:
i do this with the Python code for now, and my question is: with erf, for example, one has typically 5 domains 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, 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. for erf functions, blending of approximations, as i understand it, |
pure Python codes are posted to ASPN: http://code.activestate.com/recipes/576391/ i also rewrote the patch so that they don't |
There have been requests to add other C99 libm functions to Python: With these two, and the gamma and error functions, we're pretty close to |
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 code has gone in a new file called Modules/math_cextra.c. The old There are currently many test failures; some of these are undoubtedly To test the code on non-Windows, build with something like: CC='gcc -DTEST_GAMMA' ./configure && make the '-DTEST_GAMMA' just forces the build to use nirinA's code; without This is just work-in-progress; apart from the test failures, there's The patch is against the trunk. |
Here's a more careful patch for just the gamma function. It's a fairly I'd appreciate any feedback, especially if anyone has the chance to test This patch *doesn't* use the system gamma function, even when available. Once this is working, I'll concentrate on lgamma. Then erf and erfc. |
New patch for gamma , with some tweaks:
On my machine, testing with approx. 10**7 random samples, this version I plan to check this in in a week or two. Feedback welcome! It would be |
I'll test your patch on Windows. Are you working against the trunk or |
Thanks! The patch is against the trunk. (It doesn't quite apply cleanly Hmm. Rereading my previous comment, I seem to have a blindness for
should have been
and
should have been
|
I'm only setup to test the official Windows build setup under PCbuild. The patch against the trunk failed for PC/VC6/pythoncore.dsp. I don't Your patch built fine, with one warning about a loss of precision on the n = round(2.0*y); I recommend explicitly casting the result of round() to int, to get rid I ran the following tests, all of which passed: PCbuild/python.exe Lib/test/regrtest.py -v test_math I appreciate your work on this patch. Let me know if there's anything |
Many thanks, Daniel.
I've no idea why that would happen. A line-ending issue, perhaps? If
Done. Thanks! |
gamma5.patch. very minor changes from the last patch:
|
Gamma function added in r75117 (trunk), r75119 (py3k). |
A first attempt for lgamma, using Lanczos' formula. In general, this gives very good accuracy. As the tests show, there's one I think it's reasonable to allow an absolute error of a few times the |
I agree with your reasoning. Any chance of also getting the incomplete beta and gamma functions, |
FYI, approximately 20 of the gamma test cases fail on PPC Macs. Attached |
FYI, mysterious numeric differences on PPC are often due to the C |
For the record, these OS X installer build scripts (running on 10.5) gcc-4.0 -arch ppc -arch i386 -isysroot /Developer/SDKs/MacOSX10.4u.sdk - and: $ gcc-4.0 --version
powerpc-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 I'm suspecting an endianness issue �and/or a struct module bug. Ned, do you still get these failures with a direct single-architecture |
Adding -mno-fused-madd would be worth trying. It usually fixes PPC bugs ;-) |
It was a stupid Mark bug. I'd missed out a '<' from the struct.unpack Thanks for reporting this, Ned. |
BTW, the gamma computation shouldn't be at all fragile, unless I've messed |
Mark, you needn't bother: you found the smoking gun already! From your |
Verified that r75454 makes the gamma test errors go away on a PPC Mac. |
having the gamma function, one certainly needs the other Euler constant: the C constant is taken from GSL |
for expm1, we use the Taylor development near zero, but |
finally, for the error function, with the same kind of idea as with expm1, these patches are against r27a1:76674 if all these make sense, i'll check for NAN, infinity, test suite... |
may be some error when i send this the first time, |
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 For expm1, I was planning to use the libm function when it exists (which 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 |
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 Here's a prototype patch for the erf and erfc functions that uses a power |
Here's the C version. |
Here's a patch to add expm1. Rather than putting the code in pymath.c, I'm fairly confident that the maths and numerics are okay, but less |
I committed the patch for expm1 in r76861 (trunk) and r76863 (py3k), with 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 |
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? |
Not much, I suspect. :) I'd expect (but am mostly guessing) exp2(x) to have better accuracy than Similarly for log2: log2(n) should be a touch more accurate than But we've already got the 'bit_length' method for integers, which fills |
Will close this unless there's an outcry of support for exp2 and log2. |
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:
The text was updated successfully, but these errors were encountered: