Maybe a fix for error rounding behaviour of div_euclid

libs
1

TL;DR: 11f32.div_euclid(2.2f32) yields 5 but should be 4. I'm trying fix it.

I had asked deepseek about how to fix the rounding behaviour of div_euclid function.

Here shows what I got:

#![feature(core_float_math)]
use core::f32::math;
use std::time::Instant;
fn div_euclid(a: f32, b: f32) -> f32 {
    let res = (a / b).floor();
    if math::mul_add(b, res, -a) <= 0f32 {
        res
    } else if b > 0f32 {
        res.next_down().floor()
    } else {
        res.next_up().ceil()
    }
}
fn accurate_div_euclid(a: f32, b: f32) -> f32 {
    (a as f64).div_euclid(b as f64) as f32
}
const RANGE: std::ops::Range<i32> = -5000..10001;
fn main() {
    let tester = [
        1f32,
        1.1f32,
        1.21f32,
        1.331f32,
        1.4641f32,
        1.61051f32,
        1.771561f32,
        1.9487171f32,
        -1f32,
        -1.1f32,
        -1.21f32,
        -1.331f32,
        -1.4641f32,
        -1.61051f32,
        -1.771561f32,
        -1.9487171f32,
    ]
    .into_iter()
    .flat_map(|x| {
        RANGE.flat_map(move |y| {
            let base = y as f32 * x;
            if math::mul_add(y as f32, x, -base) >= 0f32 {
                [(base.next_down(), x), (base, x)]
            } else {
                [(base, x), (base.next_up(), x)]
            }
        })
    })
    .collect::<Vec<_>>();
    let now = Instant::now();
    let cum = tester
        .iter()
        .map(|&(x, y)| x.div_euclid(y) as i64)
        .sum::<i64>();
    println!("got {cum}, cost {:?}", now.elapsed());

    let now = Instant::now();
    let cum = tester
        .iter()
        .map(|&(x, y)| div_euclid(x, y) as i64)
        .sum::<i64>();
    println!("got {cum}, cost {:?}", now.elapsed());

    let now = Instant::now();
    let cum = tester
        .iter()
        .map(|&(x, y)| accurate_div_euclid(x, y) as i64)
        .sum::<i64>();
    println!("got {cum}, cost {:?}", now.elapsed());

    for (x, y) in tester {
        if div_euclid(x, y) as i64 != accurate_div_euclid(x, y) as i64 {
            println!(
                "{x}/{y}: left = {}, right = {}",
                div_euclid(x, y) as i64,
                accurate_div_euclid(x, y) as i64
            )
        }
    }
}

I tried tested this code and surprisingly found that convert to f64 yields faster result than call div_euclid directly.

Thus I have little confident whether my code is better than the original one.

2

2.2f32 does not exist as a number – 2.2 is not one of the values an f32 (or f64 for that matter) can represent. As such, when you write 2.2f32, the compiler will need to pick a nearby number that can be represented, such as 2.2000000476837158203125. 11 divided by 2.2000000476837158203125 is slightly less than 5, so Rust correctly rounds the result down to 4.

The only way I can think of to address this sort of problem within the compiler would be to give a warning when expressing a floating-point number as a constant that can't be exactly converted to a f32/f64 as appropriate, to let people know that they are trying to do something that (fixed-precision) floating-point numbers can't do. But most of the time, when someone writes 2.2f32 into a program, they are looking for a number approximately equal to 2.2 rather than the exact value, so such a warning would have a huge number of false positives.

As such, the best solution is probably just to avoid using floating point in situations where you need exact results: in practice, most floating-point calculations end up going via numbers that can't be exactly represented using the floating point type you're doing the calculation in, and thus end up producing slightly inaccurate or approximated results due to rounding.

3

The problem is that, here, 11f32 % 2.2f32 ~2.1999f32, not zero. This inconsistance make the div_euclid unusable in some case.

I first observed this problem in a game called Elin, when proform some action, 11 inputs will turn to 10 outputs, and the remain part have a small probability convert to a new item.

The author has to use input / 1.1 + prob((input % 1.1)/1.1) to calculate the total output count.

There might be no other simpler solution here.

Since rust has documented div_euclid works for infinitely precision, it might be better to fix the div_euclid algorithm.

4

div_euclid works at infinite precision. However, converting the text "2.2" into an f32 is not infinite precision.

If you want exact values for multiplication/division/addition/subtraction, then I'd suggest using the num_rational crate, which will have some performance cost.

5

Here, precision is not the concern. For a game (or something similar), divided by 1.1f32.next_up() or 1.1f32.next_done() are nearly identical things and should not yield significant difference.

The problem is that, a.div_euclid(b) * b +a.rem_euclid(b) is not (even approximately) equals to a.

This is why I submit such fix.

Maybe we should notice that, for many divisors, trivial division is different from the div_euclid method

6

Getting very different answers from very similar inputs is expected, and is the correct behavior.

It's like how 10.0.div_euclid(1.0) should be 10, but 10.0.div_euclid(1.0001) should be 9.