|
| 1 | +# This file is a part of Julia. License is MIT: https://julialang.org/license |
| 2 | + |
| 3 | +module RationalToFloat |
| 4 | + |
| 5 | +# Performance optimization. Unlike raw `<<` or `>>>`, this is supposed |
| 6 | +# to compile to a single instruction, because the semantics correspond |
| 7 | +# to what hardware usually provides. |
| 8 | +function machine_shift(shift::S, a::T, b) where {S,T<:Base.BitInteger} |
| 9 | + @inline begin |
| 10 | + mask = 8*sizeof(T) - 1 |
| 11 | + c = b & mask |
| 12 | + shift(a, c) |
| 13 | + end |
| 14 | +end |
| 15 | + |
| 16 | +machine_shift(::S, a::Bool, ::Any) where {S} = error("unsupported") |
| 17 | + |
| 18 | +# Fallback for `BigInt` etc. |
| 19 | +machine_shift(shift::S, a, b) where {S} = shift(a, b) |
| 20 | + |
| 21 | +# Arguments are positive integers. |
| 22 | +function div_significand_with_remainder(num, den, minimum_significand_size) |
| 23 | + clamped = x -> max(zero(x), x)::typeof(x) |
| 24 | + bw = Base.top_set_bit # bit width |
| 25 | + shift = clamped(minimum_significand_size + bw(den) - bw(num) + 0x2) |
| 26 | + t = machine_shift(<<, num, shift) |
| 27 | + (divrem(t, den, RoundToZero)..., shift) |
| 28 | +end |
| 29 | + |
| 30 | +# `divrem(n, 1<<k, RoundToZero)` |
| 31 | +function divrem_2(n, k) |
| 32 | + quo = machine_shift(>>>, n, k) |
| 33 | + tmp = machine_shift(<<, quo, k) |
| 34 | + rem = n - tmp |
| 35 | + (quo, rem) |
| 36 | +end |
| 37 | + |
| 38 | +function to_float_components_impl(num, den, precision, max_subnormal_exp) |
| 39 | + # `+1` because we need an extra, "round", bit for some rounding modes. |
| 40 | + # |
| 41 | + # TODO: as a performance optimization, only do this when required |
| 42 | + # by the rounding mode |
| 43 | + prec_p_1 = precision + true |
| 44 | + |
| 45 | + (quo0, rem0, shift) = div_significand_with_remainder(num, den, prec_p_1) |
| 46 | + width = Base.top_set_bit(quo0) |
| 47 | + excess_width = width - prec_p_1 |
| 48 | + exp = width - shift - true |
| 49 | + |
| 50 | + exp_underflow = if isnothing(max_subnormal_exp) |
| 51 | + zero(exp) |
| 52 | + else |
| 53 | + let d = max_subnormal_exp - exp, T = typeof(d), z = zero(d)::T |
| 54 | + (signbit(d) ? z : d + true)::T |
| 55 | + end |
| 56 | + end |
| 57 | + |
| 58 | + (quo1, rem1) = divrem_2(quo0, exp_underflow + excess_width) |
| 59 | + integral_significand = quo1 >>> true |
| 60 | + round_bit = quo1 % Bool |
| 61 | + sticky_bit = !iszero(rem1) | !iszero(rem0) |
| 62 | + |
| 63 | + (; integral_significand, exponent = exp, round_bit, sticky_bit) |
| 64 | +end |
| 65 | + |
| 66 | +struct RoundingIncrementHelper |
| 67 | + final_bit::Bool |
| 68 | + round_bit::Bool |
| 69 | + sticky_bit::Bool |
| 70 | +end |
| 71 | + |
| 72 | +(h::RoundingIncrementHelper)(::Base.Rounding.FinalBit) = h.final_bit |
| 73 | +(h::RoundingIncrementHelper)(::Base.Rounding.RoundBit) = h.round_bit |
| 74 | +(h::RoundingIncrementHelper)(::Base.Rounding.StickyBit) = h.sticky_bit |
| 75 | + |
| 76 | +function to_float_components_rounded(num, den, precision, max_subnormal_exp, romo, sign_bit) |
| 77 | + overflows = (x, p) -> x == machine_shift(<<, one(x), p) |
| 78 | + t = to_float_components_impl(num, den, precision, max_subnormal_exp) |
| 79 | + raw_significand = t.integral_significand |
| 80 | + rh = RoundingIncrementHelper(raw_significand % Bool, t.round_bit, t.sticky_bit) |
| 81 | + incr = Base.Rounding.correct_rounding_requires_increment(rh, romo, sign_bit) |
| 82 | + rounded = raw_significand + incr |
| 83 | + (integral_significand, exponent) = let exp = t.exponent |
| 84 | + if overflows(rounded, precision) |
| 85 | + (rounded >>> true, exp + true) |
| 86 | + else |
| 87 | + (rounded, exp) |
| 88 | + end |
| 89 | + end |
| 90 | + (; integral_significand, exponent) |
| 91 | +end |
| 92 | + |
| 93 | +function to_float_components(::Type{T}, num, den, precision, max_subnormal_exp, romo, sb) where {T} |
| 94 | + to_float_components_rounded(abs(T(num)), den, precision, max_subnormal_exp, romo, sb) |
| 95 | +end |
| 96 | + |
| 97 | +function to_floating_point_fallback(::Type{T}, ::Type{S}, num, den, rm, prec) where {T,S} |
| 98 | + num_is_zero = iszero(num) |
| 99 | + den_is_zero = iszero(den) |
| 100 | + sb = signbit(num) |
| 101 | + is_zero = num_is_zero & !den_is_zero |
| 102 | + is_inf = !num_is_zero & den_is_zero |
| 103 | + is_regular = !num_is_zero & !den_is_zero |
| 104 | + if is_regular |
| 105 | + let |
| 106 | + c = to_float_components(S, num, den, prec, nothing, rm, sb) |
| 107 | + signif = T(c.integral_significand)::T |
| 108 | + x = ldexp(signif, c.exponent - prec + true)::T |
| 109 | + sb ? -x : x |
| 110 | + end::T |
| 111 | + else |
| 112 | + if is_zero |
| 113 | + zero(T)::T |
| 114 | + elseif is_inf |
| 115 | + T(Inf)::T |
| 116 | + else |
| 117 | + T(NaN)::T |
| 118 | + end |
| 119 | + end::T |
| 120 | +end |
| 121 | + |
| 122 | +function to_floating_point_impl(::Type{T}, ::Type{S}, num, den, rm, prec) where {T,S} |
| 123 | + to_floating_point_fallback(T, S, num, den, rm, prec) |
| 124 | +end |
| 125 | + |
| 126 | +function to_floating_point_impl(::Type{T}, ::Type{S}, num, den, rm, prec) where {T<:Base.IEEEFloat,S} |
| 127 | + num_is_zero = iszero(num) |
| 128 | + den_is_zero = iszero(den) |
| 129 | + sb = signbit(num) |
| 130 | + is_zero = num_is_zero & !den_is_zero |
| 131 | + is_inf = !num_is_zero & den_is_zero |
| 132 | + is_regular = !num_is_zero & !den_is_zero |
| 133 | + Rnd = Base.Rounding |
| 134 | + (rm_is_to_zero, rm_is_from_zero) = if Rnd.rounds_to_nearest(rm) |
| 135 | + (false, false) |
| 136 | + else |
| 137 | + let from = Rnd.rounds_away_from_zero(rm, sb) |
| 138 | + (!from, from) |
| 139 | + end |
| 140 | + end::NTuple{2,Bool} |
| 141 | + exp_max = Base.exponent_max(T) |
| 142 | + exp_min = Base.exponent_min(T) |
| 143 | + ieee_repr = Base.ieee754_representation |
| 144 | + repr_zero = ieee_repr(T, sb, Val(:zero)) |
| 145 | + repr_inf = ieee_repr(T, sb, Val(:inf)) |
| 146 | + repr_nan = ieee_repr(T, sb, Val(:nan)) |
| 147 | + U = typeof(repr_zero) |
| 148 | + repr_zero::U |
| 149 | + repr_inf::U |
| 150 | + repr_nan::U |
| 151 | + |
| 152 | + ret_u = if is_regular |
| 153 | + let |
| 154 | + c = to_float_components(S, num, den, prec, exp_min - 1, rm, sb) |
| 155 | + exp = c.exponent |
| 156 | + exp_diff = exp - exp_min |
| 157 | + is_normal = 0 ≤ exp_diff |
| 158 | + exp_is_huge_p = exp_max < exp |
| 159 | + exp_is_huge_n = signbit(exp_diff + prec) |
| 160 | + rounds_to_inf = exp_is_huge_p & !rm_is_to_zero |
| 161 | + rounds_to_zero = exp_is_huge_n & !rm_is_from_zero |
| 162 | + |
| 163 | + if !rounds_to_zero & !exp_is_huge_p |
| 164 | + let signif = (c.integral_significand % U) & Base.significand_mask(T) |
| 165 | + exp_field = (max(exp_diff, zero(exp_diff)) + is_normal) % U |
| 166 | + ieee_repr(T, sb, exp_field, signif)::U |
| 167 | + end |
| 168 | + elseif rounds_to_zero |
| 169 | + repr_zero |
| 170 | + elseif rounds_to_inf |
| 171 | + repr_inf |
| 172 | + else |
| 173 | + ieee_repr(T, sb, Val(:omega)) |
| 174 | + end |
| 175 | + end |
| 176 | + else |
| 177 | + if is_zero |
| 178 | + repr_zero |
| 179 | + elseif is_inf |
| 180 | + repr_inf |
| 181 | + else |
| 182 | + repr_nan |
| 183 | + end |
| 184 | + end::U |
| 185 | + |
| 186 | + reinterpret(T, ret_u)::T |
| 187 | +end |
| 188 | + |
| 189 | +# `BigInt` is a safe default. |
| 190 | +to_float_promote_type(::Type{F}, ::Type{S}) where {F,S} = BigInt |
| 191 | + |
| 192 | +const BitIntegerOrBool = Union{Bool,Base.BitInteger} |
| 193 | + |
| 194 | +# As an optimization, use an integer type narrower than `BigInt` when possible. |
| 195 | +function to_float_promote_type(::Type{F}, ::Type{S}) where {F<:Base.IEEEFloat,S<:BitIntegerOrBool} |
| 196 | + Max = if sizeof(F) ≤ sizeof(S) |
| 197 | + S |
| 198 | + else |
| 199 | + (S <: Signed) ? Base.inttype(F) : Base.uinttype(F) |
| 200 | + end |
| 201 | + widen(Max) |
| 202 | +end |
| 203 | + |
| 204 | +function to_floating_point(::Type{F}, num::T, den::T, rm, prec) where {F,T} |
| 205 | + S = to_float_promote_type(F, T) |
| 206 | + to_floating_point_impl(F, S, num, den, rm, prec) |
| 207 | +end |
| 208 | + |
| 209 | +end |
0 commit comments