k = ((min_exp + DiyFP::SIGNIFICAND_SIZE - 1) * D_1_LOG2_10).ceil

index = ((CACHED_POWER_OFFSET + k.to_i - 1) / CACHED_EXP_STEP) + 1

pow = PowCache[index]

return DiyFP.new(pow.significand, pow.binary_exp), pow.decimal_exp.to_i

end

#

# This result is not normalized.

def -(other : DiyFP)

self.class.new(frac - other.frac, exp)

end

end

def normalize

f = frac

e = exp

end

def self.from_f(d : Float64 | Float32)

frac, exp = IEEE.frac_and_exp(d)

new(frac, exp)

end

e -= DiyFP::SIGNIFICAND_SIZE - IEEE::SIGNIFICAND_SIZE_64

DiyFP.new(f, e)

end

# Conceptually rest ~= too_high - buffer

# We need to do the following tests in this order to avoid over- and

# underflows.

while (

rest < small_distance && # Negated condition 1

unsafe_interval - rest >= ten_kappa && # Negated condition 2

# imprecision.

def digit_gen(low : DiyFP, w : DiyFP, high : DiyFP, buffer_p) : {Bool, Int32, Int32}

buffer = buffer_p.to_slice(128)

# low, w and high are imprecise, but by less than one ulp (unit in the last

# place).

# If we remove (resp. add) 1 ulp from low (resp. high) we are certain that

while kappa > 0

digit = integrals / divisor

# pp [digit, kappa]

buffer[length] = 48_u8 + digit

length += 1

integrals %= divisor

# data (like the interval or 'unit'), too.

# Note that the multiplication by 10 does not overflow, because w.e >= -60

# and thus one.e >= -60.

loop do

fractionals *= 10

unit *= 10

unsafe_interval = DiyFP.new(unsafe_interval.frac * 10, unsafe_interval.exp)

digit = (fractionals >> -one.exp).to_i

buffer[length] = 48_u8 + digit

length += 1

fractionals &= one.frac - 1

# boundary_minus and boundary_plus will round to v when convert to a float.

# Grisu3 will never output representations that lie exactly on a boundary.

boundaries = IEEE.normalized_boundaries(v)

ten_mk, mk = CachedPowers.get_cached_power_for_binary_exponent(w.exp)

# In other words: let f = scaled_w.f() and e = scaled_w.e(), then

# (f-1) * 2^e < w*10^k < (f+1) * 2^e

scaled_w = w * ten_mk

# In theory it would be possible to avoid some recomputations by computing

# the difference between w and boundary_minus/plus (a power of 2) and to

decimal_exponent = -mk + kappa

return result, decimal_exponent, length

end

# exponent as m_plus.

# Precondition: the value encoded by this Flaot must be greater than 0.

def normalized_boundaries(v : Float64)

w = DiyFP.from_f(v)

m_plus = DiyFP.new((w.frac << 1) + 1, w.exp - 1).normalize

end

def normalized_boundaries(v : Float32)

w = DiyFP.from_f(v)

m_plus = DiyFP.new((w.frac << 1) + 1, w.exp - 1).normalize

def frac_and_exp(v : Float64)

d64 = to_uint(v)

if (d64 & EXPONENT_MASK_64) == 0 # denormal float

frac = d64 & SIGNIFICAND_MASK_64

def frac_and_exp(v : Float32)

d32 = to_uint(v)

if (d32 & EXPONENT_MASK_32) == 0 # denormal float

frac = d32 & SIGNIFICAND_MASK_32

baised_e = ((d32 & EXPONENT_MASK_32) >> PHYSICAL_SIGNIFICAND_SIZE_32).to_i

baised_e - EXPONENT_BIAS_32

end