]> granicus.if.org Git - python/commitdiff
Merged revisions 76861 via svnmerge from
authorMark Dickinson <dickinsm@gmail.com>
Wed, 16 Dec 2009 20:23:42 +0000 (20:23 +0000)
committerMark Dickinson <dickinsm@gmail.com>
Wed, 16 Dec 2009 20:23:42 +0000 (20:23 +0000)
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.
........

13 files changed:
Doc/library/math.rst
Lib/test/math_testcases.txt
Lib/test/test_math.py
Misc/NEWS
Modules/Setup.dist
Modules/_math.c [new file with mode: 0644]
Modules/_math.h [new file with mode: 0644]
Modules/mathmodule.c
PC/VC6/pythoncore.dsp
PC/VS7.1/pythoncore.vcproj
PC/VS8.0/pythoncore.vcproj
PCbuild/pythoncore.vcproj
setup.py

index e903c5f9022210990b6503bff06cdcb67edbe827..39aea2903ec3d1a726d67ad399e74ae64efc3e4f 100644 (file)
@@ -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*).
index 2b1fc6947005991352ae75a1ca1ed107bb74ef09..bb9929054f96953ad88651483e94617a3450db85 100644 (file)
@@ -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
index 20524300212ebd8b61a6984c4451a9be25e054e6..5ae5d3a9b193607bf0bbf8cea71741cd35b9a477 100644 (file)
@@ -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
 
index 4ad37731d2c8bfb79aedb775be878156cb546eb1..d4121d1ee0f2d0addddbd909d44861bb59ebdc6a 100644 (file)
--- 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.
index 9f189deb4c077520f395f428d72ed49513e5e031..e1cb235ffa66dd54c16ee4f26cae414c24246abc 100644 (file)
@@ -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 (file)
index 0000000..9d330aa
--- /dev/null
@@ -0,0 +1,31 @@
+/* Definitions of some C99 math library functions, for those platforms
+   that don't implement these functions already. */
+
+#include <float.h>
+#include <math.h>
+
+/* 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 (file)
index 0000000..69c96b5
--- /dev/null
@@ -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
index ece68a742c803eb56d79bc0c823fd50351a93672..ef842989ba888ec8dfbe1f178c08e0cd3db83e91 100644 (file)
@@ -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},
index ee7e36628c69641aa85b9d207a448bdf3b14fdb5..4e598ffa28d79a8690d994dd02e3310defb35fa1 100644 (file)
@@ -157,6 +157,10 @@ SOURCE=..\..\Modules\_lsprof.c
 # End Source File\r
 # Begin Source File\r
 \r
+SOURCE=..\..\Modules\_math.c\r
+# End Source File\r
+# Begin Source File\r
+\r
 SOURCE=..\..\Modules\_pickle.c\r
 # End Source File\r
 # Begin Source File\r
index 214a690c998cef926a17f6d21cee9ec6877422d1..b5b4f582305ad5af386cbbdc7a2129f69f65c978 100644 (file)
                <File
                        RelativePath="..\..\Modules\_lsprof.c">
                </File>
+               <File
+                       RelativePath="..\..\Modules\_math.c">
+               </File>
                <File
                        RelativePath="..\..\Modules\_pickle.c">
                </File>
index 726e6301cadcca0d078d64c6b20ab1ec4bf82e47..22850138289a3a9760e3af4e0e3254c9e17ed134 100644 (file)
                                RelativePath="..\..\Modules\_lsprof.c"
                                >
                        </File>
+                       <File
+                               RelativePath="..\..\Modules\_math.c"
+                               >
+                       </File>
+                       <File
+                               RelativePath="..\..\Modules\_math.h"
+                               >
+                       </File>
                        <File
                                RelativePath="..\..\Modules\_pickle.c"
                                >
index 650d244b2f9344073c20f08fffb1b12cee0c3a3b..21594cd687786ea2d76c91a4e54340da21052b96 100644 (file)
                                RelativePath="..\Modules\_lsprof.c"
                                >
                        </File>
+                       <File
+                               RelativePath="..\Modules\_math.c"
+                               >
+                       </File>
+                       <File
+                               RelativePath="..\Modules\_math.h"
+                               >
+                       </File>
                        <File
                                RelativePath="..\Modules\_pickle.c"
                                >
index 43b6a6a6557b2fef44e5e02cf7414e6bb7a3e5a4..3666e8fbb743892a3b7639209f1f7305e73c55ac 100644 (file)
--- 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'],