/* the shared function declaration in floatobject.h??? */

extern int
float_sum(PyObject *seq, PyObject *sum, double *sum_ptr,
          double (*func)(PyObject *));  /* return item as double */


/* the math_sum function in mathmodule.c ... */

static PyObject*
math_sum(PyObject *self, PyObject *args)
{
	double s = 0.0;
	PyObject *seq, *sum = NULL;

	if (!PyArg_UnpackTuple(args, "sum", 1, 2, &seq, &sum))
		return NULL;

	if (float_sum(seq, sum, &s, PyFloat_AsDouble))
		return NULL;

	return PyFloat_FromDouble(s);
}



/* the cmath_sum function in cmathmodule.c ... */

static double
_cmath_sum_real(PyObject *item)
{
	Py_complex c = PyComplex_AsCComplex(item);
	return c.real;
}

static double
_cmath_sum_imag(PyObject *item)
{
	Py_complex c = PyComplex_AsCComplex(item);
	return c.imag;
}

static PyObject*
cmath_sum(PyObject *self, PyObject *args)
{
	Py_complex s = {0.0, 0.0};
	PyObject *seq, *sum = NULL;

	if (!PyArg_UnpackTuple(args, "sum", 1, 2, &seq, &sum))
		return NULL;

	if (float_sum(seq, sum, &s.real, _cmath_sum_real) ||
	    float_sum(seq, sum, &s.imag, _cmath_sum_imag))
		return NULL;

	return PyComplex_FromCComplex(s);
}



/* the shared functions in floatobject.c??? */

static double
_float_add(double x, double y)
{
	double s;

	errno = 0;
	s = x + y;
	/* following the IEEE 754 rules for summation:

	1. if the summands include a NaN, return a NaN

	2. else if the summands include infinities of
	   both signs, raise ValueError,

	3. else if the summands include infinities of
	   only one sign, return infinity with that sign,

	4. else (all summands are finite) if the result
	   is infinite, raise OverflowError.  The result
	   can never be a NaN if all summands are finite.
	*/
	if (Py_IS_NAN(s)) {
		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
			errno = EDOM;
		else
			errno = 0;  /* rule 1 */
	}
	else if (Py_IS_INFINITY(s)) {
		if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
			errno = ERANGE;  /* rule 4 */
		else if (Py_IS_INFINITY(x) && Py_IS_INFINITY(y) && x != y)
			errno = EDOM;  /* rule 2 */
		else
			errno = 0;  /* rule 3 */
	}
	return s;
}


int
float_sum(PyObject *seq, PyObject *sum, double *sum_ptr,
          double (*func)(PyObject *))  /* return item as double */
{
#	define FLOAT_SUM_PS  128  /* stack partials */

	int e = 1;  /* meaning error */
	Py_ssize_t i, j, n = 0, m = FLOAT_SUM_PS;
	double x, y, s, hi, lo, ps[FLOAT_SUM_PS], *p = ps;
	PyObject *iter, *item;

	iter = PyObject_GetIter(seq);
	if (iter == NULL)
		return e;

	if (sum != NULL) { 
		s = func(sum);  /* e.g. PyFloat_AsDouble() */
		if (PyErr_Occurred())
			goto float_sum_exit;
		p[n++] = s;
	} else
		p[0] = 0.0;

	PyFPE_START_PROTECT("sum", goto float_sum_exit)

	/* for x in iterable */
	for(;;) {
		/* some invariants */
		assert(e != 0);
		assert(n >= 0);
		assert(m >= n);
		assert((m == MATH_SUM_PS && p == ps) ||
		       (m >  MATH_SUM_PS && p != NULL));

		item = PyIter_Next(iter);
		if (item == NULL) {
			if (PyErr_Occurred())
				goto float_sum_error;
			else
				break;
		}
		x = func(item);  /* e.g. PyFloat_AsDouble() */
		Py_DECREF(item);
		if (PyErr_Occurred())
			goto float_sum_error;
		/* i = 0; for y in partials */
		for (i = j = 0; j < n; j++) {
			y = p[j];
			hi = _float_add(x, y);
			if (errno && is_error(hi))
				goto float_sum_error;
			/* lo = min_x_y - (hi - max_x_y) */
			lo = fabs(x) < fabs(y)
			   ? x - (hi - y)
			   : y - (hi - x);
			if (lo != 0.0)
				p[i++] = lo;
			x = hi;
		}
		/* ps[i:] = [x] */
		n = i;
		if (n >= m) {  /* extend partials */
			void *v = NULL;
			m += m / 2;  /* plus 50% */
			if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
				if (p == ps) {
					v = PyMem_Malloc(sizeof(double) * m);
					if (v != NULL)
						memcpy(v, ps, sizeof(double) * n);
				} else
					v = PyMem_Realloc(p, sizeof(double) * m);
			}
			if (v == NULL) {  /* no memory or size overflow */
				PyErr_SetString(PyExc_MemoryError, "math sum error");
				goto float_sum_error;
			}
			p = (double*) v;
		}
		p[n++] = x;
	}

	/* sum(ps, 0.0) */
	s = p[0];  /* always valid */
	for (i = 1; i < n; i++) {
		s = _float_add(p[i], s);
		if (errno && is_error(s))
			goto float_sum_error;
	}
	*sum_ptr = s;
	e = 0;

float_sum_error:
	PyFPE_END_PROTECT(s)

float_sum_exit:
	Py_DECREF(iter);
	if (p != ps)
		PyMem_Free(p);
	return e;

#	undef FLOAT_SUM_PS
}