82
82
* should return a zero on underflow, and +- HUGE_VAL on
83
83
* overflow, so testing the result for zero suffices to
84
84
* distinguish the cases).
86
* On some platforms (Ubuntu/ia64) it seems that errno can be
87
* set to ERANGE for subnormal results that do *not* underflow
88
* to zero. So to be safe, we'll ignore ERANGE whenever the
89
* function result is less than one in absolute value.
87
94
PyErr_SetString(PyExc_OverflowError,
88
95
"math range error");
93
98
/* Unexpected math error */
322
327
sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
323
328
overflow of the first partial sum.
325
Note 3: The itermediate values lo, yr, and hi are declared volatile so
326
aggressive compilers won't algebraicly reduce lo to always be exactly 0.0.
330
Note 3: The intermediate values lo, yr, and hi are declared volatile so
331
aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
327
332
Also, the volatile declaration forces the values to be stored in memory as
328
333
regular doubles instead of extended long precision (80-bit) values. This
329
prevents double rounding because any addition or substraction of two doubles
334
prevents double rounding because any addition or subtraction of two doubles
330
335
can be resolved exactly into double-sized hi and lo values. As long as the
331
336
hi value gets forced into a double before yr and lo are computed, the extra
332
337
bits in downstream extended precision operations (x87 for example) will be
336
341
Note 4: A similar implementation is in Modules/cmathmodule.c.
337
342
Be sure to update both when making changes.
339
Note 5: The signature of math.sum() differs from __builtin__.sum()
344
Note 5: The signature of math.fsum() differs from __builtin__.sum()
340
345
because the start argument doesn't make sense in the context of
341
346
accurate summation. Since the partials table is collapsed before
342
347
returning a result, sum(seq2, start=sum(seq1)) may not equal the
348
353
/* Extend the partials array p[] by doubling its size. */
349
354
static int /* non-zero on error */
350
_sum_realloc(double **p_ptr, Py_ssize_t n,
355
_fsum_realloc(double **p_ptr, Py_ssize_t n,
351
356
double *ps, Py_ssize_t *m_ptr)
407
math_sum(PyObject *self, PyObject *seq)
412
math_fsum(PyObject *self, PyObject *seq)
409
414
PyObject *item, *iter, *sum = NULL;
410
415
Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
411
416
double x, y, t, ps[NUM_PARTIALS], *p = ps;
417
double xsave, special_sum = 0.0, inf_sum = 0.0;
412
418
volatile double hi, yr, lo;
414
420
iter = PyObject_GetIter(seq);
415
421
if (iter == NULL)
418
PyFPE_START_PROTECT("sum", Py_DECREF(iter); return NULL)
424
PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
420
426
for(;;) { /* for x in iterable */
421
427
assert(0 <= n && n <= m);
449
n = i; /* ps[i:] = [x] */
456
n = i; /* ps[i:] = [x] */
451
/* If non-finite, reset partials, effectively
452
adding subsequent items without roundoff
453
and yielding correct non-finite results,
454
provided IEEE 754 rules are observed */
455
if (! Py_IS_FINITE(x))
458
if (! Py_IS_FINITE(x)) {
459
/* a nonfinite x could arise either as
460
a result of intermediate overflow, or
461
as a result of a nan or inf in the
463
if (Py_IS_FINITE(xsave)) {
464
PyErr_SetString(PyExc_OverflowError,
465
"intermediate overflow in fsum");
468
if (Py_IS_INFINITY(xsave))
470
special_sum += xsave;
457
else if (n >= m && _sum_realloc(&p, n, ps, &m))
474
else if (n >= m && _fsum_realloc(&p, n, ps, &m))
481
if (special_sum != 0.0) {
482
if (Py_IS_NAN(inf_sum))
483
PyErr_SetString(PyExc_ValueError,
484
"-inf + inf in fsum");
486
sum = PyFloat_FromDouble(special_sum);
466
if (Py_IS_FINITE(hi)) {
467
/* sum_exact(ps, hi) from the top, stop when the sum becomes inexact. */
471
assert(fabs(y) < fabs(x));
478
/* Make half-even rounding work across multiple partials. Needed
479
so that sum([1e-16, 1, 1e16]) will round-up the last digit to
480
two instead of down to zero (the 1e-16 makes the 1 slightly
481
closer to two). With a potential 1 ULP rounding error fixed-up,
482
math.sum() can guarantee commutativity. */
483
if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
484
(lo > 0.0 && p[n-1] > 0.0))) {
493
/* sum_exact(ps, hi) from the top, stop when the sum becomes
498
assert(fabs(y) < fabs(x));
492
else { /* raise exception corresponding to a special value */
493
errno = Py_IS_NAN(hi) ? EDOM : ERANGE;
505
/* Make half-even rounding work across multiple partials.
506
Needed so that sum([1e-16, 1, 1e16]) will round-up the last
507
digit to two instead of down to zero (the 1e-16 makes the 1
508
slightly closer to two). With a potential 1 ULP rounding
509
error fixed-up, math.fsum() can guarantee commutativity. */
510
if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
511
(lo > 0.0 && p[n-1] > 0.0))) {
498
519
sum = PyFloat_FromDouble(hi);
501
522
PyFPE_END_PROTECT(hi)
1002
1021
{"floor", math_floor, METH_O, math_floor_doc},
1003
1022
{"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
1004
1023
{"frexp", math_frexp, METH_O, math_frexp_doc},
1024
{"fsum", math_fsum, METH_O, math_fsum_doc},
1005
1025
{"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
1006
1026
{"isinf", math_isinf, METH_O, math_isinf_doc},
1007
1027
{"isnan", math_isnan, METH_O, math_isnan_doc},
1015
1035
{"sin", math_sin, METH_O, math_sin_doc},
1016
1036
{"sinh", math_sinh, METH_O, math_sinh_doc},
1017
1037
{"sqrt", math_sqrt, METH_O, math_sqrt_doc},
1018
{"sum", math_sum, METH_O, math_sum_doc},
1019
1038
{"tan", math_tan, METH_O, math_tan_doc},
1020
1039
{"tanh", math_tanh, METH_O, math_tanh_doc},
1021
{"trunc", math_trunc, METH_O, math_trunc_doc},
1040
{"trunc", math_trunc, METH_O, math_trunc_doc},
1022
1041
{NULL, NULL} /* sentinel */
added
removed
Modules/mathmodule.c
