From 2bfdddd4dc12d85737577f175fbc7a1a7addacb7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 10:53:21 -0400 Subject: [PATCH 1/8] Leverage Base reduction functions for sum Co-authored-by: Simon Byrne --- src/KahanSummation.jl | 84 ++++++++++++++++++++++++++++++++----------- 1 file changed, 64 insertions(+), 20 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 01f51a3..c894ae9 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -65,30 +65,74 @@ 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 - s - c + 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 + 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 From 050d49d84389a630d3d0883942bdd44481980515 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 10:55:02 -0400 Subject: [PATCH 2/8] bump version + make Julia compat just "1" --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 789c805..399d99b 100644 --- a/Project.toml +++ b/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" From 0ec75a997b6e3a8fbf39efaee9f82e13526e09f4 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 10:55:33 -0400 Subject: [PATCH 3/8] Fix typo --- src/KahanSummation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index c894ae9..1a31312 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -130,7 +130,7 @@ Base.mapreduce_first(f, ::typeof(plus_kbn), x) = TwicePrecisionN(x, zero(x)) 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(f, X; kw...) = singleprec(mapreduce(f, plus_kbn, X; kw...)) sum_kbn(X; kw...) = sum_kbn(identity, X; kw...) From a6e4bef330cf6958fe1ae5ccac8e06b3d0cd3e26 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 11:02:04 -0400 Subject: [PATCH 4/8] add comments explaining why cumsum and sum are implemented differently --- src/KahanSummation.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 1a31312..ecd9f26 100644 --- a/src/KahanSummation.jl +++ b/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) From ee2d697bedb126ca2ce5324985d1f12bb4b6f00f Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 11:11:00 -0400 Subject: [PATCH 5/8] attempt to bring coverage back --- test/runtests.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 45aaedf..86d7c2e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,3 +49,15 @@ end @test isequal(sum_kbn(1:3), 6) @test isequal(sum_kbn((i for i in [1,2,3])), 6) end + +@testset "twice-precision addition" + # 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 +end \ No newline at end of file From f50b57c90f9d45e11d7a5bc026808049dd7a937b Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 11:14:16 -0400 Subject: [PATCH 6/8] more test coverage --- test/runtests.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 86d7c2e..a53734c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,6 +48,9 @@ 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" @@ -60,4 +63,8 @@ end i12 = KahanSummation.plus_kbn(i1, i2) f12 = KahanSummation.singleprec(i12) @test f12 == 1 + + i1 = convert(KahanSummation.TwicePrecisionN{Float64}, 5) + i2 = convert(KahanSummation.TwicePrecisionN{Float64}, 3) + @test KahanSummation.singleprec(KahanSummation.plus_kbn(i1, i2)) == 5.0 + 3.0 end \ No newline at end of file From b735870688cc533157c29a30459a73e22366b397 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 11:15:36 -0400 Subject: [PATCH 7/8] fix typo --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index a53734c..100ff07 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,7 +53,7 @@ end @test sum_kbn(x -> x + 1, Float64[]) == zero(Float64) end -@testset "twice-precision addition" +@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 @@ -64,7 +64,9 @@ end 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 \ No newline at end of file From c659cf8ce8f1b74618bc42653e3a4d8259750de7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 10 May 2025 11:18:38 -0400 Subject: [PATCH 8/8] forgot that it's stored in negated form --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 100ff07..cfbe1bb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,11 +62,11 @@ end i2 = KahanSummation.TwicePrecisionN{Float64}(-1e100, 1) i12 = KahanSummation.plus_kbn(i1, i2) f12 = KahanSummation.singleprec(i12) - @test f12 == 1 + @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 \ No newline at end of file +end