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. (EDIT: The correct answer is 4, but apparently Rust is producing 5 instead, so my statement that the compiler is correct was wrong.)

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.

9 Likes
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.

4 Likes
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.

5 Likes
7

Though 10.0_f32.rem_euclid(1.0001) returns something which is almost 1.

The concern is that 11f32.div_euclid(2.2f32) and 11f32.rem_euclid(2.2f32) disagree.

Still, floats have different precisions at different scales, so I strongly suspect that quirks like this are inevitable.

8

How about using a.rem_euclid(b), and (a - a.rem_euclid(b)) / b instead, if f(a, b) * b + a.rem_euclid(b) ≈ a is the important part?

The std docs claiming that a ≈ a.div_euclid(b) * b + a.rem_euclid(b) is unfortunate, since I think that approximation and the property about infinite precision might not both be able to be satisfied.

9

Ah, this is a known issue. f64::div_euclid and f64::rem_euclid yield inconsistent results · Issue #107904 · rust-lang/rust · GitHub

Someone just needs to put in the work to revive the fix PR at https://github.com/rust-lang/rust/pull/134145

6 Likes
10

The proper solution is to tell the CPU to use the round-down rounding mode instead of the round-to-nearest rounding mode for the division. The problem is that the Rust developers decided that using CPU rounding modes other than round-to-nearest would be instant UB.

If you have an AVX-512 CPU, then it can be implemented using _mm_div_round_ss. Otherwise, you need to use inline assembly, since the function to do this on other x86 CPUs specifies that using it to fix this "leads to immediate Undefined Behavior".

If you're not on x86, then you'll probably need to use inline assembly.

11

The phrasing of the documentation suggests to me that this is because of limitations in LLVM, probably boiling down to "nobody has ever bothered to implement full support for #pragma STDC FENV_ACCESS ON in LLVM."

(My initial guess was concurrency issues, but C2011 (which is the oldest edition of the C standard that has a concept of threads) clearly requires the floating-point environment to be thread-specific.)

Your best bet for getting this fixed is to find someone with deep pockets who cares as much about floating point correctness as William Kahan did, and get them to fund the LLVM-level fixes. Or -- this may actually be easier -- convince Intel and ARM to add per-instruction rounding direction controls to every FP instruction (RISC-V already has this IIRC), which would enable round-to-±infinity arithmetic to be supported via methods, e.g. x.add_round_up(y). Once everybody has the latest hardware.

1 Like
12

We (t-opsem) do reserve the right to remove any UB, so it's possible that we would remove this UB were it to become practical to support. But, even ignoring whether LLVM's fp-strict mode is sufficient, there are some fundamental problems that would need to be addressed:

  • Assuming the default fpenv is required to do such basic things as fold and reorder floating point operations. With true strict fpenv, floating point operations basically need to be treated as little black-box operations we know nothing about beyond the processor-defined call ABI.
  • As such, most functions are compiled assuming the default fpenv. If you are in strict fpenv mode you must not call a function that isn't in strict fpenv mode. This includes the standard library.
    • Furthermore, some algorithms might have been written with the assumption of the default fpenv and would be incorrect under a different fpenv. So it's not possible to just recompile things under strict fpenv if necessary.
  • The Rust Abstract Machine as currently proposed does not have a concept of the fpenv as a mutable bit of execution state. This would be the easiest to fix, but doing so is the source of all the other complications.
  • Enabling strict fpenv would need to be handled like conditional target feature is handled today (i.e. prevent inlining into code that might not require it be enabled) but also handle the idea that strict fpenv can be turned off again, which isn't the case for target feature support (once it's seen the entire program is known to have it available; that's what allows for caching feature tests).
  • Most everything else is "just" implementation limitations rather than a fundamental complication of strict fpenv support.
    • We recognize the desire to and benefit of being able to write things like this in Rust syntax using intrinsics instead of inline assembly. But it's complicated, and Rust fundamentally can't be used as a "fancy macro assembler" like some of the more exotic vendor intrinsics end up implicitly assuming C can (i.e. when they're defined as emitting some specific machine code rather than behavior-derived "as-if" storytelling).
13

What's stopping us from implementing our own library solution to this using inline assembly? I know typically we only like to add intrinsics that have a direct representation in LLVM. But couldn't we add a fn div_round_up (or div_with_rounding) with inline assembly implementations on tier 1&2 targets, with slow backup code for other cases?

14

This is kinda tangential, but: in at least one of Kahan's "it sucks that nobody wants to implement the full scope of my original vision for IEEE 754" essays, he talks about how the point of having the rounding mode be dynamically scoped (in the Lisp sense) is, you're supposed to be able to run your entire program three times without recompiling it all, once each in round-to-nearest, round-to-negative-infinity, and round-to-positive-infinity mode, and if the results are more than ε different from each other you have a bug.

I am not convinced this actually makes sense as a debugging tactic, but it clearly does mean Kahan thinks a compiler shouldn't need to care what the rounding mode is, even if it is aggressively optimizing of floating-point code. And he also evidently thinks that all, or almost all, individual functions can and should be written such that they give correct results regardless of rounding mode.

1 Like