Issue #3366: Add expm1 function to math module. Thanks Eric Smith for

testing on Windows.

4 files changed, 47 insertions, 1 deletions

diff --git a/Modules/Setup.dist b/Modules/Setup.dist
index 5e69a60..3a6c973 100644
--- a/Modules/Setup.dist
+++ b/Modules/Setup.dist
@@ -169,7 +169,7 @@ GLHACK=-Dclear=__GLclear
 
 #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 0000000..9d330aa
--- /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 <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
index 0000000..69c96b5
--- /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 b1e4521..469fdf8 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
@@ -686,6 +687,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.")
 FUNC1(floor, floor, 0,
@@ -1420,6 +1425,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},