From 19525fc3af042aa959ea128c40c930a43b5f9202 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hamza=20Yusuf=20=C3=87ak=C4=B1r?= Date: Tue, 12 Dec 2023 23:58:27 +0100 Subject: [PATCH 1/5] Avoid mpfr `ccall`s for more types --- base/mpfr.jl | 145 +++++++++++++++++++++++++++++++++++---------------- test/mpfr.jl | 1 + 2 files changed, 101 insertions(+), 45 deletions(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index 276fd430ff1e0..a36938cd72ea7 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -237,10 +237,48 @@ function _duplicate(x::BigFloat) return z end +function is_bigfloat_special(x::Real) + # x is special if it is a zero, an Inf, or a NaN + return iszero(x) || !isfinite(x) +end + +function BigFloat_with_special_values(x::Real, precision::Integer=DEFAULT_PRECISION[]) + z = BigFloat(; precision=precision) + z.sign = 1-2*signbit(x) + + if iszero(x) + z.exp = mpfr_special_exponent_zero + elseif isinf(x) + z.exp = mpfr_special_exponent_inf + else + z.exp = mpfr_special_exponent_nan + end + + return z +end + +function store_significand_in_limbs!(z::BigFloat, significand::UInt64) + nlimbs = (z.prec + 8*Core.sizeof(Limb) - 1) ÷ (8*Core.sizeof(Limb)) + + # Limb is a CLong which is a UInt32 on Windows which makes this more complicated and slower. + if Limb === UInt64 + for i in 1:nlimbs-1 + unsafe_store!(z.d, 0x0, i) + end + unsafe_store!(z.d, significand, nlimbs) + else + for i in 1:nlimbs-2 + unsafe_store!(z.d, 0x0, i) + end + unsafe_store!(z.d, significand % UInt32, nlimbs-1) + unsafe_store!(z.d, (significand >> 32) % UInt32, nlimbs) + end +end + # convert to BigFloat for (fJ, fC) in ((:si,:Clong), (:ui,:Culong)) @eval begin - function BigFloat(x::($fC), r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) + function BigFloat_mpfr_fallback(x::($fC), r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) z = BigFloat(;precision=precision) ccall(($(string(:mpfr_set_,fJ)), libmpfr), Int32, (Ref{BigFloat}, $fC, MPFRRoundingMode), z, x, r) return z @@ -248,47 +286,61 @@ for (fJ, fC) in ((:si,:Clong), (:ui,:Culong)) end end -function BigFloat(x::Float64, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) - z = BigFloat(;precision) +function BigFloat(x::Base.BitInteger64, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) + if is_bigfloat_special(x) + return BigFloat_with_special_values(x, precision)::BigFloat + end + + unsigned_x = unsigned(abs(x)) + exp_x = (sizeof(x) << 3) - leading_zeros(unsigned_x) + + requiredprecision = exp_x - 1 + + if requiredprecision < precision + z = BigFloat(; precision=precision)::BigFloat + z.exp = exp_x + z.sign = 1 - 2*signbit(x) + significand = UInt64(unsigned_x) << leading_zeros(UInt64(unsigned_x)) + store_significand_in_limbs!(z, significand) + + return z + else + return BigFloat_mpfr_fallback(x, r; precision=precision)::BigFloat + end +end + +function BigFloat_mpfr_fallback(x::Float64, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) + z = BigFloat(;precision=precision) + ccall((:mpfr_set_d, libmpfr), Int32, (Ref{BigFloat}, Float64, MPFRRoundingMode), z, x, r) + return z +end + +function BigFloat_mpfr_fallback(x::Union{Float16, Float32}, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) + BigFloat_mpfr_fallback(Float64(x), r; precision=precision) +end + +function BigFloat(x::Base.IEEEFloat, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) + if is_bigfloat_special(x) + return BigFloat_with_special_values(x, precision)::BigFloat + end + # punt on the hard case where we might have to deal with rounding # we could use this path in all cases, but mpfr_set_d has a lot of overhead. - if precision <= Base.significand_bits(Float64) - ccall((:mpfr_set_d, libmpfr), Int32, (Ref{BigFloat}, Float64, MPFRRoundingMode), z, x, r) - if isnan(x) && signbit(x) != signbit(z) - z.sign = -z.sign - end - return z + if precision <= Base.significand_bits(typeof(x)) + return BigFloat_mpfr_fallback(x, r; precision=precision) end + + z = BigFloat(;precision) z.sign = 1-2*signbit(x) - if iszero(x) || !isfinite(x) - if isinf(x) - z.exp = mpfr_special_exponent_inf - elseif isnan(x) - z.exp = mpfr_special_exponent_nan - else - z.exp = mpfr_special_exponent_zero - end - return z - end z.exp = 1 + exponent(x) + # BigFloat doesn't have an implicit bit - val = reinterpret(UInt64, significand(x))<<11 | typemin(Int64) - nlimbs = (precision + 8*Core.sizeof(Limb) - 1) ÷ (8*Core.sizeof(Limb)) + val = UInt64(reinterpret(Base.uinttype(typeof(x)), significand(x))) + val <<= (sizeof(UInt64) << 3 - sizeof(x) << 3) + Base.exponent_bits(typeof(x)) + val |= typemin(Int64) - # Limb is a CLong which is a UInt32 on windows (thank M$) which makes this more complicated and slower. - if Limb === UInt64 - for i in 1:nlimbs-1 - unsafe_store!(z.d, 0x0, i) - end - unsafe_store!(z.d, val, nlimbs) - else - for i in 1:nlimbs-2 - unsafe_store!(z.d, 0x0, i) - end - unsafe_store!(z.d, val % UInt32, nlimbs-1) - unsafe_store!(z.d, (val >> 32) % UInt32, nlimbs) - end - z + store_significand_in_limbs!(z, val) + return z end function BigFloat(x::BigInt, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) @@ -302,18 +354,21 @@ BigFloat(x::Integer; precision::Integer=DEFAULT_PRECISION[]) = BigFloat(x::Integer, r::MPFRRoundingMode; precision::Integer=DEFAULT_PRECISION[]) = BigFloat(BigInt(x)::BigInt, r; precision=precision) -BigFloat(x::Union{Bool,Int8,Int16,Int32}, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) = - BigFloat(convert(Clong, x), r; precision=precision) -BigFloat(x::Union{UInt8,UInt16,UInt32}, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) = - BigFloat(convert(Culong, x), r; precision=precision) - -BigFloat(x::Union{Float16,Float32}, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) = - BigFloat(Float64(x), r; precision=precision) +BigFloat(x::Bool, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) = + BigFloat(convert(Clong, x), r; precision=precision) function BigFloat(x::Rational, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) - setprecision(BigFloat, precision) do - setrounding_raw(BigFloat, r) do - BigFloat(numerator(x))::BigFloat / BigFloat(denominator(x))::BigFloat + if is_bigfloat_special(x) + BigFloat_with_special_values(x, precision)::BigFloat + elseif ispow2(denominator(x)) + z = BigFloat(numerator(x), r; precision)::BigFloat + z.exp -= trailing_zeros(denominator(x)) + z + else + setprecision(BigFloat, precision) do + setrounding_raw(BigFloat, r) do + BigFloat(numerator(x))::BigFloat / BigFloat(denominator(x))::BigFloat + end end end end diff --git a/test/mpfr.jl b/test/mpfr.jl index a0dd15d97f70c..0a2871fba93a6 100644 --- a/test/mpfr.jl +++ b/test/mpfr.jl @@ -27,6 +27,7 @@ import Base.MPFR @test typeof(BigFloat(typemax(UInt64))) == BigFloat @test typeof(BigFloat(typemax(UInt128))) == BigFloat + @test typeof(BigFloat(floatmax(Float16))) == BigFloat @test typeof(BigFloat(floatmax(Float32))) == BigFloat @test typeof(BigFloat(floatmax(Float64))) == BigFloat From 40db873522af4f6dde65f29a26a87f5af2f726a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hamza=20Yusuf=20=C3=87ak=C4=B1r?= Date: Thu, 14 Dec 2023 00:10:57 +0100 Subject: [PATCH 2/5] Fix fallback on Windows --- base/mpfr.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index a36938cd72ea7..90cb076fddfa9 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -305,7 +305,10 @@ function BigFloat(x::Base.BitInteger64, r::MPFRRoundingMode=ROUNDING_MODE[]; pre return z else - return BigFloat_mpfr_fallback(x, r; precision=precision)::BigFloat + # If Clong is 32 bit, MPFR does not have a method for (U)Int64 on Windows + return sizeof(x) <= sizeof(Clong) ? + BigFloat_mpfr_fallback(x, r; precision)::BigFloat : + BigFloat(BigInt(x), r; precision) end end From 0e69aa31aefaa6b98d37bca22caad90b79f7bf19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hamza=20Yusuf=20=C3=87ak=C4=B1r?= Date: Thu, 14 Dec 2023 01:21:56 +0100 Subject: [PATCH 3/5] Remove trailing whitespace --- base/mpfr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index 90cb076fddfa9..d21cc18df83d0 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -306,7 +306,7 @@ function BigFloat(x::Base.BitInteger64, r::MPFRRoundingMode=ROUNDING_MODE[]; pre return z else # If Clong is 32 bit, MPFR does not have a method for (U)Int64 on Windows - return sizeof(x) <= sizeof(Clong) ? + return sizeof(x) <= sizeof(Clong) ? BigFloat_mpfr_fallback(x, r; precision)::BigFloat : BigFloat(BigInt(x), r; precision) end From fbe3e0d66704b7df168450d72e3271f619a7f8e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hamza=20Yusuf=20=C3=87ak=C4=B1r?= Date: Sat, 13 Jan 2024 01:36:47 +0100 Subject: [PATCH 4/5] Stop creating a `BigFloat` for denominator --- base/mpfr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index d21cc18df83d0..d8e8fedf3aed8 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -370,7 +370,7 @@ function BigFloat(x::Rational, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::I else setprecision(BigFloat, precision) do setrounding_raw(BigFloat, r) do - BigFloat(numerator(x))::BigFloat / BigFloat(denominator(x))::BigFloat + BigFloat(numerator(x))::BigFloat / denominator(x) end end end From 770900211a157eabf44e1cf8f5bb74a67fbba4a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hamza=20Yusuf=20=C3=87ak=C4=B1r?= Date: Sat, 13 Jan 2024 02:06:28 +0100 Subject: [PATCH 5/5] Multiply instead of bitshift when computing bitsize --- base/mpfr.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index d8e8fedf3aed8..4da7e52f40b8e 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -292,7 +292,7 @@ function BigFloat(x::Base.BitInteger64, r::MPFRRoundingMode=ROUNDING_MODE[]; pre end unsigned_x = unsigned(abs(x)) - exp_x = (sizeof(x) << 3) - leading_zeros(unsigned_x) + exp_x = 8 * sizeof(x) - leading_zeros(unsigned_x) requiredprecision = exp_x - 1 @@ -339,7 +339,7 @@ function BigFloat(x::Base.IEEEFloat, r::MPFRRoundingMode=ROUNDING_MODE[]; precis # BigFloat doesn't have an implicit bit val = UInt64(reinterpret(Base.uinttype(typeof(x)), significand(x))) - val <<= (sizeof(UInt64) << 3 - sizeof(x) << 3) + Base.exponent_bits(typeof(x)) + val <<= (8 * (sizeof(UInt64) - sizeof(x))) + Base.exponent_bits(typeof(x)) val |= typemin(Int64) store_significand_in_limbs!(z, val)