Compute a float using an extended-precision representation.
Fast conversion of a the significant digits and decimal exponent a float to an extended representation with a binary float. This algorithm will accurately parse the vast majority of cases, and uses a 128-bit representation (with a fallback 192-bit representation).
This algorithm scales the exponent by the decimal exponent using pre-computed powers-of-5, and calculates if the representation can be unambiguously rounded to the nearest machine float. Near-halfway cases are not handled here, and are represented by a negative, biased binary exponent.
The algorithm is described in detail in "Daniel Lemire, Number Parsing at a Gigabyte per Second" in section 5, "Fast Algorithm", and section 6, "Exact Numbers And Ties", available online: https://arxiv.org/abs/2101.11408.pdf.
pub fn convertEiselLemire(comptime T: type, q: i64, w_: u64) ?BiasedFp(f64)pub fn convertEiselLemire(comptime T: type, q: i64, w_: u64) ?BiasedFp(f64) {
std.debug.assert(T == f16 or T == f32 or T == f64);
var w = w_;
const float_info = FloatInfo.from(T);
// Short-circuit if the value can only be a literal 0 or infinity.
if (w == 0 or q < float_info.smallest_power_of_ten) {
return BiasedFp(f64).zero();
} else if (q > float_info.largest_power_of_ten) {
return BiasedFp(f64).inf(T);
}
// Normalize our significant digits, so the most-significant bit is set.
const lz = @clz(@as(u64, @bitCast(w)));
w = math.shl(u64, w, lz);
const r = computeProductApprox(q, w, float_info.mantissa_explicit_bits + 3);
if (r.lo == 0xffff_ffff_ffff_ffff) {
// If we have failed to approximate w x 5^-q with our 128-bit value.
// Since the addition of 1 could lead to an overflow which could then
// round up over the half-way point, this can lead to improper rounding
// of a float.
//
// However, this can only occur if q ∈ [-27, 55]. The upper bound of q
// is 55 because 5^55 < 2^128, however, this can only happen if 5^q > 2^64,
// since otherwise the product can be represented in 64-bits, producing
// an exact result. For negative exponents, rounding-to-even can
// only occur if 5^-q < 2^64.
//
// For detailed explanations of rounding for negative exponents, see
// <https://arxiv.org/pdf/2101.11408.pdf#section.9.1>. For detailed
// explanations of rounding for positive exponents, see
// <https://arxiv.org/pdf/2101.11408.pdf#section.8>.
const inside_safe_exponent = q >= -27 and q <= 55;
if (!inside_safe_exponent) {
return null;
}
}
const upper_bit = @as(i32, @intCast(r.hi >> 63));
var mantissa = math.shr(u64, r.hi, upper_bit + 64 - @as(i32, @intCast(float_info.mantissa_explicit_bits)) - 3);
var power2 = power(@as(i32, @intCast(q))) + upper_bit - @as(i32, @intCast(lz)) - float_info.minimum_exponent;
if (power2 <= 0) {
if (-power2 + 1 >= 64) {
// Have more than 64 bits below the minimum exponent, must be 0.
return BiasedFp(f64).zero();
}
// Have a subnormal value.
mantissa = math.shr(u64, mantissa, -power2 + 1);
mantissa += mantissa & 1;
mantissa >>= 1;
power2 = @intFromBool(mantissa >= (1 << float_info.mantissa_explicit_bits));
return BiasedFp(f64){ .f = mantissa, .e = power2 };
}
// Need to handle rounding ties. Normally, we need to round up,
// but if we fall right in between and and we have an even basis, we
// need to round down.
//
// This will only occur if:
// 1. The lower 64 bits of the 128-bit representation is 0.
// IE, 5^q fits in single 64-bit word.
// 2. The least-significant bit prior to truncated mantissa is odd.
// 3. All the bits truncated when shifting to mantissa bits + 1 are 0.
//
// Or, we may fall between two floats: we are exactly halfway.
if (r.lo <= 1 and
q >= float_info.min_exponent_round_to_even and
q <= float_info.max_exponent_round_to_even and
mantissa & 3 == 1 and
math.shl(u64, mantissa, (upper_bit + 64 - @as(i32, @intCast(float_info.mantissa_explicit_bits)) - 3)) == r.hi)
{
// Zero the lowest bit, so we don't round up.
mantissa &= ~@as(u64, 1);
}
// Round-to-even, then shift the significant digits into place.
mantissa += mantissa & 1;
mantissa >>= 1;
if (mantissa >= 2 << float_info.mantissa_explicit_bits) {
// Rounding up overflowed, so the carry bit is set. Set the
// mantissa to 1 (only the implicit, hidden bit is set) and
// increase the exponent.
mantissa = 1 << float_info.mantissa_explicit_bits;
power2 += 1;
}
// Zero out the hidden bit
mantissa &= ~(@as(u64, 1) << float_info.mantissa_explicit_bits);
if (power2 >= float_info.infinite_power) {
// Exponent is above largest normal value, must be infinite
return BiasedFp(f64).inf(T);
}
return BiasedFp(f64){ .f = mantissa, .e = power2 };
}