computeProductApprox

fn computeProductApprox(q: i64, w: u64, comptime precision: usize) U128

File

Code

fn computeProductApprox(q: i64, w: u64, comptime precision: usize) U128 {
    std.debug.assert(q >= eisel_lemire_smallest_power_of_five);
    std.debug.assert(q <= eisel_lemire_largest_power_of_five);
    std.debug.assert(precision <= 64);

    const mask = if (precision < 64)
        0xffff_ffff_ffff_ffff >> precision
    else
        0xffff_ffff_ffff_ffff;

    // 5^q < 2^64, then the multiplication always provides an exact value.
    // That means whenever we need to round ties to even, we always have
    // an exact value.
    const index = @as(usize, @intCast(q - @as(i64, @intCast(eisel_lemire_smallest_power_of_five))));
    const pow5 = eisel_lemire_table_powers_of_five_128[index];

    // Only need one multiplication as long as there is 1 zero but
    // in the explicit mantissa bits, +1 for the hidden bit, +1 to
    // determine the rounding direction, +1 for if the computed
    // product has a leading zero.
    var first = U128.mul(w, pow5.lo);
    if (first.hi & mask == mask) {
        // Need to do a second multiplication to get better precision
        // for the lower product. This will always be exact
        // where q is < 55, since 5^55 < 2^128. If this wraps,
        // then we need to need to round up the hi product.
        const second = U128.mul(w, pow5.hi);

        first.lo +%= second.hi;
        if (second.hi > first.lo) {
            first.hi += 1;
        }
    }

    return .{ .lo = first.lo, .hi = first.hi };
}