mrkn	2016-12-06 22:40:31 +0900 (Tue, 06 Dec 2016)

  New Revision: 57001

  https://svn.ruby-lang.org/cgi-bin/viewvc.cgi?view=revision&revision=57001

  Log:
    array.c, enum.c: change sum algorithm
    
    * array.c (rb_ary_sum): change the algorithm to Kahan-Babuska balancing
      summation to be more precise.
      [Feature #12871] [ruby-core:77771]
    
    * enum.c (sum_iter, enum_sum): ditto.
    
    * test_array.rb, test_enum.rb: add an assertion for the above change.

  Modified files:
    trunk/array.c
    trunk/enum.c
    trunk/test/ruby/test_array.rb
    trunk/test/ruby/test_enum.rb
Index: enum.c
===================================================================
--- enum.c	(revision 57000)
+++ enum.c	(revision 57001)
@@ -3644,8 +3644,11 @@ sum_iter(VALUE i, struct enum_sum_memo * https://github.com/ruby/ruby/blob/trunk/enum.c#L3644
         }
     }
     else if (RB_FLOAT_TYPE_P(v)) {
-        /* Kahan's compensated summation algorithm */
-        double x, y, t;
+        /*
+         * Kahan-Babuska balancing compensated summation algorithm
+         * See http://link.springer.com/article/10.1007/s0????-???-????-x
+         */
+        double x, t;
 
       float_value:
         if (RB_FLOAT_TYPE_P(i))
@@ -3662,9 +3665,11 @@ sum_iter(VALUE i, struct enum_sum_memo * https://github.com/ruby/ruby/blob/trunk/enum.c#L3665
             goto some_value;
         }
 
-        y = x - c;
-        t = f + y;
-        c = (t - f) - y;
+        t = f + x;
+        if (fabs(f) >= fabs(x))
+            c += ((f - t) + x);
+        else
+            c += ((x - t) + f);
         f = t;
     }
     else {
@@ -3788,7 +3793,7 @@ enum_sum(int argc, VALUE* argv, VALUE ob https://github.com/ruby/ruby/blob/trunk/enum.c#L3793
         rb_block_call(obj, id_each, 0, 0, enum_sum_i, (VALUE)&memo);
 
     if (memo.float_value) {
-        return DBL2NUM(memo.f);
+        return DBL2NUM(memo.f + memo.c);
     }
     else {
         if (memo.n != 0)
Index: test/ruby/test_array.rb
===================================================================
--- test/ruby/test_array.rb	(revision 57000)
+++ test/ruby/test_array.rb	(revision 57001)
@@ -2821,6 +2821,7 @@ class TestArray < Test::Unit::TestCase https://github.com/ruby/ruby/blob/trunk/test/ruby/test_array.rb#L2821
     assert_float_equal(large_number+(small_number*10), [large_number, *[small_number]*10].sum)
     assert_float_equal(large_number+(small_number*10), [large_number/1r, *[small_number]*10].sum)
     assert_float_equal(large_number+(small_number*11), [small_number, large_number/1r, *[small_number]*10].sum)
+    assert_float_equal(small_number, [large_number, small_number, -large_number].sum)
 
     assert_equal("abc", ["a", "b", "c"].sum(""))
     assert_equal([1, [2], 3], [[1], [[2]], [3]].sum([]))
Index: test/ruby/test_enum.rb
===================================================================
--- test/ruby/test_enum.rb	(revision 57000)
+++ test/ruby/test_enum.rb	(revision 57001)
@@ -903,6 +903,7 @@ class TestEnumerable < Test::Unit::TestC https://github.com/ruby/ruby/blob/trunk/test/ruby/test_enum.rb#L903
     assert_float_equal(large_number+(small_number*10), [large_number, *[small_number]*10].each.sum)
     assert_float_equal(large_number+(small_number*10), [large_number/1r, *[small_number]*10].each.sum)
     assert_float_equal(large_number+(small_number*11), [small_number, large_number/1r, *[small_number]*10].each.sum)
+    assert_float_equal(small_number, [large_number, small_number, -large_number].each.sum)
 
     assert_equal("abc", ["a", "b", "c"].each.sum(""))
     assert_equal([1, [2], 3], [[1], [[2]], [3]].each.sum([]))
Index: array.c
===================================================================
--- array.c	(revision 57000)
+++ array.c	(revision 57001)
@@ -5796,14 +5796,17 @@ rb_ary_sum(int argc, VALUE *argv, VALUE https://github.com/ruby/ruby/blob/trunk/array.c#L5796
     }
 
     if (RB_FLOAT_TYPE_P(e)) {
-        /* Kahan's compensated summation algorithm */
+        /*
+         * Kahan-Babuska balancing compensated summation algorithm
+         * See http://link.springer.com/article/10.1007/s0????-???-????-x
+         */
         double f, c;
 
         f = NUM2DBL(v);
         c = 0.0;
         goto has_float_value;
         for (; i < RARRAY_LEN(ary); i++) {
-            double x, y, t;
+            double x, t;
             e = RARRAY_AREF(ary, i);
             if (block_given)
                 e = rb_yield(e);
@@ -5819,11 +5822,14 @@ rb_ary_sum(int argc, VALUE *argv, VALUE https://github.com/ruby/ruby/blob/trunk/array.c#L5822
             else
                 goto not_float;
 
-            y = x - c;
-            t = f + y;
-            c = (t - f) - y;
+            t = f + x;
+            if (fabs(f) >= fabs(x))
+                c += ((f - t) + x);
+            else
+                c += ((x - t) + f);
             f = t;
         }
+        f += c;
         return DBL2NUM(f);
 
       not_float:

--
ML: ruby-changes@q...
Info: http://www.atdot.net/~ko1/quickml/