Skip to content

Commit 664b511

Browse files
mdickinsonmdickinson
committed
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 python#3366: Add expm1 function to math module. Thanks Eric Smith for testing on Windows. ........
1 parent ef1992b commit 664b511

File tree

13 files changed

+162
-10
lines changed

13 files changed

+162
-10
lines changed

‎Doc/library/math.rst‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,20 @@ Power and logarithmic functions
148148
Return ``e**x``.
149149

150150

151+
.. function:: expm1(x)
152+
153+
Return ``e**x - 1``. For small floats *x*, the subtraction in
154+
``exp(x) - 1`` can result in a significant loss of precision; the
155+
:func:`expm1` function provides a way to compute this quantity to
156+
full precision::
157+
158+
>>> from math import exp, expm1
159+
>>> exp(1e-5) - 1 # gives result accurate to 11 places
160+
1.0000050000069649e-05
161+
>>> expm1(1e-5) # result accurate to full precision
162+
1.0000050000166668e-05
163+
164+
151165
.. function:: log(x[, base])
152166

153167
With one argument, return the natural logarithm of *x* (to base *e*).

‎Lib/test/math_testcases.txt‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -249,3 +249,73 @@ gam0132 gamma -4503599627370495.5 -> 0.0
249249
-- thanks to loss of accuracy in 1-x
250250
gam0140 gamma -63.349078729022985 -> 4.1777971677761880e-88
251251
gam0141 gamma -127.45117632943295 -> 1.1831110896236810e-214
252+
253+
-----------------------------------------------------------
254+
-- expm1: exp(x) - 1, without precision loss for small x --
255+
-----------------------------------------------------------
256+
257+
-- special values
258+
expm10000 expm1 0.0 -> 0.0
259+
expm10001 expm1 -0.0 -> -0.0
260+
expm10002 expm1 inf -> inf
261+
expm10003 expm1 -inf -> -1.0
262+
expm10004 expm1 nan -> nan
263+
264+
-- expm1(x) ~ x for tiny x
265+
expm10010 expm1 5e-324 -> 5e-324
266+
expm10011 expm1 1e-320 -> 1e-320
267+
expm10012 expm1 1e-300 -> 1e-300
268+
expm10013 expm1 1e-150 -> 1e-150
269+
expm10014 expm1 1e-20 -> 1e-20
270+
271+
expm10020 expm1 -5e-324 -> -5e-324
272+
expm10021 expm1 -1e-320 -> -1e-320
273+
expm10022 expm1 -1e-300 -> -1e-300
274+
expm10023 expm1 -1e-150 -> -1e-150
275+
expm10024 expm1 -1e-20 -> -1e-20
276+
277+
-- moderate sized values, where direct evaluation runs into trouble
278+
expm10100 expm1 1e-10 -> 1.0000000000500000e-10
279+
expm10101 expm1 -9.9999999999999995e-08 -> -9.9999995000000163e-8
280+
expm10102 expm1 3.0000000000000001e-05 -> 3.0000450004500034e-5
281+
expm10103 expm1 -0.0070000000000000001 -> -0.0069755570667648951
282+
expm10104 expm1 -0.071499208740094633 -> -0.069002985744820250
283+
expm10105 expm1 -0.063296004180116799 -> -0.061334416373633009
284+
expm10106 expm1 0.02390954035597756 -> 0.024197665143819942
285+
expm10107 expm1 0.085637352649044901 -> 0.089411184580357767
286+
expm10108 expm1 0.5966174947411006 -> 0.81596588596501485
287+
expm10109 expm1 0.30247206212075139 -> 0.35319987035848677
288+
expm10110 expm1 0.74574727375889516 -> 1.1080161116737459
289+
expm10111 expm1 0.97767512926555711 -> 1.6582689207372185
290+
expm10112 expm1 0.8450154566787712 -> 1.3280137976535897
291+
expm10113 expm1 -0.13979260323125264 -> -0.13046144381396060
292+
expm10114 expm1 -0.52899322039643271 -> -0.41080213643695923
293+
expm10115 expm1 -0.74083261478900631 -> -0.52328317124797097
294+
expm10116 expm1 -0.93847766984546055 -> -0.60877704724085946
295+
expm10117 expm1 10.0 -> 22025.465794806718
296+
expm10118 expm1 27.0 -> 532048240600.79865
297+
expm10119 expm1 123 -> 2.6195173187490626e+53
298+
expm10120 expm1 -12.0 -> -0.99999385578764666
299+
expm10121 expm1 -35.100000000000001 -> -0.99999999999999944
300+
301+
-- extreme negative values
302+
expm10201 expm1 -37.0 -> -0.99999999999999989
303+
expm10200 expm1 -38.0 -> -1.0
304+
expm10210 expm1 -710.0 -> -1.0
305+
-- the formula expm1(x) = 2 * sinh(x/2) * exp(x/2) doesn't work so
306+
-- well when exp(x/2) is subnormal or underflows to zero; check we're
307+
-- not using it!
308+
expm10211 expm1 -1420.0 -> -1.0
309+
expm10212 expm1 -1450.0 -> -1.0
310+
expm10213 expm1 -1500.0 -> -1.0
311+
expm10214 expm1 -1e50 -> -1.0
312+
expm10215 expm1 -1.79e308 -> -1.0
313+
314+
-- extreme positive values
315+
expm10300 expm1 300 -> 1.9424263952412558e+130
316+
expm10301 expm1 700 -> 1.0142320547350045e+304
317+
expm10302 expm1 709.78271289328393 -> 1.7976931346824240e+308
318+
expm10303 expm1 709.78271289348402 -> inf overflow
319+
expm10304 expm1 1000 -> inf overflow
320+
expm10305 expm1 1e50 -> inf overflow
321+
expm10306 expm1 1.79e308 -> inf overflow

