Mirror
/
cpython
mirror of https://github.com/python/cpython
4
1
Fork 0
Code Issues Releases Wiki Activity

Address double-rounding scenarios by setting all variables to long doubles.

This commit is contained in:
Raymond Hettinger 2008-06-09 09:29:17 +00:00
parent ee4bcad68e
commit 7b1ed66372
1 changed files with 18 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
Modules/mathmodule.c
View File

@ -324,17 +324,12 @@ FUNC1(tanh, tanh, 0,
Note 3: The itermediate values lo, yr, and hi are declared volatile so
aggressive compilers won't algebraicly reduce lo to always be exactly 0.0.
Also, the volatile declaration forces the values to be stored in memory as
regular doubles instead of extended long precision (80-bit) values. This
prevents double rounding because any addition or substraction of two doubles
can be resolved exactly into double-sized hi and lo values. As long as the
hi value gets forced into a double before yr and lo are computed, the extra
bits in downstream extended precision operations (x87 for example) will be
exactly zero and therefore can be losslessly stored back into a double,
thereby preventing double rounding.
Note 4: A similar implementation is in Modules/cmathmodule.c.
Be sure to update both when making changes.
Note 4: Intermediate values and partial sums are declared as long doubles
as a way to eliminate double rounding environments where the operations
are carried-out in extended precision but stored in double precision
variables. In some cases, this doesn't help because the compiler
treats long doubles as doubles (i.e. the MS compiler for Win32 builds).
Note 5: The signature of math.sum() differs from __builtin__.sum()
because the start argument doesn't make sense in the context of
@ -347,28 +342,28 @@ FUNC1(tanh, tanh, 0,
/* Extend the partials array p[] by doubling its size. */
static int /* non-zero on error */
_sum_realloc(double **p_ptr, Py_ssize_t n,
double *ps, Py_ssize_t *m_ptr)
_sum_realloc(long double **p_ptr, Py_ssize_t n,
long double *ps, Py_ssize_t *m_ptr)
{
void *v = NULL;
Py_ssize_t m = *m_ptr;
m += m; /* double */
if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
double *p = *p_ptr;
m += m; /* long double */
if (n < m && m < (PY_SSIZE_T_MAX / sizeof(long double))) {
long double *p = *p_ptr;
if (p == ps) {
v = PyMem_Malloc(sizeof(double) * m);
v = PyMem_Malloc(sizeof(long double) * m);
if (v != NULL)
memcpy(v, ps, sizeof(double) * n);
memcpy(v, ps, sizeof(long double) * n);
}
else
v = PyMem_Realloc(p, sizeof(double) * m);
v = PyMem_Realloc(p, sizeof(long double) * m);
}
if (v == NULL) { /* size overflow or no memory */
PyErr_SetString(PyExc_MemoryError, "math sum partials");
return 1;
}
*p_ptr = (double*) v;
*p_ptr = (long double*) v;
*m_ptr = m;
return 0;
}
@ -408,8 +403,8 @@ math_sum(PyObject *self, PyObject *seq)
{
PyObject *item, *iter, *sum = NULL;
Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
double x, y, t, ps[NUM_PARTIALS], *p = ps;
volatile double hi, yr, lo;
long double x, y, t, ps[NUM_PARTIALS], *p = ps;
volatile long double hi, yr, lo;
iter = PyObject_GetIter(seq);
if (iter == NULL)
@ -428,7 +423,7 @@ math_sum(PyObject *self, PyObject *seq)
goto _sum_error;
break;
}
x = PyFloat_AsDouble(item);
x = (long double)PyFloat_AsDouble(item);
Py_DECREF(item);
if (PyErr_Occurred())
goto _sum_error;
@ -495,7 +490,7 @@ math_sum(PyObject *self, PyObject *seq)
goto _sum_error;
}
}
sum = PyFloat_FromDouble(hi);
sum = PyFloat_FromDouble((double)hi);
_sum_error:
PyFPE_END_PROTECT(hi)