From 664b511c0acdfdecdec92d2255ffd94c4e6d5f7a Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Wed, 16 Dec 2009 20:23:42 +0000 Subject: [PATCH] Merged revisions 76861 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r76861 | mark.dickinson | 2009-12-16 20:13:40 +0000 (Wed, 16 Dec 2009) | 3 lines Issue #3366: Add expm1 function to math module. Thanks Eric Smith for testing on Windows. ........ --- Doc/library/math.rst | 14 ++++++++ Lib/test/math_testcases.txt | 70 +++++++++++++++++++++++++++++++++++++ Lib/test/test_math.py | 13 ++++--- Misc/NEWS | 2 +- Modules/Setup.dist | 2 +- Modules/_math.c | 31 ++++++++++++++++ Modules/_math.h | 9 +++++ Modules/mathmodule.c | 6 ++++ PC/VC6/pythoncore.dsp | 4 +++ PC/VS7.1/pythoncore.vcproj | 3 ++ PC/VS8.0/pythoncore.vcproj | 8 +++++ PCbuild/pythoncore.vcproj | 8 +++++ setup.py | 2 +- 13 files changed, 162 insertions(+), 10 deletions(-) create mode 100644 Modules/_math.c create mode 100644 Modules/_math.h diff --git a/Doc/library/math.rst b/Doc/library/math.rst index e903c5f9022..39aea2903ec 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -148,6 +148,20 @@ Power and logarithmic functions Return ``e**x``. +.. function:: expm1(x) + + Return ``e**x - 1``. For small floats *x*, the subtraction in + ``exp(x) - 1`` can result in a significant loss of precision; the + :func:`expm1` function provides a way to compute this quantity to + full precision:: + + >>> from math import exp, expm1 + >>> exp(1e-5) - 1 # gives result accurate to 11 places + 1.0000050000069649e-05 + >>> expm1(1e-5) # result accurate to full precision + 1.0000050000166668e-05 + + .. function:: log(x[, base]) With one argument, return the natural logarithm of *x* (to base *e*). diff --git a/Lib/test/math_testcases.txt b/Lib/test/math_testcases.txt index 2b1fc694700..bb9929054f9 100644 --- a/Lib/test/math_testcases.txt +++ b/Lib/test/math_testcases.txt @@ -249,3 +249,73 @@ gam0132 gamma -4503599627370495.5 -> 0.0 -- thanks to loss of accuracy in 1-x gam0140 gamma -63.349078729022985 -> 4.1777971677761880e-88 gam0141 gamma -127.45117632943295 -> 1.1831110896236810e-214 + +----------------------------------------------------------- +-- expm1: exp(x) - 1, without precision loss for small x -- +----------------------------------------------------------- + +-- special values +expm10000 expm1 0.0 -> 0.0 +expm10001 expm1 -0.0 -> -0.0 +expm10002 expm1 inf -> inf +expm10003 expm1 -inf -> -1.0 +expm10004 expm1 nan -> nan + +-- expm1(x) ~ x for tiny x +expm10010 expm1 5e-324 -> 5e-324 +expm10011 expm1 1e-320 -> 1e-320 +expm10012 expm1 1e-300 -> 1e-300 +expm10013 expm1 1e-150 -> 1e-150 +expm10014 expm1 1e-20 -> 1e-20 + +expm10020 expm1 -5e-324 -> -5e-324 +expm10021 expm1 -1e-320 -> -1e-320 +expm10022 expm1 -1e-300 -> -1e-300 +expm10023 expm1 -1e-150 -> -1e-150 +expm10024 expm1 -1e-20 -> -1e-20 + +-- moderate sized values, where direct evaluation runs into trouble +expm10100 expm1 1e-10 -> 1.0000000000500000e-10 +expm10101 expm1 -9.9999999999999995e-08 -> -9.9999995000000163e-8 +expm10102 expm1 3.0000000000000001e-05 -> 3.0000450004500034e-5 +expm10103 expm1 -0.0070000000000000001 -> -0.0069755570667648951 +expm10104 expm1 -0.071499208740094633 -> -0.069002985744820250 +expm10105 expm1 -0.063296004180116799 -> -0.061334416373633009 +expm10106 expm1 0.02390954035597756 -> 0.024197665143819942 +expm10107 expm1 0.085637352649044901 -> 0.089411184580357767 +expm10108 expm1 0.5966174947411006 -> 0.81596588596501485 +expm10109 expm1 0.30247206212075139 -> 0.35319987035848677 +expm10110 expm1 0.74574727375889516 -> 1.1080161116737459 +expm10111 expm1 0.97767512926555711 -> 1.6582689207372185 +expm10112 expm1 0.8450154566787712 -> 1.3280137976535897 +expm10113 expm1 -0.13979260323125264 -> -0.13046144381396060 +expm10114 expm1 -0.52899322039643271 -> -0.41080213643695923 +expm10115 expm1 -0.74083261478900631 -> -0.52328317124797097 +expm10116 expm1 -0.93847766984546055 -> -0.60877704724085946 +expm10117 expm1 10.0 -> 22025.465794806718 +expm10118 expm1 27.0 -> 532048240600.79865 +expm10119 expm1 123 -> 2.6195173187490626e+53 +expm10120 expm1 -12.0 -> -0.99999385578764666 +expm10121 expm1 -35.100000000000001 -> -0.99999999999999944 + +-- extreme negative values +expm10201 expm1 -37.0 -> -0.99999999999999989 +expm10200 expm1 -38.0 -> -1.0 +expm10210 expm1 -710.0 -> -1.0 +-- the formula expm1(x) = 2 * sinh(x/2) * exp(x/2) doesn't work so +-- well when exp(x/2) is subnormal or underflows to zero; check we're +-- not using it! +expm10211 expm1 -1420.0 -> -1.0 +expm10212 expm1 -1450.0 -> -1.0 +expm10213 expm1 -1500.0 -> -1.0 +expm10214 expm1 -1e50 -> -1.0 +expm10215 expm1 -1.79e308 -> -1.0 + +-- extreme positive values +expm10300 expm1 300 -> 1.9424263952412558e+130 +expm10301 expm1 700 -> 1.0142320547350045e+304 +expm10302 expm1 709.78271289328393 -> 1.7976931346824240e+308 +expm10303 expm1 709.78271289348402 -> inf overflow +expm10304 expm1 1000 -> inf overflow +expm10305 expm1 1e50 -> inf overflow +expm10306 expm1 1.79e308 -> inf overflow diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 20524300212..5ae5d3a9b19 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -984,17 +984,16 @@ class MathTests(unittest.TestCase): if math.isnan(expected) and math.isnan(got): continue if not math.isnan(expected) and not math.isnan(got): - # we use different closeness criteria for - # different functions. - if fn == 'gamma': - accuracy_failure = ulps_check(expected, got, 20) - elif fn == 'lgamma': + if fn == 'lgamma': + # we use a weaker accuracy test for lgamma; + # lgamma only achieves an absolute error of + # a few multiples of the machine accuracy, in + # general. accuracy_failure = acc_check(expected, got, rel_err = 5e-15, abs_err = 5e-15) else: - raise ValueError("don't know how to check accuracy " - "for this function") + accuracy_failure = ulps_check(expected, got, 20) if accuracy_failure is None: continue diff --git a/Misc/NEWS b/Misc/NEWS index 4ad37731d2c..d4121d1ee0f 100644 --- a/Misc/NEWS +++ b/Misc/NEWS @@ -457,7 +457,7 @@ Extension Modules - Issue #7078: Set struct.__doc__ from _struct.__doc__. -- Issue #3366: Add gamma, lgamma functions to math module. +- Issue #3366: Add expm1, gamma, lgamma functions to math module. - Issue #6877: It is now possible to link the readline extension to the libedit readline emulation on OSX 10.5 or later. diff --git a/Modules/Setup.dist b/Modules/Setup.dist index 9f189deb4c0..e1cb235ffa6 100644 --- a/Modules/Setup.dist +++ b/Modules/Setup.dist @@ -158,7 +158,7 @@ _symtable symtablemodule.c #array arraymodule.c # array objects #cmath cmathmodule.c # -lm # complex math library functions -#math mathmodule.c # -lm # math library functions, e.g. sin() +#math mathmodule.c _math.c # -lm # math library functions, e.g. sin() #_struct _struct.c # binary structure packing/unpacking #time timemodule.c # -lm # time operations and variables #operator operator.c # operator.add() and similar goodies diff --git a/Modules/_math.c b/Modules/_math.c new file mode 100644 index 00000000000..9d330aa0d59 --- /dev/null +++ b/Modules/_math.c @@ -0,0 +1,31 @@ +/* Definitions of some C99 math library functions, for those platforms + that don't implement these functions already. */ + +#include +#include + +/* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed + to avoid the significant loss of precision that arises from direct + evaluation of the expression exp(x) - 1, for x near 0. */ + +double +_Py_expm1(double x) +{ + /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this + also works fine for infinities and nans. + + For smaller x, we can use a method due to Kahan that achieves close to + full accuracy. + */ + + if (fabs(x) < 0.7) { + double u; + u = exp(x); + if (u == 1.0) + return x; + else + return (u - 1.0) * x / log(u); + } + else + return exp(x) - 1.0; +} diff --git a/Modules/_math.h b/Modules/_math.h new file mode 100644 index 00000000000..69c96b5ab76 --- /dev/null +++ b/Modules/_math.h @@ -0,0 +1,9 @@ +double _Py_expm1(double x); + +#ifdef HAVE_EXPM1 +#define m_expm1 expm1 +#else +/* if the system doesn't have expm1, use the substitute + function defined in Modules/_math.c. */ +#define m_expm1 _Py_expm1 +#endif diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index ece68a742c8..ef842989ba8 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -53,6 +53,7 @@ raised for division by zero and mod by zero. */ #include "Python.h" +#include "_math.h" #include "longintrepr.h" /* just for SHIFT */ #ifdef _OSF_SOURCE @@ -722,6 +723,10 @@ FUNC1(cosh, cosh, 1, "cosh(x)\n\nReturn the hyperbolic cosine of x.") FUNC1(exp, exp, 1, "exp(x)\n\nReturn e raised to the power of x.") +FUNC1(expm1, m_expm1, 1, + "expm1(x)\n\nReturn exp(x)-1.\n" + "This function avoids the loss of precision involved in the direct " + "evaluation of exp(x)-1 for small x.") FUNC1(fabs, fabs, 0, "fabs(x)\n\nReturn the absolute value of the float x.") @@ -1493,6 +1498,7 @@ static PyMethodDef math_methods[] = { {"cosh", math_cosh, METH_O, math_cosh_doc}, {"degrees", math_degrees, METH_O, math_degrees_doc}, {"exp", math_exp, METH_O, math_exp_doc}, + {"expm1", math_expm1, METH_O, math_expm1_doc}, {"fabs", math_fabs, METH_O, math_fabs_doc}, {"factorial", math_factorial, METH_O, math_factorial_doc}, {"floor", math_floor, METH_O, math_floor_doc}, diff --git a/PC/VC6/pythoncore.dsp b/PC/VC6/pythoncore.dsp index ee7e36628c6..4e598ffa28d 100644 --- a/PC/VC6/pythoncore.dsp +++ b/PC/VC6/pythoncore.dsp @@ -157,6 +157,10 @@ SOURCE=..\..\Modules\_lsprof.c # End Source File # Begin Source File +SOURCE=..\..\Modules\_math.c +# End Source File +# Begin Source File + SOURCE=..\..\Modules\_pickle.c # End Source File # Begin Source File diff --git a/PC/VS7.1/pythoncore.vcproj b/PC/VS7.1/pythoncore.vcproj index 214a690c998..b5b4f582305 100644 --- a/PC/VS7.1/pythoncore.vcproj +++ b/PC/VS7.1/pythoncore.vcproj @@ -408,6 +408,9 @@ + + diff --git a/PC/VS8.0/pythoncore.vcproj b/PC/VS8.0/pythoncore.vcproj index 726e6301cad..22850138289 100644 --- a/PC/VS8.0/pythoncore.vcproj +++ b/PC/VS8.0/pythoncore.vcproj @@ -1042,6 +1042,14 @@ RelativePath="..\..\Modules\_lsprof.c" > + + + + diff --git a/PCbuild/pythoncore.vcproj b/PCbuild/pythoncore.vcproj index 650d244b2f9..21594cd6877 100644 --- a/PCbuild/pythoncore.vcproj +++ b/PCbuild/pythoncore.vcproj @@ -1011,6 +1011,14 @@ RelativePath="..\Modules\_lsprof.c" > + + + + diff --git a/setup.py b/setup.py index 43b6a6a6557..3666e8fbb74 100644 --- a/setup.py +++ b/setup.py @@ -398,7 +398,7 @@ class PyBuildExt(build_ext): libraries=math_libs) ) # math library functions, e.g. sin() - exts.append( Extension('math', ['mathmodule.c'], + exts.append( Extension('math', ['mathmodule.c', '_math.c'], libraries=math_libs) ) # time operations and variables exts.append( Extension('time', ['timemodule.c'],