inspired by this thread and my comment about the implementation of a sum function from the math library of Python, I converted this function to Purebasic.
Here are a few more links:
- Source of original function: https://github.com/python/cpython/blob/ ... le.c#L1356
- More detailed explanation and alternatives: https://code.activestate.com/recipes/393090/
In the end there is a small example summing up 40000 specific values which sum up to the correct value 20000.0 when using fsum but sum up to 0.0 with the naive implementation.
Feel free to use it where you need it:
Code: Select all
Macro assert(expression)
CompilerIf #PB_Compiler_Debugger
If Not (expression)
DebuggerError("assert failed at " + Str(#PB_Compiler_Line - 2))
EndIf
CompilerEndIf
EndMacro
Macro IsFinite(float)
(Not (IsNAN(float) Or IsInfinity(float)))
EndMacro
;#define NUM_PARTIALS 32 /* initial partials array size, on stack */
#NUM_PARTIALS = 32
; Source: https://github.com/python/cpython/blob/main/Modules/mathmodule.c#L1356
;Static PyObject *
;math_fsum(PyObject *Module, PyObject *seq)
;/*[clinic End generated code: output=ba5c672b87fe34fc input=c51b7d8caf6f6e82]*/
;{
Procedure.d fsum(Array seq.d(1))
;PyObject *item, *iter, *sum = NULL;
Protected item.i = 0, seq_size.i, sum.d
;Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
Protected i.i, j.i, n.i = 0, m.i = #NUM_PARTIALS
;double x, y, t, ps[NUM_PARTIALS], *p = ps;
Protected Dim ps.d(#NUM_PARTIALS - 1)
Protected x.d, y.d
; t is not needed because of Swap command
; *p = ps is not needed because we can access ps() directly
;double xsave, special_sum = 0.0, inf_sum = 0.0;
Protected.d xsave, special_sum = 0.0, inf_sum = 0.0
;double hi, yr, lo = 0.0;
Protected.d hi, yr, lo = 0.0
;
;iter = PyObject_GetIter(seq);
seq_size = ArraySize(seq()) + 1
;if (iter == NULL)
If seq_size = 0
;return NULL;
ProcedureReturn 0
EndIf
;
;for(;;) { /* for x in iterable */
While item < seq_size
;assert(0 <= n && n <= m);
assert(0 <= n And n <= m)
;assert((m == NUM_PARTIALS && p == ps) ||
;(m > NUM_PARTIALS && p != NULL));
;
;item = PyIter_Next(iter);
;if (item == NULL) {
; If (PyErr_Occurred())
; goto _fsum_error;
; break;
;}
;ASSIGN_DOUBLE(x, item, error_with_item);
x = seq(item)
;Py_DECREF(item);
;
;xsave = x;
xsave = x
;For (i = j = 0; j < n; j++) { /* for y in partials */
i = 0
j = 0
While j < n
;y = p[j];
y = ps(j)
;if (fabs(x) < fabs(y)) {
If Abs(x) < Abs(y)
;t = x; x = y; y = t;
Swap x, y
;}
EndIf
;hi = x + y;
hi = x + y
;yr = hi - x;
yr = hi - x
;lo = y - yr;
lo = y - yr
;If (lo != 0.0)
If lo <> 0.0
;p[i++] = lo;
ps(i) = lo
i + 1
EndIf
;x = hi;
x = hi
;}
j + 1
Wend
;
;n = i; /* ps[i:] = [x] */
n = i
;if (x != 0.0) {
If x <> 0.0
;if (! Py_IS_FINITE(x)) {
If Not IsFinite(x)
;/* a nonfinite x could arise either as a result of intermediate overflow, or as a result of a nan or inf in the summands */
;if (Py_IS_FINITE(xsave)) {
If IsFinite(xsave)
;PyErr_SetString(PyExc_OverflowError, "intermediate overflow in fsum");
DebuggerError("intermediate overflow in fsum")
;goto _fsum_error;
Goto _fsum_error
;}
EndIf
;if (Py_IS_INFINITY(xsave))
If IsInfinity(xsave)
; inf_sum += xsave;
inf_sum + xsave
EndIf
;special_sum += xsave;
special_sum + xsave
;/* reset partials */
;n = 0;
n = 0
;}
;else if (n >= m && _fsum_realloc(&p, n, ps, &m))
ElseIf (n >= m) ; here we can just use ReDim and don't care about too less memory
;goto _fsum_error;
m + m
ReDim ps(m - 1)
;else
Else
;p[n++] = x;
ps(n) = x
n + 1
EndIf
;}
EndIf
;}
item + 1
Wend
;
;If (special_sum != 0.0) {
If special_sum <> 0.0
;if (Py_IS_NAN(inf_sum))
If IsNAN(inf_sum)
;PyErr_SetString(PyExc_ValueError, "-inf + inf in fsum");
DebuggerError("-inf + inf in fsum")
;else
Else
;sum = PyFloat_FromDouble(special_sum);
sum = special_sum
EndIf
;goto _fsum_error;
Goto _fsum_error
;}
EndIf
;
;hi = 0.0;
hi = 0.0
;If (n > 0) {
If n > 0
;hi = p[--n];
n - 1
hi = ps(n)
;/* sum_exact(ps, hi) from the top, stop when the sum becomes inexact. */
;While (n > 0) {
While n > 0
;x = hi;
x = hi
;y = p[--n];
n - 1
y = ps(n)
;assert(fabs(y) < fabs(x));
assert(Abs(y) < Abs(x))
;hi = x + y;
hi = x + y
;yr = hi - x;
yr = hi - x
;lo = y - yr;
lo = y - yr
;If (lo != 0.0)
If lo <> 0.0
;Break;
Break
EndIf
;}
Wend
;/* Make half-even rounding work across multiple partials.
; Needed so that sum([1e-16, 1, 1e16]) will round-up the last
; digit to two instead of down To zero (the 1e-16 makes the 1
; slightly closer to two). With a potential 1 ULP rounding
; error fixed-up, math.fsum() can guarantee commutativity. */
;If (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || (lo > 0.0 && p[n-1] > 0.0))) {
If n > 0 And ((lo < 0.0 And ps(n - 1) < 0.0) Or (lo > 0.0 And ps(n - 1) > 0.0))
;y = lo * 2.0;
y = lo * 2.0
;x = hi + y;
x = hi + y
;yr = x - hi;
yr = x - hi
;If (y == yr)
If y = yr
;hi = x;
hi = x
EndIf
;}
EndIf
;}
EndIf
;sum = PyFloat_FromDouble(hi);
sum = hi
;
;_fsum_error:
_fsum_error:
;Py_DECREF(iter);
;If (p != ps)
;PyMem_Free(p);
;Return sum;
ProcedureReturn sum
;
; error_with_item:
; Py_DECREF(item);
; Goto _fsum_error;
; }
EndProcedure
; Create test case, see https://code.activestate.com/recipes/393090/
Debug "Create array of 40000 summands in the form [1, 1e100, 1, -1e100] * 10000"
Define i.i
Dim summands.d(39999)
For i = 0 To 39999 Step 4
summands(i + 0) = 1
summands(i + 1) = 1e100
summands(i + 2) = 1
summands(i + 3) = -1e100
Next
For i = 0 To 3
Debug "summands(" + i + ") = " + summands(i)
Next
Debug "..."
Debug ""
; fsum
Debug "Sum using fsum:"
Debug fsum(summands())
Debug ""
; simple sum
Define sum.d = 0
For i = 0 To 39999
sum + summands(i)
Next
Debug "Naive sum:"
Debug sum
