From 3fd5ce452b70c35ed65234eaefe5e8859a05d753 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 12 Dec 2017 09:47:49 -0800 Subject: [PATCH 1/4] leverage Base reduction functions --- src/KahanSummation.jl | 81 ++++++++++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 28 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index fde6700..6837dd4 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -29,6 +29,59 @@ else promote_sys_size_add(x::T) where {T} = Base.r_promote(+, zero(T)::T) end + +import Base.TwicePrecision + + +function plus_kbn(x::T, y::T) where {T} + hi = x + y + lo = abs(x) > abs(y) ? (x - hi) + y : (y - hi) + x + TwicePrecision(hi, lo) +end +function plus_kbn(x::T, y::TwicePrecision{T}) where {T} + hi = x + y.hi + if abs(x) > abs(y.hi) + lo = ((x - hi) + y.hi) + y.lo + else + lo = ((y.hi - hi) + x) + y.lo + end + TwicePrecision(hi, lo) +end +plus_kbn(x::TwicePrecision{T}, y::T) where {T} = plus_kbn(y, x) + +function plus_kbn(x::TwicePrecision{T}, y::TwicePrecision{T}) where {T} + hi = x.hi + y.hi + if abs(x.hi) > abs(y.hi) + lo = (((x.hi - hi) + y.hi) + y.lo) + x.lo + else + lo = (((y.hi - hi) + x.hi) + x.lo) + y.lo + end + TwicePrecision(hi, lo) +end + +Base.r_promote_type(::typeof(plus_kbn), ::Type{T}) where {T<:AbstractFloat} = + TwicePrecision{T} + +Base.mr_empty(f, ::typeof(plus_kbn), T) = TwicePrecision{T} + +singleprec(x::TwicePrecision{T}) where {T} = convert(T, x) + + +""" + 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) = singleprec(mapreduce(f, plus_kbn, X)) +sum_kbn(X) = sum_kbn(identity, X) + + + + + + + """ cumsum_kbn(A, dim::Integer) @@ -85,32 +138,4 @@ function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat return r end -""" - sum_kbn(A) - -Return 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 - end # module From 3e65147375556f6bbfae70dbdb183972d0a405db Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 12 Dec 2017 10:49:47 -0800 Subject: [PATCH 2/4] use own TwicePrecision type to solve -0.0 problems --- src/KahanSummation.jl | 46 ++++++++++++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 16 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 6837dd4..26c45dd 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -29,42 +29,56 @@ else promote_sys_size_add(x::T) where {T} = Base.r_promote(+, zero(T)::T) end +""" + TwicePrecisionN{T} -import Base.TwicePrecision +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`. +""" +struct TwicePrecisionN{T} + hi::T + nlo::T +end function plus_kbn(x::T, y::T) where {T} hi = x + y - lo = abs(x) > abs(y) ? (x - hi) + y : (y - hi) + x - TwicePrecision(hi, lo) + nlo = abs(x) > abs(y) ? (hi - x ) - y : (hi - y) - x + TwicePrecisionN(hi, lo) end -function plus_kbn(x::T, y::TwicePrecision{T}) where {T} +function plus_kbn(x::T, y::TwicePrecisionN{T}) where {T} hi = x + y.hi if abs(x) > abs(y.hi) - lo = ((x - hi) + y.hi) + y.lo + nlo = ((hi - x) - y.hi) + y.nlo else - lo = ((y.hi - hi) + x) + y.lo + nlo = ((hi - y.hi) - x) + y.nlo end - TwicePrecision(hi, lo) + TwicePrecisionN(hi, nlo) end -plus_kbn(x::TwicePrecision{T}, y::T) where {T} = plus_kbn(y, x) +plus_kbn(x::TwicePrecisionN{T}, y::T) where {T} = plus_kbn(y, x) -function plus_kbn(x::TwicePrecision{T}, y::TwicePrecision{T}) where {T} +function plus_kbn(x::TwicePrecisionN{T}, y::TwicePrecisionN{T}) where {T} hi = x.hi + y.hi if abs(x.hi) > abs(y.hi) - lo = (((x.hi - hi) + y.hi) + y.lo) + x.lo + nlo = (((hi - x.hi) - y.hi) + y.nlo) + x.nlo else - lo = (((y.hi - hi) + x.hi) + x.lo) + y.lo + nlo = (((hi - y.hi) - x.hi) + x.nlo) + y.nlo end - TwicePrecision(hi, lo) + TwicePrecisionN(hi, nlo) end -Base.r_promote_type(::typeof(plus_kbn), ::Type{T}) where {T<:AbstractFloat} = - TwicePrecision{T} +Base.r_promote_type(::typeof(plus_kbn), ::Type{T}) where {T} = + TwicePrecisionN{T} + +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) -Base.mr_empty(f, ::typeof(plus_kbn), T) = TwicePrecision{T} +Base.mr_empty(f, ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T)) -singleprec(x::TwicePrecision{T}) where {T} = convert(T, x) +singleprec(x::TwicePrecisionN{T}) where {T} = convert(T, x) """ From a1b998b574ff51ed633779bd5413c6684bd61fdf Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 12 Dec 2017 14:35:07 -0800 Subject: [PATCH 3/4] inlining --- src/KahanSummation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 26c45dd..071b80d 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -42,12 +42,12 @@ struct TwicePrecisionN{T} end -function plus_kbn(x::T, y::T) where {T} +@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, lo) end -function plus_kbn(x::T, y::TwicePrecisionN{T}) where {T} +@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 @@ -56,9 +56,9 @@ function plus_kbn(x::T, y::TwicePrecisionN{T}) where {T} end TwicePrecisionN(hi, nlo) end -plus_kbn(x::TwicePrecisionN{T}, y::T) where {T} = plus_kbn(y, x) +@inline plus_kbn(x::TwicePrecisionN{T}, y::T) where {T} = plus_kbn(y, x) -function plus_kbn(x::TwicePrecisionN{T}, y::TwicePrecisionN{T}) where {T} +@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 From e210c0704fe84b158b5dd53020631b4c47754976 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 12 Dec 2017 15:44:52 -0800 Subject: [PATCH 4/4] a few tweaks, prep for Base PR --- src/KahanSummation.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/KahanSummation.jl b/src/KahanSummation.jl index 071b80d..42fc7aa 100644 --- a/src/KahanSummation.jl +++ b/src/KahanSummation.jl @@ -45,7 +45,7 @@ end @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, lo) + TwicePrecisionN(hi, nlo) end @inline function plus_kbn(x::T, y::TwicePrecisionN{T}) where {T} hi = x + y.hi @@ -68,15 +68,20 @@ end TwicePrecisionN(hi, nlo) end -Base.r_promote_type(::typeof(plus_kbn), ::Type{T}) where {T} = - TwicePrecisionN{T} - 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) -Base.mr_empty(f, ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T)) +@static if VERSION >= v"0.7.0-" + 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_single(f, ::typeof(plus_kbn), x) = TwicePrecisionN(x, zero(x)) +else + Base.r_promote_type(::typeof(plus_kbn), ::Type{T}) where {T} = + TwicePrecisionN{T} + Base.mr_empty(f, ::typeof(plus_kbn), T) = TwicePrecisionN(zero(T),zero(T)) +end singleprec(x::TwicePrecisionN{T}) where {T} = convert(T, x)