diff --git a/src/implementations/BigFloat.jl b/src/implementations/BigFloat.jl index 442abed..ba8d414 100644 --- a/src/implementations/BigFloat.jl +++ b/src/implementations/BigFloat.jl @@ -9,12 +9,20 @@ mutability(::Type{BigFloat}) = IsMutable() -# Copied from `deepcopy_internal` implementation in Julia: -# https://github.com/JuliaLang/julia/blob/7d41d1eb610cad490cbaece8887f9bbd2a775021/base/mpfr.jl#L1041-L1050 -function mutable_copy(x::BigFloat) - d = x._d - d′ = GC.@preserve d unsafe_string(pointer(d), sizeof(d)) # creates a definitely-new String - return Base.MPFR._BigFloat(x.prec, x.sign, x.exp, d′) +# These methods are copied from `deepcopy_internal` in `base/mpfr.jl`. We don't +# use `mutable_copy(x) = deepcopy(x)` because this creates an empty `IdDict()` +# which costs some extra allocations. We don't need the IdDict case because we +# never call `mutable_copy` recursively. +@static if VERSION >= v"1.12.0-DEV.1343" + mutable_copy(x::BigFloat) = Base.MPFR._BigFloat(copy(getfield(x, :d))) +else + function mutable_copy(x::BigFloat) + d = x._d + GC.@preserve d begin + d′ = unsafe_string(pointer(d), sizeof(d)) + return Base.MPFR._BigFloat(x.prec, x.sign, x.exp, d′) + end + end end const _MPFRRoundingMode = Base.MPFR.MPFRRoundingMode @@ -297,12 +305,12 @@ function operate_to!( end struct DotBuffer{F<:Real} - compensation::F - summation_temp::F - multiplication_temp::F - inner_temp::F + c::F + t::F + input::F + tmp::F - DotBuffer{F}() where {F<:Real} = new{F}(ntuple(i -> F(), Val{4}())...) + DotBuffer{F}() where {F<:Real} = new{F}(zero(F), zero(F), zero(F), zero(F)) end function buffer_for( @@ -366,87 +374,61 @@ end # # function KahanBabushkaNeumaierSum(input) # sum = 0.0 -# # # A running compensation for lost low-order bits. # c = 0.0 -# # for i ∈ eachindex(input) # t = sum + input[i] -# # if abs(input[i]) ≤ abs(sum) -# c += (sum - t) + input[i] +# tmp = (sum - t) + input[i] # else -# c += (input[i] - t) + sum +# tmp = (input[i] - t) + sum # end -# +# c += tmp # sum = t # end -# -# # The result, with the correction only applied once in the very -# # end. +# # The result, with the correction only applied once in the very end. # sum + c # end -function buffered_operate_to!( - buf::DotBuffer{F}, - sum::F, - ::typeof(LinearAlgebra.dot), - x::AbstractVector{F}, - y::AbstractVector{F}, -) where {F<:BigFloat} - set! = (o, i) -> operate_to!(o, copy, i) - local swap! = function (x::BigFloat, y::BigFloat) - ccall((:mpfr_swap, :libmpfr), Cvoid, (Ref{BigFloat}, Ref{BigFloat}), x, y) - return nothing +# Returns abs(x) <= abs(y) without allocating. +function _abs_lte_abs(x::BigFloat, y::BigFloat) + x_is_neg, y_is_neg = signbit(x), signbit(y) + if x_is_neg != y_is_neg + operate!(-, x) end - - # Returns abs(x) <= abs(y) without allocating. - local abs_lte_abs = function (x::F, y::F) - local x_is_neg = signbit(x) - local y_is_neg = signbit(y) - - local x_neg = x_is_neg != y_is_neg - - x_neg && operate!(-, x) - - local ret = if y_is_neg - y <= x - else - x <= y - end - - x_neg && operate!(-, x) - - return ret + ret = y_is_neg ? y <= x : x <= y + if x_is_neg != y_is_neg + operate!(-, x) end + return ret +end - operate!(zero, sum) - operate!(zero, buf.compensation) - - for i in 0:(length(x)-1) - set!(buf.multiplication_temp, x[begin+i]) - operate!(*, buf.multiplication_temp, y[begin+i]) - - operate!(zero, buf.summation_temp) - operate_to!(buf.summation_temp, +, buf.multiplication_temp, sum) - - if abs_lte_abs(buf.multiplication_temp, sum) - set!(buf.inner_temp, sum) - operate!(-, buf.inner_temp, buf.summation_temp) - operate!(+, buf.inner_temp, buf.multiplication_temp) +function buffered_operate_to!( + buf::DotBuffer{BigFloat}, + sum::BigFloat, + ::typeof(LinearAlgebra.dot), + x::AbstractVector{BigFloat}, + y::AbstractVector{BigFloat}, +) # See pseudocode description + operate!(zero, sum) # sum = 0 + operate!(zero, buf.c) # c = 0 + for (xi, yi) in zip(x, y) # for i in eachindex(input) + operate_to!(buf.input, copy, xi) # input = x[i] + operate!(*, buf.input, yi) # input = x[i] * y[i] + operate_to!(buf.t, +, sum, buf.input) # t = sum + input + if _abs_lte_abs(buf.input, sum) # if |input| < |sum| + operate_to!(buf.tmp, copy, sum) # tmp = sum + operate!(-, buf.tmp, buf.t) # tmp = sum - t + operate!(+, buf.tmp, buf.input) # tmp = (sum - t) + input else - set!(buf.inner_temp, buf.multiplication_temp) - operate!(-, buf.inner_temp, buf.summation_temp) - operate!(+, buf.inner_temp, sum) + operate_to!(buf.tmp, copy, buf.input) # tmp = input + operate!(-, buf.tmp, buf.t) # tmp = input - t + operate!(+, buf.tmp, sum) # tmp = (input - t) + sum end - - operate!(+, buf.compensation, buf.inner_temp) - - swap!(sum, buf.summation_temp) + operate!(+, buf.c, buf.tmp) # c += tmp + operate_to!(sum, copy, buf.t) # sum = t end - - operate!(+, sum, buf.compensation) - + operate!(+, sum, buf.c) # sum += c return sum end diff --git a/test/bigfloat_dot.jl b/test/bigfloat_dot.jl index 37de53a..e783610 100644 --- a/test/bigfloat_dot.jl +++ b/test/bigfloat_dot.jl @@ -4,71 +4,6 @@ # v.2.0. If a copy of the MPL was not distributed with this file, You can obtain # one at http://mozilla.org/MPL/2.0/. -backup_bigfloats(v::AbstractVector{BigFloat}) = map(MA.copy_if_mutable, v) - -absolute_error(accurate::Real, approximate::Real) = abs(accurate - approximate) - -function relative_error(accurate::Real, approximate::Real) - return absolute_error(accurate, approximate) / abs(accurate) -end - -function dotter(x::V, y::V) where {V<:AbstractVector{<:Real}} - let x = x, y = y - () -> LinearAlgebra.dot(x, y) - end -end - -function reference_dot(x::V, y::V) where {F<:Real,V<:AbstractVector{F}} - return setprecision(dotter(x, y), F, 8 * precision(F)) -end - -function dot_test_relative_error(x::V, y::V) where {V<:AbstractVector{BigFloat}} - buf = MA.buffer_for(LinearAlgebra.dot, V, V) - - input = (x, y) - backup = map(backup_bigfloats, input) - - output = BigFloat() - - MA.buffered_operate_to!!(buf, output, LinearAlgebra.dot, input...) - - @test input == backup - - return relative_error(reference_dot(input...), output) -end - -subtracter(s::Real) = - let s = s - x -> x - s - end - -our_rand(n::Int, bias::Real) = map(subtracter(bias), rand(BigFloat, n)) - -function rand_dot_rel_err(size::Int, bias::Real) - x = our_rand(size, bias) - y = our_rand(size, bias) - return dot_test_relative_error(x, y) -end - -function max_rand_dot_rel_err(size::Int, bias::Real, iter_cnt::Int) - max_rel_err = zero(BigFloat) - for i in 1:iter_cnt - rel_err = rand_dot_rel_err(size, bias) - <(max_rel_err, rel_err) && (max_rel_err = rel_err) - end - return max_rel_err -end - -function max_rand_dot_ulps(size::Int, bias::Real, iter_cnt::Int) - return max_rand_dot_rel_err(size, bias, iter_cnt) / eps(BigFloat) -end - -function ulper(size::Int, bias::Real, iter_cnt::Int) - let s = size, b = bias, c = iter_cnt - () -> max_rand_dot_ulps(s, b, c) - end -end - @testset "prec:$prec size:$size bias:$bias" for (prec, size, bias) in Iterators.product( # These precisions (in bits) are most probably smaller than what @@ -78,11 +13,9 @@ end # precision (except when vector lengths are really huge with # respect to the precision). (32, 64), - # Compensated summation should be accurate even for very large # input vectors, so test that. (10000,), - # The zero "bias" signifies that the input will be entirely # nonnegative (drawn from the interval [0, 1]), while a positive # bias shifts that interval towards negative infinity. We want to @@ -91,8 +24,36 @@ end # no guarantee on the relative error in that case. (0.0, 2^-2, 2^-2 + 2^-3 + 2^-4), ) - iter_cnt = 10 - err = setprecision(ulper(size, bias, iter_cnt), BigFloat, prec) + err = setprecision(BigFloat, prec) do + maximum_relative_error = mapreduce(max, 1:10) do _ + # Generate some random vectors for dot(x, y) input. + x = rand(BigFloat, size) .- bias + y = rand(BigFloat, size) .- bias + # Copy x and y so that we can check we haven't mutated them after + # the fact. + old_x, old_y = MA.copy_if_mutable(x), MA.copy_if_mutable(y) + # Compute output = dot(x, y) + buf = MA.buffer_for( + LinearAlgebra.dot, + Vector{BigFloat}, + Vector{BigFloat}, + ) + output = BigFloat() + MA.buffered_operate_to!!(buf, output, LinearAlgebra.dot, x, y) + # Check that we haven't mutated x or y + @test old_x == x + @test old_y == y + # Compute dot(x, y) in larger precision. This will be used to + # compare with our `dot`. + accurate = setprecision(BigFloat, 8 * precision(BigFloat)) do + return LinearAlgebra.dot(x, y) + end + # Compute the relative error + return abs(accurate - output) / abs(accurate) + end + # Return estimate for ULP + return maximum_relative_error / eps(BigFloat) + end @test 0 <= err < 1 end