get pairwise-sum accuracy for cumsum (fix #9648)

Author: stevengj

base/array.jl

@@ -1439,20 +1439,24 @@ symdiff(a, b, rest...) = symdiff(a, symdiff(b, rest...))
 _cumsum_type{T<:Number}(v::AbstractArray{T}) = typeof(+zero(T))
 _cumsum_type(v) = typeof(v[1]+v[1])
 
-for (f, fp, op) = ((:cumsum, :cumsum_pairwise, :+),
-                   (:cumprod, :cumprod_pairwise, :*) )
+for (f, fp, op) = ((:cumsum, :cumsum_pairwise!, :+),
+                   (:cumprod, :cumprod_pairwise!, :*) )
     # in-place cumsum of c = s+v(i1:n), using pairwise summation as for sum
-    @eval function ($fp)(v::AbstractVector, c::AbstractVector, s, i1, n)
+    @eval function ($fp){T}(v::AbstractVector, c::AbstractVector{T}, s, i1, n)
+        local s_::T # for sum(v(i1:n)), i.e. sum without s
         if n < 128
-            @inbounds c[i1] = ($op)(s, v[i1])
+            @inbounds s_ = v[i1]
+            @inbounds c[i1] = ($op)(s, s_)
             for i = i1+1:i1+n-1
-                @inbounds c[i] = $(op)(c[i-1], v[i])
+                @inbounds s_ = $(op)(s_, v[i])
+                @inbounds c[i] = $(op)(s, s_)
             end
         else
             n2 = div(n,2)
-            ($fp)(v, c, s, i1, n2)
-            ($fp)(v, c, c[(i1+n2)-1], i1+n2, n-n2)
+            s_ = ($fp)(v, c, s, i1, n2)
+            s_ = $(op)(s_, ($fp)(v, c, s + s_, i1+n2, n-n2))
         end
+        return s_
     end
 
     @eval function ($f)(v::AbstractVector)

test/arrayops.jl

@@ -974,3 +974,8 @@ a = [ [ 1 0 0 ], [ 0 0 0 ] ]
 @test rotl90(a,4) == a
 @test rotr90(a,4) == a
 @test rot180(a,2) == a
+
+# issue #9648
+let x = fill(1.5f0, 10^7)
+    @test abs(1.5f7 - cumsum(x)[end]) < 3*eps(1.5f7)
+end