1const std = @import("std");
  2const math = std.math;
  3const common = @import("common.zig");
  4const BiasedFp = common.BiasedFp;
  5const Decimal = @import("decimal.zig").Decimal;
  6const mantissaType = common.mantissaType;
  7
  8const max_shift = 60;
  9const num_powers = 19;
 10const powers = [_]u8{ 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, 33, 36, 39, 43, 46, 49, 53, 56, 59 };
 11
 12pub fn getShift(n: usize) usize {
 13    return if (n < num_powers) powers[n] else max_shift;
 14}
 15
 16/// Parse the significant digits and biased, binary exponent of a float.
 17///
 18/// This is a fallback algorithm that uses a big-integer representation
 19/// of the float, and therefore is considerably slower than faster
 20/// approximations. However, it will always determine how to round
 21/// the significant digits to the nearest machine float, allowing
 22/// use to handle near half-way cases.
 23///
 24/// Near half-way cases are halfway between two consecutive machine floats.
 25/// For example, the float `16777217.0` has a bitwise representation of
 26/// `100000000000000000000000 1`. Rounding to a single-precision float,
 27/// the trailing `1` is truncated. Using round-nearest, tie-even, any
 28/// value above `16777217.0` must be rounded up to `16777218.0`, while
 29/// any value before or equal to `16777217.0` must be rounded down
 30/// to `16777216.0`. These near-halfway conversions therefore may require
 31/// a large number of digits to unambiguously determine how to round.
 32///
 33/// The algorithms described here are based on "Processing Long Numbers Quickly",
 34/// available here: <https://arxiv.org/pdf/2101.11408.pdf#section.11>.
 35///
 36/// Note that this function needs a lot of stack space and is marked
 37/// cold to hint against inlining into the caller.
 38pub fn convertSlow(comptime T: type, s: []const u8) BiasedFp(T) {
 39    @branchHint(.cold);
 40
 41    const MantissaT = mantissaType(T);
 42    const min_exponent = -(1 << (math.floatExponentBits(T) - 1)) + 1;
 43    const infinite_power = (1 << math.floatExponentBits(T)) - 1;
 44    const fractional_bits = math.floatFractionalBits(T);
 45
 46    var d = Decimal(T).parse(s); // no need to recheck underscores
 47    if (d.num_digits == 0 or d.decimal_point < Decimal(T).min_exponent) {
 48        return BiasedFp(T).zero();
 49    } else if (d.decimal_point >= Decimal(T).max_exponent) {
 50        return BiasedFp(T).inf(T);
 51    }
 52
 53    var exp2: i32 = 0;
 54    // Shift right toward (1/2 .. 1]
 55    while (d.decimal_point > 0) {
 56        const n = @as(usize, @intCast(d.decimal_point));
 57        const shift = getShift(n);
 58        d.rightShift(shift);
 59        if (d.decimal_point < -Decimal(T).decimal_point_range) {
 60            return BiasedFp(T).zero();
 61        }
 62        exp2 += @as(i32, @intCast(shift));
 63    }
 64    //  Shift left toward (1/2 .. 1]
 65    while (d.decimal_point <= 0) {
 66        const shift = blk: {
 67            if (d.decimal_point == 0) {
 68                break :blk switch (d.digits[0]) {
 69                    5...9 => break,
 70                    0, 1 => @as(usize, 2),
 71                    else => 1,
 72                };
 73            } else {
 74                const n = @as(usize, @intCast(-d.decimal_point));
 75                break :blk getShift(n);
 76            }
 77        };
 78        d.leftShift(shift);
 79        if (d.decimal_point > Decimal(T).decimal_point_range) {
 80            return BiasedFp(T).inf(T);
 81        }
 82        exp2 -= @as(i32, @intCast(shift));
 83    }
 84    // We are now in the range [1/2 .. 1] but the binary format uses [1 .. 2]
 85    exp2 -= 1;
 86    while (min_exponent + 1 > exp2) {
 87        var n = @as(usize, @intCast((min_exponent + 1) - exp2));
 88        if (n > max_shift) {
 89            n = max_shift;
 90        }
 91        d.rightShift(n);
 92        exp2 += @as(i32, @intCast(n));
 93    }
 94    if (exp2 - min_exponent >= infinite_power) {
 95        return BiasedFp(T).inf(T);
 96    }
 97
 98    // Shift the decimal to the hidden bit, and then round the value
 99    // to get the high mantissa+1 bits.
100    d.leftShift(fractional_bits + 1);
101    var mantissa = d.round();
102    if (mantissa >= (@as(MantissaT, 1) << (fractional_bits + 1))) {
103        // Rounding up overflowed to the carry bit, need to
104        // shift back to the hidden bit.
105        d.rightShift(1);
106        exp2 += 1;
107        mantissa = d.round();
108        if ((exp2 - min_exponent) >= infinite_power) {
109            return BiasedFp(T).inf(T);
110        }
111    }
112    var power2 = exp2 - min_exponent;
113    if (mantissa < (@as(MantissaT, 1) << fractional_bits)) {
114        power2 -= 1;
115    }
116    // Zero out all the bits above the mantissa bits.
117    mantissa &= (@as(MantissaT, 1) << math.floatMantissaBits(T)) - 1;
118    return .{ .f = mantissa, .e = power2 };
119}