diff --git a/base/Base.jl b/base/Base.jl index 0ca13265adc4f1..46bca3120268d2 100644 --- a/base/Base.jl +++ b/base/Base.jl @@ -253,6 +253,7 @@ include("rounding.jl") include("div.jl") include("rawbigints.jl") include("float.jl") +include("rational_to_float.jl") include("twiceprecision.jl") include("complex.jl") include("rational.jl") diff --git a/base/mpfr.jl b/base/mpfr.jl index 276fd430ff1e09..c92768392a46ce 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -19,7 +19,8 @@ import isone, big, _string_n, decompose, minmax, sinpi, cospi, sincospi, tanpi, sind, cosd, tand, asind, acosd, atand, uinttype, exponent_max, exponent_min, ieee754_representation, significand_mask, - RawBigIntRoundingIncrementHelper, truncated, RawBigInt + RawBigIntRoundingIncrementHelper, truncated, RawBigInt, unsafe_rational, + RationalToFloat, rational_to_floating_point using .Base.Libc @@ -310,12 +311,51 @@ BigFloat(x::Union{UInt8,UInt16,UInt32}, r::MPFRRoundingMode=ROUNDING_MODE[]; pre BigFloat(x::Union{Float16,Float32}, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) = BigFloat(Float64(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 +function set_2exp!(z::BigFloat, n::BigInt, exp::Int, rm::MPFRRoundingMode) + ccall( + (:mpfr_set_z_2exp, libmpfr), + Int32, + (Ref{BigFloat}, Ref{BigInt}, Int, MPFRRoundingMode), + z, n, exp, rm, + ) + nothing +end + +function RationalToFloat.to_floating_point_impl(::Type{BigFloat}, ::Type{BigInt}, num, den, romo, prec) + num_is_zero = iszero(num) + den_is_zero = iszero(den) + s = Int8(sign(num)) + sb = signbit(s) + is_zero = num_is_zero & !den_is_zero + is_inf = !num_is_zero & den_is_zero + is_regular = !num_is_zero & !den_is_zero + + if is_regular + let rtfc = RationalToFloat.to_float_components + c = rtfc(BigInt, num, den, prec, nothing, romo, sb) + ret = BigFloat(precision = prec) + mpfr_romo = convert(MPFRRoundingMode, romo) + set_2exp!(ret, s * c.integral_significand, Int(c.exponent - prec + true), mpfr_romo) + ret end - end + else + if is_zero + BigFloat(false, MPFRRoundToZero, precision = prec) + elseif is_inf + BigFloat(s * Inf, MPFRRoundToZero, precision = prec) + else + BigFloat(precision = prec) + end + end::BigFloat +end + +function BigFloat(x::Rational, r::RoundingMode; precision::Integer = DEFAULT_PRECISION[]) + rational_to_floating_point(BigFloat, x, r, precision) +end + +function BigFloat(x::Rational, r::MPFRRoundingMode = ROUNDING_MODE[]; + precision::Integer = DEFAULT_PRECISION[]) + rational_to_floating_point(BigFloat, x, r, precision) end function tryparse(::Type{BigFloat}, s::AbstractString; base::Integer=0, precision::Integer=DEFAULT_PRECISION[], rounding::MPFRRoundingMode=ROUNDING_MODE[]) diff --git a/base/rational.jl b/base/rational.jl index 1b060384d95b20..1102db2bdc75f8 100644 --- a/base/rational.jl +++ b/base/rational.jl @@ -140,10 +140,21 @@ Bool(x::Rational) = x==0 ? false : x==1 ? true : (::Type{T})(x::Rational) where {T<:Integer} = (isinteger(x) ? convert(T, x.num)::T : throw(InexactError(nameof(T), T, x))) -AbstractFloat(x::Rational) = (float(x.num)/float(x.den))::AbstractFloat -function (::Type{T})(x::Rational{S}) where T<:AbstractFloat where S - P = promote_type(T,S) - convert(T, convert(P,x.num)/convert(P,x.den))::T +function numerator_denominator_promoted(x) + y = unsafe_rational(numerator(x), denominator(x)) + (numerator(y), denominator(y)) +end + +function rational_to_floating_point(::Type{F}, x, rm = RoundNearest, prec = precision(F)) where {F} + nd = numerator_denominator_promoted(x) + RationalToFloat.to_floating_point(F, nd..., rm, prec)::F +end + +(::Type{F})(x::Rational) where {F<:AbstractFloat} = rational_to_floating_point(F, x)::F + +function AbstractFloat(x::Q) where {Q<:Rational} + T = float(Q) + T(x)::T::AbstractFloat end function Rational{T}(x::AbstractFloat) where T<:Integer diff --git a/base/rational_to_float.jl b/base/rational_to_float.jl new file mode 100644 index 00000000000000..3751b2095022f2 --- /dev/null +++ b/base/rational_to_float.jl @@ -0,0 +1,226 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +module RationalToFloat + +const Rnd = Base.Rounding + +# Performance optimization. Unlike raw `<<` or `>>>`, this is supposed +# to compile to a single instruction, because the semantics correspond +# to what hardware usually provides. +function machine_shift(shift::S, a::T, b) where {S,T<:Base.BitInteger} + @inline begin + mask = 8*sizeof(T) - 1 + c = b & mask + shift(a, c) + end +end + +machine_shift(::S, a::Bool, ::Any) where {S} = error("unsupported") + +# Fallback for `BigInt` etc. +machine_shift(shift::S, a, b) where {S} = shift(a, b) + +# Arguments are positive integers. +function div_significand_with_remainder(num, den, minimum_significand_size) + clamped = x -> max(zero(x), x)::typeof(x) + bw = Base.top_set_bit # bit width + shift = clamped(minimum_significand_size + bw(den) - bw(num) + 0x2) + t = machine_shift(<<, num, shift) + (divrem(t, den, RoundToZero)..., shift) +end + +# `divrem(n, 1<>>, n, k) + tmp = machine_shift(<<, quo, k) + rem = n - tmp + (quo, rem) +end + +function to_float_components_impl(num, den, precision, max_subnormal_exp) + # `+1` because we need an extra, "round", bit for some rounding modes. + # + # TODO: as a performance optimization, only do this when required + # by the rounding mode + prec_p_1 = precision + true + + (quo0, rem0, shift) = div_significand_with_remainder(num, den, prec_p_1) + width = Base.top_set_bit(quo0) + excess_width = width - prec_p_1 + exp = width - shift - true + + exp_underflow = if isnothing(max_subnormal_exp) + zero(exp) + else + let d = max_subnormal_exp - exp, T = typeof(d), z = zero(d)::T + (signbit(d) ? z : d + true)::T + end + end + + (quo1, rem1) = divrem_2(quo0, exp_underflow + excess_width) + (quo2, rem2) = (quo1 >>> true, quo1 % Bool) + + round_bit = rem2 + sticky_bit = !iszero(rem1) | !iszero(rem0) + + (; integral_significand = quo2, exponent = exp, round_bit, sticky_bit) +end + +struct RoundingIncrementHelper + final_bit::Bool + round_bit::Bool + sticky_bit::Bool +end + +(h::RoundingIncrementHelper)(::Rnd.FinalBit) = h.final_bit +(h::RoundingIncrementHelper)(::Rnd.RoundBit) = h.round_bit +(h::RoundingIncrementHelper)(::Rnd.StickyBit) = h.sticky_bit + +function to_float_components_rounded(num, den, precision, max_subnormal_exp, romo, sign_bit) + # Helps ensure type stability + fun = function(f::F, v) where {F} + r = f(v) + T = typeof(r) + (T(v)::T, r)::NTuple{2,T} + end + + fun_sig = x -> x >>> true + fun_exp = x -> x + true + overflows = (x, p) -> x == machine_shift(<<, one(x), p) + + t = to_float_components_impl(num, den, precision, max_subnormal_exp) + raw_significand = t.integral_significand + rh = RoundingIncrementHelper(raw_significand % Bool, t.round_bit, t.sticky_bit) + incr = Rnd.correct_rounding_requires_increment(rh, romo, sign_bit) + rounded = raw_significand + incr + (integral_significand, exponent) = let exp = t.exponent, s01, e01, se0, se1, T + s01 = fun(fun_sig, rounded) + e01 = fun(fun_exp, exp) + se0 = (first(s01), first(e01)) + se1 = ( last(s01), last(e01)) + T = typeof(se0) + (overflows(rounded, precision) ? se1 : se0)::T + end + (; integral_significand, exponent) +end + +function to_float_components(::Type{T}, num, den, precision, max_subnormal_exp, romo, sb) where {T} + to_float_components_rounded(abs(T(num)), den, precision, max_subnormal_exp, romo, sb) +end + +function to_floating_point_fallback(::Type{T}, ::Type{S}, num, den, rm, prec) where {T,S} + num_is_zero = iszero(num) + den_is_zero = iszero(den) + sb = signbit(num) + is_zero = num_is_zero & !den_is_zero + is_inf = !num_is_zero & den_is_zero + is_regular = !num_is_zero & !den_is_zero + if is_regular + let + c = to_float_components(S, num, den, prec, nothing, rm, sb) + exp = c.exponent + signif = T(c.integral_significand)::T + let x = ldexp(signif, exp - prec + true)::T + sb ? -x : x + end::T + end + else + if is_zero + zero(T)::T + elseif is_inf + T(Inf)::T + else + T(NaN)::T + end + end::T +end + +function to_floating_point_impl(::Type{T}, ::Type{S}, num, den, rm, prec) where {T,S} + to_floating_point_fallback(T, S, num, den, rm, prec) +end + +function to_floating_point_impl(::Type{T}, ::Type{S}, num, den, rm, prec) where {T<:Base.IEEEFloat,S} + num_is_zero = iszero(num) + den_is_zero = iszero(den) + sb = signbit(num) + is_zero = num_is_zero & !den_is_zero + is_inf = !num_is_zero & den_is_zero + is_regular = !num_is_zero & !den_is_zero + (rm_is_to_zero, rm_is_from_zero) = if Rnd.rounds_to_nearest(rm) + (false, false) + else + let from = Rnd.rounds_away_from_zero(rm, sb) + (!from, from) + end + end::NTuple{2,Bool} + exp_max = Base.exponent_max(T) + exp_min = Base.exponent_min(T) + ieee_repr = Base.ieee754_representation + repr_zero = ieee_repr(T, sb, Val(:zero)) + repr_inf = ieee_repr(T, sb, Val(:inf)) + repr_nan = ieee_repr(T, sb, Val(:nan)) + U = typeof(repr_zero) + repr_zero::U + repr_inf::U + repr_nan::U + + ret_u = if is_regular + let + c = let e = exp_min - 1 + to_float_components(S, num, den, prec, e, rm, sb) + end + exp = c.exponent + exp_diff = exp - exp_min + is_normal = 0 ≤ exp_diff + exp_is_huge_p = exp_max < exp + exp_is_huge_n = signbit(exp_diff + prec) + rounds_to_inf = exp_is_huge_p & !rm_is_to_zero + rounds_to_zero = exp_is_huge_n & !rm_is_from_zero + + if !rounds_to_zero & !exp_is_huge_p + let signif = (c.integral_significand % U) & Base.significand_mask(T) + exp_field = (max(exp_diff, zero(exp_diff)) + is_normal) % U + ieee_repr(T, sb, exp_field, signif)::U + end + elseif rounds_to_zero + repr_zero + elseif rounds_to_inf + repr_inf + else + ieee_repr(T, sb, Val(:omega)) + end + end + else + if is_zero + repr_zero + elseif is_inf + repr_inf + else + repr_nan + end + end::U + + reinterpret(T, ret_u)::T +end + +# `BigInt` is a safe default. +to_float_promote_type(::Type{F}, ::Type{S}) where {F,S} = BigInt + +const BitIntegerOrBool = Union{Bool,Base.BitInteger} + +# As an optimization, use an integer type narrower than `BigInt` when possible. +function to_float_promote_type(::Type{F}, ::Type{S}) where {F<:Base.IEEEFloat,S<:BitIntegerOrBool} + Max = if sizeof(F) ≤ sizeof(S) + S + else + (S <: Signed) ? Base.inttype(F) : Base.uinttype(F) + end + widen(Max) +end + +function to_floating_point(::Type{F}, num::T, den::T, rm, prec) where {F,T} + S = to_float_promote_type(F, T) + to_floating_point_impl(F, S, num, den, rm, prec) +end + +end diff --git a/test/choosetests.jl b/test/choosetests.jl index beed4e15a58df3..e04eaec411ad2c 100644 --- a/test/choosetests.jl +++ b/test/choosetests.jl @@ -11,7 +11,7 @@ const TESTNAMES = [ "char", "strings", "triplequote", "unicode", "intrinsics", "dict", "hashing", "iobuffer", "staged", "offsetarray", "arrayops", "tuple", "reduce", "reducedim", "abstractarray", - "intfuncs", "simdloop", "vecelement", "rational", + "intfuncs", "simdloop", "vecelement", "rational", "rational_to_float", "bitarray", "copy", "math", "fastmath", "functional", "iterators", "operators", "ordering", "path", "ccall", "parse", "loading", "gmp", "sorting", "spawn", "backtrace", "exceptions", diff --git a/test/rational.jl b/test/rational.jl index 118f2e62892285..6c705f6164c16c 100644 --- a/test/rational.jl +++ b/test/rational.jl @@ -745,3 +745,303 @@ end @test rationalize(Int64, nextfloat(0.1) * im; tol=0) == precise_next * im @test rationalize(0.1im; tol=eps(0.1)) == rationalize(0.1im) end + +naive_rational_to_bigfloat(x) = BigFloat(numerator(x))/BigFloat(denominator(x)) + +function test_rational_to_float_systematic_with_naive_impl(::Type{F}, q) where {F} + no_alloc_rat = let + f(::Rational{Bool}) = true + f(x) = isbits(widen(x)) + f + end + want_no_alloc = no_alloc_rat(q) && isbitstype(widen(F)) + want_no_alloc && @test iszero(@allocated F(q)) + @test F(naive_rational_to_bigfloat(q)) == @inferred F(q) +end + +function test_rational_to_float_systematic_with_naive_impl(::Type{BigFloat}, q, rm, prec) + reference = BigFloat(naive_rational_to_bigfloat(q), rm, precision = prec) + @test precision(reference) == prec + t = let reference = reference, prec = prec + function(res) + @test res == reference + @test precision(res) == prec + end + end + + t(BigFloat(q, rm, precision = prec)) + + let + res = setprecision(BigFloat, prec) do + BigFloat(q, rm) + end + t(res) + end + + let + if rm isa RoundingMode + res = setrounding(BigFloat, rm) do + BigFloat(q, precision = prec) + end + t(res) + end + end +end + +function test_rational_to_float_systematic_with_naive_impl(::Type{BigFloat}, q, rm) + for prec ∈ Base.OneTo(10) + test_rational_to_float_systematic_with_naive_impl(BigFloat, q, rm, prec) + end +end + +function test_rational_to_float_systematic_with_naive_impl(::Type{BigFloat}, q) + M = Base.MPFR + rms = ( + RoundDown, RoundFromZero, RoundNearest, RoundToZero, RoundUp, + # XXX: these could be tested, but would require more effort to + # test, as MPFR doesn't support them: + # + # ```julia + # RoundNearestTiesAway, RoundNearestTiesUp, + # ``` + M.MPFRRoundDown, M.MPFRRoundUp, M.MPFRRoundFromZero, M.MPFRRoundToZero, + M.MPFRRoundNearest + ) + for rm ∈ rms + test_rational_to_float_systematic_with_naive_impl(BigFloat, q, rm) + end +end + +function min_no_promotion_impl(isless::F, a::T, b) where {F, T} + if isless(b, a) + T(b)::T + else + a + end +end +min_no_promotion(a::T, b) where {T} = min_no_promotion_impl(isless, a, b) +max_no_promotion(a::T, b) where {T} = min_no_promotion_impl(((l, r) -> isless(r, l)), a, b) + +safe_typemin(::Type{T}) where {T} = typemin(T)::T +safe_typemax(::Type{T}) where {T} = typemax(T)::T + +safe_typemin(::Type{BigInt}) = BigInt(-10000) +safe_typemax(::Type{BigInt}) = BigInt( 10000) + +function test_rational_to_float_systematic_with_naive(::Type{F}, ::Type{T}) where {F,T} + range_typemin = function(x) + l = safe_typemin(T) + u = safe_typemax(T) + b = min_no_promotion(u, l + x) + l:b + end + + range_typemax = function(x) + l = safe_typemin(T) + u = safe_typemax(T) + a = max_no_promotion(l, u - x) + a:u + end + + range_zero = function(x) + l = safe_typemin(T) + u = safe_typemax(T) + a = max_no_promotion(l, -x) + b = min_no_promotion(u, x) + a:b + end + + k = 20 + nums = Vector{T}() + dens = Vector{T}() + + append!(nums, range_typemin(k)) + + for x ∈ range_zero(k) + push!(nums, x) + (0 < x) && push!(dens, x) + end + + for x ∈ range_typemax(k) + push!(nums, x) + (0 < x) && push!(dens, x) + end + + sort!(nums) + sort!(dens) + unique!(nums) + unique!(dens) + + for num ∈ nums + for den ∈ dens + # XXX: shift `n` and `d` to vary the FP exponent more + test_rational_to_float_systematic_with_naive_impl(F, num//den) + end + end +end + +function test_rational_to_float_ldexp_fallback(::Type{F}, ::Type{T}) where {F,T} + l = max_no_promotion(safe_typemin(T), -10) + u = min_no_promotion(safe_typemax(T), 10) + prec = precision(F) + rm = RoundNearest + rtfp = Base.RationalToFloat.to_floating_point_fallback + @testset "num: $num" for num ∈ l:u + @testset "den: $den" for den ∈ true:u + @test F(num//den) == rtfp(F, BigInt, num, den, rm, prec) + end + end +end + +function test_rational_to_float_correctly_rounded_up_down(x, f, ::Val{:down}) + @test f ≤ x < nextfloat(f) +end + +function test_rational_to_float_correctly_rounded_up_down(x, f, ::Val{:up}) + @test prevfloat(f) < x ≤ f +end + +function test_rational_to_float_simplest_rms_up_down(::Type{F}, ::Type{T}, rm, sb, ::V) where {F,T,V<:Val} + range_typemax = function(x) + l = one(T)::T + u = safe_typemax(T) + a = max_no_promotion(l, u - x) + a:u + end + + range_zero = function(x) + u = safe_typemax(T) + b = min_no_promotion(u, x) + true:b + end + + k = 50 + nds = Vector{T}() + append!(nds, range_zero(k)) + append!(nds, range_typemax(k)) + sort!(nds) + unique!(nds) + + rtfp = Base.RationalToFloat.to_floating_point + prec = precision(F) + + @testset "num: $num" for num ∈ nds + @testset "den: $den" for den ∈ nds + x_pos = num // den + x = sb ? -x_pos : x_pos + f = rtfp(F, numerator(x), denominator(x), rm, prec) + test_rational_to_float_correctly_rounded_up_down(x, f, V()) + end + end +end + +function test_rational_to_float_simplest_rms(::Type{F}, ::Type{T}) where {F,T} + M = Base.MPFR + rms_up = (RoundUp, M.MPFRRoundUp) + rms_down = (RoundDown, M.MPFRRoundDown) + rms_fz = (RoundFromZero, M.MPFRRoundFromZero) + rms_tz = (RoundToZero, M.MPFRRoundToZero) + + if T <: Signed + @testset "negatives" begin + @testset "round down" begin + @testset "rm: $rm" for rm ∈ (rms_fz..., rms_down...) + test_rational_to_float_simplest_rms_up_down(F, T, rm, true, Val(:down)) + end + end + @testset "round up" begin + @testset "rm: $rm" for rm ∈ (rms_tz..., rms_up...) + test_rational_to_float_simplest_rms_up_down(F, T, rm, true, Val(:up)) + end + end + end + end + + @testset "positives" begin + @testset "round down" begin + @testset "rm: $rm" for rm ∈ (rms_tz..., rms_down...) + test_rational_to_float_simplest_rms_up_down(F, T, rm, false, Val(:down)) + end + end + @testset "round up" begin + @testset "rm: $rm" for rm ∈ (rms_fz..., rms_up...) + test_rational_to_float_simplest_rms_up_down(F, T, rm, false, Val(:up)) + end + end + end +end + +struct AbstractFloatFallback <: AbstractFloat + x::Float32 + AbstractFloatFallback(x::Float32) = new(x) +end +AbstractFloatFallback(n::Integer) = AbstractFloatFallback(Float32(n)) +Base.precision(::Type{AbstractFloatFallback}) = 10 +Base.ldexp(x::AbstractFloatFallback, n::Integer) = AbstractFloatFallback(ldexp(x.x, n)) +Base.:-(x::AbstractFloatFallback) = AbstractFloatFallback(-x.x) + +struct IntegerFallback <: Integer + n::Int + IntegerFallback(n::Int) = new(n) +end +Base.signbit(n::IntegerFallback) = signbit(n.n) +Base.BigInt(n::IntegerFallback) = BigInt(n.n) +Base.oneunit(n::IntegerFallback) = IntegerFallback(oneunit(n.n)) +Base.:+(l::IntegerFallback, r::IntegerFallback) = IntegerFallback(l.n + r.n) +Base.log2(n::IntegerFallback) = log2(n.n) + +@testset "rational number to floating-point number" begin + @testset "defined for unknown types" begin + Floats = (Float16, AbstractFloatFallback) + rationals = (3//1, Base.unsafe_rational(IntegerFallback(3), IntegerFallback(1))) + @testset "F: $F" for F ∈ Floats + @testset "q: $q" for q ∈ rationals + if F <: AbstractFloatFallback + @test F(q).x == 3 + else + @test F(q) == 3 + end + end + end + end + + Us = (UInt8, UInt16, UInt32, UInt64, UInt128) + Ss = (Int8, Int16, Int32, Int64, Int128) + Ts = (Bool, Us..., Ss..., BigInt) + Fs = (Float16, Float32, Float64, BigFloat) + + @testset "the simplest rounding modes; up, down" begin + @testset "F: $F" for F ∈ Fs + @testset "T: $T" for T ∈ Ts + test_rational_to_float_simplest_rms(F, T) + end + end + end + + @testset "with naive conversion via `BigFloat`" begin + @testset "F: $F" for F ∈ Fs + @testset "T: $T" for T ∈ Ts + test_rational_to_float_systematic_with_naive(F, T) + end + end + end + + @testset "`ldexp`-based fallback" begin + @testset "F: $F" for F ∈ Fs + @testset "T: $T" for T ∈ Ts + test_rational_to_float_ldexp_fallback(F, T) + end + end + end + + @testset "special" begin + @test Float16(9//70000) === Float16(0.0001286) + @test Float32(16777216//16777217) < 1 + @test float(-9223372036854775808//9223372036854775807) == -1 + + # issue #52394 + @test Float16(10^8 // (10^9 + 1)) === convert(Float16, 10^8 // (10^9 + 1)) === Float16(0.1) + @test Float16((typemax(UInt128)-0x01) // typemax(UInt128)) === Float16(1.0) + @test Float32((typemax(UInt128)-0x01) // typemax(UInt128)) === Float32(1.0) + end +end diff --git a/test/rational_to_float.jl b/test/rational_to_float.jl new file mode 100644 index 00000000000000..9c5226b6616db7 --- /dev/null +++ b/test/rational_to_float.jl @@ -0,0 +1,160 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +using Test + +function div_significand_with_remainder_naive(num, den, minimum_significand_size) + dr = let den = den + n -> divrem(n, den, RoundToZero) + end + + shift = 0 + (quo, rem) = dr(num << shift) + + # bit width + bw = Base.top_set_bit + + while bw(quo) < minimum_significand_size + shift += true + (quo, rem) = dr(num << shift) + end + + # suboptimal but allowed results + (q1, r1) = dr(num << (shift + 1)) + (q2, r2) = dr(num << (shift + 2)) + (q3, r3) = dr(num << (shift + 3)) + + ( + (quo, rem, shift), # optimal result + (q1, r1, shift + 1), + (q2, r2, shift + 2), + (q3, r3, shift + 3) + ) +end + +function test_div_significand_with_remainder(n, d, mss) + (0 < n) || error("unsupported") + (0 < d) || error("unsupported") + (0 < mss) || error("unsupported") + bw = Base.top_set_bit + wid = x -> widen(widen(x)) + t = let n = wid(n), d = wid(d) + function(q_narrow, r_narrow, sh) + q = wid(q_narrow) + r = wid(r_narrow) + @test q + r/d ≈ n*(2.0^sh)/d + @test floor(Int, log2(abs(q))) == exponent(q) == bw(q) - 1 == sh + exponent(n/d) + end + end + dswr = Base.RationalToFloat.div_significand_with_remainder + a = @inferred dswr(n, d, mss) + (a0, a1, a2, a3) = div_significand_with_remainder_naive(n, d, mss) + @test a ∈ (a0, a1, a2, a3) + t(a...) + t(a0...) + t(a1...) + t(a2...) + t(a3...) + nothing +end + +function test_div_significand_with_remainder_iters(::Type{T}, nds, msss) where {T} + @testset "mss: $mss" for mss ∈ msss + @testset "num: $num" for num ∈ nds + @testset "den: $den" for den ∈ nds + # XXX: shift `n` and `d` to vary the FP exponent more + n = T(num)::T + d = T(den)::T + test_div_significand_with_remainder(widen(n), d, mss) + end + end + end +end + +max_mss_for_type(::Type{T}) where {T<:Unsigned} = sizeof(T)*8 - 3 +max_mss_for_type(::Type{T}) where {T<:Signed} = sizeof(T)*8 - 4 +max_mss_for_type(::Type{BigInt}) = 200 # just a big value + +max_nd_for_type(::Type{T}) where {T} = typemax(T) +max_nd_for_type(::Type{BigInt}) = BigInt(100000) # just a big value + +function min_no_promotion_impl(isless::F, a::T, b) where {F, T} + if isless(b, a) + T(b)::T + else + a + end +end +min_no_promotion(a::T, b) where {T} = min_no_promotion_impl(isless, a, b) +max_no_promotion(a::T, b) where {T} = min_no_promotion_impl(((l, r) -> isless(r, l)), a, b) + +function test_div_significand_with_remainder_iters(::Type{T}) where {T} + range_typemax = function(x) + l = one(T)::T + u = max_nd_for_type(T)::T + a = max_no_promotion(l, u - x) + a:u + end + + range_zero = function(x) + u = max_nd_for_type(T)::T + b = min_no_promotion(u, x) + Base.OneTo(b) + end + + k = 50 + nds = Vector{T}() + append!(nds, range_zero(k)) + append!(nds, range_typemax(k)) + sort!(nds) + unique!(nds) + + msss = Vector{Int}() + max_mss = max_mss_for_type(T) + append!(msss, intersect(Base.OneTo(6), Base.OneTo(max_mss))) + push!(msss, max_mss) + sort!(msss) + unique!(msss) + + test_div_significand_with_remainder_iters(T, nds, msss) +end + +function test_divrem_2(::Type{T}) where {T} + dr = Base.RationalToFloat.divrem_2 + @testset "n: $n" for n ∈ false:T(typemax(Int8)) + @testset "k: $k" for k ∈ 0:6 + @test dr(n, k) == divrem(n, 1 << k, RoundToZero) + end + end +end + +function test_machine_shift(shift::S, ::Type{T}) where {S,T} + @testset "k: $k" for k ∈ 0:(sizeof(T)*8 - 1) + ms = Base.RationalToFloat.machine_shift + o = one(T) + @test ms(shift, o, k) == shift(o, k) + end +end + +const Us = (UInt8, UInt16, UInt32, UInt64, UInt128) +const Ss = (Int8, Int16, Int32, Int64, Int128) +const Ts = (Us..., Ss..., BigInt) + +@testset "machine_shift" begin + @testset "T: $T" for T ∈ Ts + @testset "shift: $shift" for shift ∈ (<<, >>>) + test_machine_shift(shift, T) + end + end +end + +@testset "divrem_2" begin + @testset "T: $T" for T ∈ Ts + test_divrem_2(T) + end +end + +@testset "div_significand_with_remainder" begin + @testset "T: $T" for T ∈ Ts + test_div_significand_with_remainder_iters(T) + end +end