Compute the exact inverse of this "simple" floating-point function

Question

I have the following function:

float int_to_qty(unsigned x) {
    const float MAX = 8.5f;
    const float MIN = .001f;
    return ((MAX-MIN) / (float)(1<<24)) * x + MIN;
}

This compiles (with reasonable options, on x86) to the following:

.LCPI0_0:
        .long   0x3507fbe7                      # float 5.06579852E-7
.LCPI0_1:
        .long   0x3a83126f                      # float 0.00100000005
int_to_qty:                             # @int_to_qty
        mov     eax, edi
        cvtsi2ss        xmm0, rax
        mulss   xmm0, dword ptr [rip + .LCPI0_0]
        addss   xmm0, dword ptr [rip + .LCPI0_1]
        ret

I consider the assembly to be the "canonical" version of the function: Convert the int to a float, multiply by a constant at 32-bit precision, add another constant at 32-bit precision, that's the result.

I want to find the exact inverse of this function. Specifically, a function unsigned qty_to_int(float qty) that will pass the following test:

int test() {
    for (unsigned i = 0; i < (1 << 24); ++i) {
        float qty = int_to_qty(i);
        if (int_to_qty(qty_to_int(qty)) != qty) {
            return 0;
        }
    }
    return 1;
}

Notes:

• In the range 4 ≤ int_to_qty(x) < 8, the returned values primarily differ by 1 ulp, which is what makes this challenging.
• In the range 8 ≤ int_to_qty(x) < 8.5, the function stops being one-to-one. In this case either answer is fine for the inverse, it doesn't have to be consistently the lowest or the highest.

Answer 1

After wrestling for a long time, I finally came up with a solution that passes the tests. (In Rust, but the translation to C is straightforward.)

pub fn qty_to_int(qty: f64) -> u32 {
    const MAX: f32 = 8.5;
    const MIN: f32 = 0.001;
    let size_inv = f64::from(1<<24) / f64::from(MAX - MIN);

    // We explicitly shrink the precision to f32 and then pop back to f64, because we *need* to
    // perform that rounding step to properly reverse the addition at the end of int_to_qty.
    // We could do the whole thing at f32 precision, except that our input is f64 so the
    // subtraction needs to be done at f64 precision.
    let fsqueezed: f32 = (qty - f64::from(MIN)) as f32;
    // The squeezed subtraction is a one-to-one operation across most of our range. *However*,
    // in the border areas where our input is just above an exponent breakpoint, but
    // subtraction will bring it below, we hit an issue: The addition in int_to_qty() has
    // irreversibly lost a bit in the lowest place! This causes issues when we invert the
    // multiply, since we are counting on the error to be centered when we round at the end.
    //
    // To solve this, we need to re-bias the error by subtracting 0.25 ulp for these cases.
    // Technically this should be applied for ranges like [2.0,2.001) as well, but they
    // don't need it since they already round correctly (due to having more headroom).
    let adj = if qty > 4.0 && qty < 4.001 {
        0.5 - 2.0 / 8.499
    } else if qty > 8.0 && qty < 8.001 {
        0.5 - 4.0 / 8.499
    } else {
        0.5
    };
    // Multiply and round, taking into account the possible adjustments.
    let fresult = f64::from(fsqueezed) * size_inv + adj;
    unsafe { fresult.to_int_unchecked() }
}