Updated cumsum to use K-B-N summation for float arrays.

base/array.jl

@@ -1448,6 +1448,64 @@ function sum{T<:FloatingPoint}(A::StridedArray{T})
     s + c
 end
 
+# Uses K-B-N summation
+function cumsum{T<:FloatingPoint}(v::StridedVector{T})
+    n = length(v)
+    r = similar(v, n)
+    if n == 0; return r; end
+
+    s = r[1] = v[1]
+    c = zero(T)
+    for i=2:n
+        vi = v[i]
+        t = s + vi
+        if abs(s) >= abs(vi)
+            c += ((s-t) + vi)
+        else
+            c += ((vi-t) + s)
+        end
+        s = t
+        r[i] = s+c
+    end
+    return r
+end
+
+# Uses K-B-N summation
+function cumsum{T<:FloatingPoint}(A::StridedArray{T}, axis::Integer)
+    dimsA = size(A)
+    ndimsA = ndims(A)
+    axis_size = dimsA[axis]
+    axis_stride = 1
+    for i = 1:(axis-1)
+        axis_stride *= size(A,i)
+    end
+
+    if axis_size <= 1
+        return A
+    end
+
+    B = similar(A)
+    C = similar(A)
+
+    for i = 1:length(A)
+        if div(i-1, axis_stride) % axis_size == 0
+            B[i] = A[i]
+            C[i] = zero(T)
+        else
+            s = B[i-axis_stride]
+            Ai = A[i]
+            B[i] = t = s + Ai
+            if abs(s) >= abs(Ai)
+                C[i] = C[i-axis_stride] + ((s-t) + Ai)
+            else
+                C[i] = C[i-axis_stride] + ((Ai-t) + s)
+            end
+        end
+    end
+
+    return B + C
+end
+
 function prod{T}(A::StridedArray{T})
     if isempty(A)
         return one(T)

test/arrayops.jl

@@ -217,7 +217,7 @@ d = drand(10,10)
 ## cumsum
 
 begin
-    local A, A1, A2, A3
+    local A, A1, A2, A3, v, v2, cv, cv2, c
     A = ones(Int,2,3,4)
     A1 = reshape(repmat([1,2],1,12),2,3,4)
     A2 = reshape(repmat([1 2 3],2,4),2,3,4)
@@ -241,4 +241,28 @@ begin
         @assert c[1,2] == A[1,1]+A[1,3]
         @assert c[2,2] == A[3,1]+A[3,3]
     end
+
+    v   = [1,1e100,1,-1e100]*1000
+    v2  = [1,-1e100,1,1e100]*1000
+
+    cv  = [1,1e100,1e100,2]*1000
+    cv2 = [1,-1e100,-1e100,2]*1000
+
+    @assert isequal(cumsum(v), cv)
+    @assert isequal(cumsum(v2), cv2)
+
+    A = [v reverse(v) v2 reverse(v2)]
+
+    c = cumsum(A, 1)
+
+    @assert isequal(c[:,1], cv)
+    @assert isequal(c[:,3], cv2)
+    @assert isequal(c[4,:], [2.0 2.0 2.0 2.0]*1000)
+
+    c = cumsum(A, 2)
+
+    @assert isequal(c[1,:], cv2')
+    @assert isequal(c[3,:], cv')
+    @assert isequal(c[:,4], [2.0,2.0,2.0,2.0]*1000)
+
 end