‎Lib/test/test_math.py‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -984,17 +984,16 @@ def test_mtestfile(self):
984984
if math.isnan(expected) and math.isnan(got):
985985
continue
986986
if not math.isnan(expected) and not math.isnan(got):
987-
# we use different closeness criteria for
988-
# different functions.
989-
if fn == 'gamma':
990-
accuracy_failure = ulps_check(expected, got, 20)
991-
elif fn == 'lgamma':
987+
if fn == 'lgamma':
988+
# we use a weaker accuracy test for lgamma;
989+
# lgamma only achieves an absolute error of
990+
# a few multiples of the machine accuracy, in
991+
# general.
992992
accuracy_failure = acc_check(expected, got,
993993
rel_err = 5e-15,
994994
abs_err = 5e-15)
995995
else:
996-
raise ValueError("don't know how to check accuracy "
997-
"for this function")
996+
accuracy_failure = ulps_check(expected, got, 20)
998997
if accuracy_failure is None:
999998
continue
1000999

‎Misc/NEWS‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -457,7 +457,7 @@ Extension Modules
457457

458458
- Issue #7078: Set struct.__doc__ from _struct.__doc__.
459459

460-
- Issue #3366: Add gamma, lgamma functions to math module.
460+
- Issue #3366: Add expm1, gamma, lgamma functions to math module.
461461

462462
- Issue #6877: It is now possible to link the readline extension to the
463463
libedit readline emulation on OSX 10.5 or later.

‎Modules/Setup.dist‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -158,7 +158,7 @@ _symtable symtablemodule.c
158158

159159
#array arraymodule.c # array objects
160160
#cmath cmathmodule.c # -lm # complex math library functions
161-
#math mathmodule.c # -lm # math library functions, e.g. sin()
161+
#math mathmodule.c _math.c # -lm # math library functions, e.g. sin()
162162
#_struct _struct.c # binary structure packing/unpacking
163163
#time timemodule.c # -lm # time operations and variables
164164
#operator operator.c # operator.add() and similar goodies

‎Modules/_math.c‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
/* Definitions of some C99 math library functions, for those platforms
2+
that don't implement these functions already. */
3+
4+
#include <float.h>
5+
#include <math.h>
6+
7+
/* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed
8+
to avoid the significant loss of precision that arises from direct
9+
evaluation of the expression exp(x) - 1, for x near 0. */
10+
11+
double
12+
_Py_expm1(double x)
13+
{
14+
/* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this
15+
also works fine for infinities and nans.
16+
17+
For smaller x, we can use a method due to Kahan that achieves close to
18+
full accuracy.
19+
*/
20+
21+
if (fabs(x) < 0.7) {
22+
double u;
23+
u = exp(x);
24+
if (u == 1.0)
25+
return x;
26+
else
27+
return (u - 1.0) * x / log(u);
28+
}
29+
else
30+
return exp(x) - 1.0;
31+
}

‎Modules/_math.h‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
double _Py_expm1(double x);
2+
3+
#ifdef HAVE_EXPM1
4+
#define m_expm1 expm1
5+
#else
6+
/* if the system doesn't have expm1, use the substitute
7+
function defined in Modules/_math.c. */
8+
#define m_expm1 _Py_expm1
9+
#endif

‎Modules/mathmodule.c‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@ raised for division by zero and mod by zero.
5353
*/
5454

5555
#include "Python.h"
56+
#include "_math.h"
5657
#include "longintrepr.h" /* just for SHIFT */
5758

5859
#ifdef _OSF_SOURCE
@@ -722,6 +723,10 @@ FUNC1(cosh, cosh, 1,
722723
"cosh(x)\n\nReturn the hyperbolic cosine of x.")
723724
FUNC1(exp, exp, 1,
724725
"exp(x)\n\nReturn e raised to the power of x.")
726+
FUNC1(expm1, m_expm1, 1,
727+
"expm1(x)\n\nReturn exp(x)-1.\n"
728+
"This function avoids the loss of precision involved in the direct "
729+
"evaluation of exp(x)-1 for small x.")
725730
FUNC1(fabs, fabs, 0,
726731
"fabs(x)\n\nReturn the absolute value of the float x.")
727732

@@ -1493,6 +1498,7 @@ static PyMethodDef math_methods[] = {
14931498
{"cosh", math_cosh, METH_O, math_cosh_doc},
14941499
{"degrees", math_degrees, METH_O, math_degrees_doc},
14951500
{"exp", math_exp, METH_O, math_exp_doc},
1501+
{"expm1", math_expm1, METH_O, math_expm1_doc},
14961502
{"fabs", math_fabs, METH_O, math_fabs_doc},
14971503
{"factorial", math_factorial, METH_O, math_factorial_doc},
14981504
{"floor", math_floor, METH_O, math_floor_doc},

‎PC/VC6/pythoncore.dsp‎

Lines changed: 4 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎PC/VS7.1/pythoncore.vcproj‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -408,6 +408,9 @@
408408
<File
409409
RelativePath="..\..\Modules\_lsprof.c">
410410
</File>
411+
<File
412+
RelativePath="..\..\Modules\_math.c">
413+
</File>
411414
<File
412415
RelativePath="..\..\Modules\_pickle.c">
413416
</File>

0 commit comments

Comments
 (0)