math.extras: adding a more exact sum for floats.

extra/math/extras/extras-tests.factor

@@ -130,6 +130,9 @@ tools.test ;
 { { 0 1 2 3 0 0 1 } } [ { 1 2 3 3 2 1 2 } [ <= ] monotonic-count ] unit-test
 { 4 } [ { 1 2 3 1 2 3 4 5 } [ < ] max-monotonic-count ] unit-test
 
+{ 4.0 } [ { 1e-30 1 3 -1e-30 } sum-floats ] unit-test
+{ 1.0000000000000002e16 } [ { 1e-16 1 1e16 } sum-floats ] unit-test
+
 { 2470 } [ 20 <iota> sum-squares ] unit-test
 { 2470 } [ 20 <iota> >array sum-squares ] unit-test
 

extra/math/extras/extras.factor

@@ -290,6 +290,53 @@ PRIVATE>
     [ 0.0 0.0 ] 2dip [ 2dip rot kahan+ ] curry
     [ -rot ] prepose each nip ; inline
 
+<PRIVATE
+
+! Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates
+! www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
+
+: sort-partial ( x y -- x' y' )
+    2dup [ abs ] bi@ < [ swap ] when ; inline
+
+:: partial+ ( x y -- hi lo )
+    x y + dup x - :> yr y yr - ; inline
+
+:: partial-sums ( seq -- seq' )
+    V{ } clone :> partials
+    seq [
+        0 partials [
+            swapd sort-partial partial+ swapd
+            [ over partials set-nth 1 + ] unless-zero
+        ] each :> i
+        i partials shorten
+        [ i partials set-nth ] unless-zero
+    ] each partials ;
+
+:: sum-exact ( partials -- n )
+    partials empty? [ 0.0 ] [
+        ! sum from the top, stop when sum becomes inexact
+        0.0 0.0 partials [
+            nip partial+ dup 0.0 = not
+        ] find-last drop :> ( lo n )
+
+        ! make half-even rounding work across multiple partials
+        n [ 0 > ] [ f ] if* [
+            n 1 - partials nth
+            [ 0.0 < lo 0.0 < and ]
+            [ 0.0 > lo 0.0 > and ] bi or [
+                lo 2.0 * :> y
+                dup y + :> x
+                x over - :> yr
+                y yr = [ drop x ] when
+            ] when
+        ] when
+    ] if ;
+
+PRIVATE>
+
+: sum-floats ( seq -- n )
+    partial-sums sum-exact ;
+
 ! SYNTAX: .. dup pop scan-object [a,b) suffix! ;
 ! SYNTAX: ... dup pop scan-object [a,b] suffix! ;