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

diff --git a/extra/math/extras/extras-tests.factor b/extra/math/extras/extras-tests.factor
index 8daa9726dfc24674af98096b726ded91047eb8d1..1d10a680d9a2b5e1c66332ae71d1866c350d9f20 100644 (file)
--- a/extra/math/extras/extras-tests.factor
+++ b/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
 

diff --git a/extra/math/extras/extras.factor b/extra/math/extras/extras.factor
index 7e99b35364d775caefc0b7e1fe779185c8d6e584..8136bdbd555797bf61004dd483c3e940a10d6ef3 100644 (file)
--- a/extra/math/extras/extras.factor
+++ b/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! ;