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

1 files changed, 19 insertions, 24 deletions

diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 7ace997..16d91e9 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -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)