diff --git a/base/abstractarraymath.jl b/base/abstractarraymath.jl index 352d438525a26..d7496360e2f5c 100644 --- a/base/abstractarraymath.jl +++ b/base/abstractarraymath.jl @@ -238,73 +238,6 @@ function circshift(a::AbstractArray, shiftamt) circshift!(similar(a), a, map(Integer, (shiftamt...,))) end -# Uses K-B-N summation -function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat - r = similar(v) - if isempty(v); return r; end - - inds = indices(v, 1) - i1 = first(inds) - s = r[i1] = v[i1] - c = zero(T) - for i=i1+1:last(inds) - 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 -# TODO: Needs a separate IndexCartesian method, this is only fast for IndexLinear - -""" - cumsum_kbn(A, dim::Integer) - -Cumulative sum along a dimension, using the Kahan-Babuska-Neumaier compensated summation -algorithm for additional accuracy. -""" -function cumsum_kbn(A::AbstractArray{T}, axis::Integer) where T<:AbstractFloat - 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 - ## Other array functions ## """ diff --git a/base/deprecated.jl b/base/deprecated.jl index 82a44c952a6e4..68960bb291104 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -2136,7 +2136,6 @@ end @deprecate RowVector{T}(n::Tuple{Int,Int}) where {T} RowVector{T}(uninitialized, n) @deprecate cumsum(A::AbstractArray) cumsum(A, 1) -@deprecate cumsum_kbn(A::AbstractArray) cumsum_kbn(A, 1) @deprecate cumprod(A::AbstractArray) cumprod(A, 1) # issue #16307 @@ -2183,6 +2182,10 @@ end @deprecate SSHCredentials SSHCredential false end +# issue #24804 +@deprecate_moved sum_kbn "KahanSummation" +@deprecate_moved cumsum_kbn "KahanSummation" + # END 0.7 deprecations # BEGIN 1.0 deprecations diff --git a/base/exports.jl b/base/exports.jl index f99dde7d85e12..ed3f5d2b34b49 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -433,7 +433,6 @@ export cumsum!, accumulate, accumulate!, - cumsum_kbn, eachindex, extrema, fill!, @@ -527,7 +526,6 @@ export sub2ind, sum!, sum, - sum_kbn, to_indices, vcat, vec, diff --git a/base/reduce.jl b/base/reduce.jl index eb9d3f804801f..e1fb0989aa669 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -391,38 +391,6 @@ julia> sum(1:20) sum(a) = mapreduce(promote_sys_size_add, +, a) sum(a::AbstractArray{Bool}) = count(a) - -# Kahan (compensated) summation: O(1) error growth, at the expense -# of a considerable increase in computational expense. - -""" - sum_kbn(A) - -Returns the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated -summation algorithm for additional accuracy. -""" -function sum_kbn(A) - T = @default_eltype(typeof(A)) - c = promote_sys_size_add(zero(T)::T) - i = start(A) - if done(A, i) - return c - end - Ai, i = next(A, i) - s = Ai - c - while !(done(A, i)) - Ai, i = next(A, i) - t = s + Ai - if abs(s) >= abs(Ai) - c -= ((s-t) + Ai) - else - c -= ((Ai-t) + s) - end - s = t - end - s - c -end - ## prod """ prod(f, itr) diff --git a/doc/src/manual/style-guide.md b/doc/src/manual/style-guide.md index e586df28e48cd..6135d734acf21 100644 --- a/doc/src/manual/style-guide.md +++ b/doc/src/manual/style-guide.md @@ -159,7 +159,7 @@ uses (e.g. `a[i]::Int`) than to try to pack many alternatives into one type. * functions are lowercase ([`maximum`](@ref), [`convert`](@ref)) and, when readable, with multiple words squashed together ([`isequal`](@ref), [`haskey`](@ref)). When necessary, use underscores as word separators. Underscores are also used to indicate a combination of concepts ([`remotecall_fetch`](@ref) - as a more efficient implementation of `fetch(remotecall(...))`) or as modifiers ([`sum_kbn`](@ref)). + as a more efficient implementation of `fetch(remotecall(...))`) or as modifiers. * conciseness is valued, but avoid abbreviation ([`indexin`](@ref) rather than `indxin`) as it becomes difficult to remember whether and how particular words are abbreviated. diff --git a/doc/src/stdlib/arrays.md b/doc/src/stdlib/arrays.md index 9ec8c580cc1fa..c39546121a7c9 100644 --- a/doc/src/stdlib/arrays.md +++ b/doc/src/stdlib/arrays.md @@ -147,7 +147,6 @@ Base.cumprod Base.cumprod! Base.cumsum Base.cumsum! -Base.cumsum_kbn Base.LinAlg.diff Base.repeat(::AbstractArray) Base.rot180 @@ -156,7 +155,6 @@ Base.rotr90 Base.reducedim Base.mapreducedim Base.mapslices -Base.sum_kbn ``` ## Combinatorics diff --git a/test/arrayops.jl b/test/arrayops.jl index 7d6e178c73471..f3a4a0015d1f7 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -870,29 +870,6 @@ end @test 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 - - @test isequal(cumsum_kbn(v), cv) - @test isequal(cumsum_kbn(v2), cv2) - - A = [v reverse(v) v2 reverse(v2)] - - c = cumsum_kbn(A, 1) - - @test isequal(c[:,1], cv) - @test isequal(c[:,3], cv2) - @test isequal(c[4,:], [2.0, 2.0, 2.0, 2.0]*1000) - - c = cumsum_kbn(A, 2) - - @test isequal(c[1,:], cv2) - @test isequal(c[3,:], cv) - @test isequal(c[:,4], [2.0,2.0,2.0,2.0]*1000) - @test repeat(BitMatrix(Matrix(I, 2, 2)), inner = (2,1), outer = (1,2)) == repeat(Matrix(I, 2, 2), inner = (2,1), outer = (1,2)) end diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 67c91db8f3c84..fc8682e93f7c0 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -383,14 +383,20 @@ I,J,N = findnz(z) @test vecnorm(A) ≈ vecnorm(parent(A)) @test vecdot(v, v) ≈ vecdot(v0, v0) -v = OffsetArray([1,1e100,1,-1e100], (-3,))*1000 -v2 = OffsetArray([1,-1e100,1,1e100], (5,))*1000 -@test isa(v, OffsetArray) -cv = OffsetArray([1,1e100,1e100,2], (-3,))*1000 -cv2 = OffsetArray([1,-1e100,-1e100,2], (5,))*1000 -@test isequal(cumsum_kbn(v), cv) -@test isequal(cumsum_kbn(v2), cv2) -@test isequal(sum_kbn(v), sum_kbn(parent(v))) +# Prior to its removal from Base, cumsum_kbn was used here. To achieve the same level of +# accuracy in the tests, we need to use BigFloats with enlarged precision. +@testset "high-precision array reduction" begin + setprecision(BigFloat, 500) do + v = OffsetArray(BigFloat[1,1e100,1,-1e100], (-3,)) .* 1000 + v2 = OffsetArray(BigFloat[1,-1e100,1,1e100], ( 5,)) .* 1000 + @test isa(v, OffsetArray) + cv = OffsetArray(BigFloat[1, 1e100, 1e100,2], (-3,)) .* 1000 + cv2 = OffsetArray(BigFloat[1,-1e100,-1e100,2], ( 5,)) .* 1000 + @test cumsum(v) ≈ cv + @test cumsum(v2) ≈ cv2 + @test sum(v) ≈ sum(parent(v)) + end +end io = IOBuffer() writedlm(io, A) diff --git a/test/reduce.jl b/test/reduce.jl index 265c47ce97c37..50bd4ac693e43 100644 --- a/test/reduce.jl +++ b/test/reduce.jl @@ -135,18 +135,6 @@ for f in (sum3, sum4, sum7, sum8) end @test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7])) -@test sum_kbn([1,1e100,1,-1e100]) === 2.0 -@test sum_kbn(Float64[]) === 0.0 -@test sum_kbn(i for i=1.0:1.0:10.0) === 55.0 -@test sum_kbn(i for i=1:1:10) === 55 -@test sum_kbn([1 2 3]) === 6 -@test sum_kbn([2+im 3-im]) === 5+0im -@test sum_kbn([1+im 2+3im]) === 3+4im -@test sum_kbn([7 8 9]) === sum_kbn([8 9 7]) -@test sum_kbn(i for i=1:1:10) === sum_kbn(i for i=10:-1:1) -@test sum_kbn([-0.0]) === -0.0 -@test sum_kbn([-0.0,-0.0]) === -0.0 - # check sum(abs, ...) for support of empty collections @testset "sum(abs, [])" begin @test @inferred(sum(abs, Float64[])) === 0.0 @@ -362,14 +350,13 @@ end ## cumsum, cummin, cummax z = rand(10^6) -let es = sum_kbn(z), es2 = sum_kbn(z[1:10^5]) +let es = sum(BigFloat.(z)), es2 = sum(BigFloat.(z[1:10^5])) @test (es - sum(z)) < es * 1e-13 cs = cumsum(z) @test (es - cs[end]) < es * 1e-13 @test (es2 - cs[10^5]) < es2 * 1e-13 end - @test sum(collect(map(UInt8,0:255))) == 32640 @test sum(collect(map(UInt8,254:255))) == 509