Leverage Base reduction functions for sum

Project.toml

@@ -1,9 +1,9 @@
 name = "KahanSummation"
 uuid = "8e2b3108-d4c1-50be-a7a2-16352aec75c3"
-version = "0.3.1"
+version = "0.4.0"
 
 [compat]
-julia = "1.0"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/KahanSummation.jl

@@ -13,6 +13,12 @@ summation algorithm for additional accuracy.
 """
 function cumsum_kbn end
 
+# Note that the implementation for cumsum_kbn will stay as is,
+# since the use of a temporary variable for accumulation in `accumulate!`
+# is an implementation detail and not guaranteed!
+
+# sum is implemented using reducing functions but cumsum cannot be.
+
 cumsum_kbn(x::AbstractArray; dims=:) = _cumsum_kbn(x, dims)
 cumsum_kbn(x; dims=:) = _cumsum_kbn(collect(x), dims)
 
@@ -65,31 +71,75 @@ function _cumsum_kbn(v::AbstractArray{T}, ::Colon) where {T}
 end
 
 """
-    sum_kbn(A)
+    TwicePrecisionN{T}(number)
+    TwicePrecisionN{T}(hi, nlo)
 
-Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated
-summation algorithm for additional accuracy.
+Represents an extended precision number as `x.hi - x.nlo`.
+We store the lower order component as the negation to avoid problems when `x.hi == -0.0`.
+
+This does not subtype Number or Real, being meant primarily for internal use
+in KahanSummation.jl.
+
+Convert a `TwicePrecisionN{T}` back to a `T` by calling `singleprec(tp)`.
 """
-function sum_kbn(A)
-    T = Base.@default_eltype(A)
-    c = Base.reduce_empty(+, T)
-    it = iterate(A)
-    it === nothing && return c
-    Ai, i = it
-    s = Ai - c
-    while (it = iterate(A, i)) !== nothing
-        Ai, i = it::Tuple{T, Int}
-        t = s + Ai
-        if abs(s) >= abs(Ai)
-            c -= ((s-t) + Ai)
-        else
-            c -= ((Ai-t) + s)
-        end
-        s = t
+struct TwicePrecisionN{T}
+    hi::T
+    nlo::T
+end
+
+singleprec(x::TwicePrecisionN{T}) where {T} = convert(T, x)
+
+# Implement Base methods
+Base.convert(::Type{TwicePrecisionN{T}}, x::Number) where {T} =
+    TwicePrecisionN{T}(convert(T, x), zero(T))
+Base.convert(::Type{T}, x::TwicePrecisionN) where {T} =
+    convert(T, x.hi - x.nlo)
+
+# Two-sum implementation
+@inline function plus_kbn(x::T, y::T) where {T}
+    hi = x + y
+    nlo = abs(x) > abs(y) ? (hi - x ) - y : (hi - y) - x
+    TwicePrecisionN(hi, nlo)
+end
+@inline function plus_kbn(x::T, y::TwicePrecisionN{T}) where {T}
+    hi = x + y.hi
+    if abs(x) > abs(y.hi)
+        nlo = ((hi - x) - y.hi) + y.nlo
+    else
+        nlo = ((hi - y.hi) - x) + y.nlo
+    end
+    TwicePrecisionN(hi, nlo)
+end
+@inline plus_kbn(x::TwicePrecisionN{T}, y::T) where {T} = plus_kbn(y, x)
+
+@inline function plus_kbn(x::TwicePrecisionN{T}, y::TwicePrecisionN{T}) where {T}
+    hi = x.hi + y.hi
+    if abs(x.hi) > abs(y.hi)
+        nlo = (((hi - x.hi) - y.hi) + y.nlo) + x.nlo
+    else
+        nlo = (((hi - y.hi) - x.hi) + x.nlo) + y.nlo
     end
-    s - c
+    TwicePrecisionN(hi, nlo)
 end
 
+# Implement methods for accumulators, specifically mapreduce
+Base.mapreduce_empty(f, ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T))
+Base.mapreduce_empty(::typeof(identity), ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T)) # disambiguate
+Base.mapreduce_first(f, ::typeof(plus_kbn), x) = TwicePrecisionN(x, zero(x))
+
+# Finally, the implementation of `sum_kbn` is trivial, dispatching to `mapreduce`.  
+# Most of the work happens in `plus_kbn`.
+
+"""
+    sum_kbn([f,] A)
+
+Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated
+summation algorithm for additional accuracy.
+"""
+sum_kbn(f, X; kw...) = singleprec(mapreduce(f, plus_kbn, X; kw...))
+sum_kbn(X; kw...) = sum_kbn(identity, X; kw...)
+
+
 ### Deprecations
 
 Base.@deprecate cumsum_kbn(A, axis) cumsum_kbn(A; dims=axis)

test/runtests.jl

@@ -48,4 +48,25 @@ end
     @test sum_kbn(Iterators.filter(isodd, 1:10)) == 25
     @test isequal(sum_kbn(1:3), 6)
     @test isequal(sum_kbn((i for i in [1,2,3])), 6)
+    # also test `sum_kbn(f, x)`
+    @test sum_kbn(x -> x + 1, [7 8 9]) == sum_kbn([7 8 9] .+ 1)
+    @test sum_kbn(x -> x + 1, Float64[]) == zero(Float64)
+end
+
+@testset "twice-precision addition" begin
+    # Note: the functions and types used here are internal
+
+    # The intent of this test is to make sure that the two-sum
+    # works on two twiceprecision numbers correctly, nothing else.
+    i1 = convert(KahanSummation.TwicePrecisionN{Float64}, 1e100)
+    i2 = KahanSummation.TwicePrecisionN{Float64}(-1e100, 1)
+    i12 = KahanSummation.plus_kbn(i1, i2)
+    f12 = KahanSummation.singleprec(i12)
+    @test f12 == -1
+
+    # test the bottom if statement in sum_kbn too, or try to at least
+    i1 = convert(KahanSummation.TwicePrecisionN{Float64}, 5)
+    i2 = convert(KahanSummation.TwicePrecisionN{Float64}, 3)
+    @test KahanSummation.singleprec(KahanSummation.plus_kbn(i1, i2)) == 5.0 + 3.0
+    @test KahanSummation.singleprec(KahanSummation.plus_kbn(i2, i1)) == 5.0 + 3.0
 end