powi: Different rounding behaviour between debug and release mode, and across platforms

[tspiteri]

Note: This was deemed a docs issue; see here for more details.

The code below tries to print the minimum positive subnormal number in different ways.

fn main() {
    println!("{:e}", f64::from_bits(1));
    let divided_minimum_normal = 2f64.powi(-1022) * 2f64.powi(-52);
    println!("{:e}", divided_minimum_normal);
    let direct_pow = 2f64.powi(-1022 - 52);
    println!("{:e}", direct_pow);
}

In release mode, this output is given:

5e-324
5e-324
5e-324

In debug mode, this output is given:

5e-324
5e-324
0e0

That is 2f64.powi(-1074) gives different results for release and debug mode. The difference seems to come from LLVM's optimizations, as the IR looks like it is using llvm.powi.f64(2.0, -1074).

I think this is fine. The difference is only in the subnormal range, and I don't think there are any documented guarantees anywhere that are being broken, but I'm posting this issue to make sure.

(This is on both stable 1.42.0 and on nightly.)

[RalfJung]

Cc @hanna-kruppe @joshtriplett

This is surprising to me. Does IEEE not specify powi behavior? Would be interesting to know what the exact LLVM transformation is that is changing behavior here.

Thanks for bringing this up, at the least this needs to be documented better in the reference (like other floating point subtleties, e.g. #55131).

[RalfJung]

FWIW, Miri prints

5e-324
5e-324
0e0

As expected, it matches the debug behavior. However, Miri implements powi using host floats when it really should use soft-float emulation, so this result does not actually mean that much. apfloat does not seem to support soft-float emulation of powi, though.

[hanna-kruppe]

LLVM performs constant folding of this intrinsic using the host libm, but perhaps it performs a different computations? If I'm reading cppreference right, since C++11 std::pow now longer has a pow(double, int) overload but promotes it to pow(double, double), which might use a very different algorithm than a function specialized for integer exponents. Indeed, if I didn't mess this up then godbolt shows the call used by LLVM turning into pow(base, (double)iexp) in C++11 but calling the same function as Rust's powi in C++98: https://godbolt.org/z/7dqUfM

This might just be an LLVM bug accidentally introduced in the move to C++11.

[RalfJung]

Oh, so LLVM should not use std::pow there, it should use some form of powi?

[hanna-kruppe]

I assume that was the intent of the particular casts it already uses, but I don't know if "some form of powi" is available in C++11 or if it would have to be re-implemented.

[tspiteri]

In {GCC,Clang} {C,C++} you could also write __builtin_powi(f, i).

[tspiteri]

Though at least in GCC, the documentation says that __builtin_powi does not have correct rounding guarantees: “Unlike the pow function no guarantees about precision and rounding are made.”

[tspiteri]

Maybe a comment could be added in the documentation of powi in the same spirit of the comment for log, changing

“Using this function is generally faster than using powf”

to

“Using this function is generally faster than using powf, but owing to implementation details the result may not be correctly rounded.”

[RalfJung]

This looks like an LLVM bug, I don't think we should adjust our docs/spec to that.

[ndmitchell]

I've just wasted a few hours, because I followed the documentation and changed a powf to powi, only to find it failed in some corner case tests later on. I don't think it's a great idea to declare that the docs only have to work with a theoretical compiler backend that doesn't exist - if you are going to reference machine/backend specific characteristics like performance, I think it's only fair to include vastly more important factors like correctness.

I would probably move further and make powi an alias for powf until it got fixed...

[steffahn]

This issue seems concerning to me, and it's existed for a while. In my opinion, unless there exists an LLVM issue for this with any chance of getting fixed in a timely manner, we should either clearly document this problem, or perhaps even implement powi in terms of powf.

[RalfJung]

Has anyone tried bringing this up with the LLVM folks?

OTOH their docs for powi say

This function returns the first value raised to the second power with an unspecified sequence of rounding operations.

For powf:

Return the same value as a corresponding libm ‘pow’ function but without trapping or setting errno.

So possibly this isn't an LLVM bug after all.

I am not sure what would be gained by changing the implementation, but I now agree that documenting that powi and powf could have different rounding behavior makes sense.

[RalfJung]

#73920 has some more examples of powi differences between platforms and build configurations -- not just for subnormals.

I think this is similar to how e.g. sin may return different results depending on how precise the platform approximates the real sine function. Debug vs release mode differences can occur when the optimizer uses a different implementation than would be used at runtime. That all seems completely fine for me. As far as I know, there is no spec that says that powi needs to give a maximally precise answer. (Basic float operations like + and * do give a maximally precise answer, as defined by IEEE 754. But I don't think that applies to powi.)

IOW, I don't think there is a bug here. Maybe we should explicitly document which functions give precise results and which functions give approximations -- and that the quality of the approximation may differ between platforms and build configurations.

[quaternic]

I think this is similar to how e.g. sin may return different results depending on how precise the platform approximates the real sine function. Debug vs release mode differences can occur when the optimizer uses a different implementation than would be used at runtime. That all seems completely fine for me. As far as I know, there is no spec that says that powi needs to give a maximally precise answer. (Basic float operations like + and * do give a maximally precise answer, as defined by IEEE 754. But I don't think that applies to powi.)

The IEEE 754 standard does have a "float to integer power" (pown) under recommended operations, along with "float to float power" (pow), which should both return correctly rounded results. Of course, nothing in Rust documentation claims to provide those, and the current note on powi does imply that there will be rounding error. What it does not point out is how significant that error may be. For example:

let p = 1 << 18;
let x = std::hint::black_box(1.0002812_f32);
let mut x32 = x;
let mut x64 = x as f64;
// compute x.powi(p) with repeated squaring using both f32 and f64
for _ in 0..18 {
    x32 *= x32;
    x64 *= x64;
}
// results are not dependent on platform or optimizations
assert_eq!(x32, 103689867361950008742584042651648.0);
assert_eq!(x64, 102599571306949897080137924476928.0);
assert_eq!(x32 / x64 as f32, 1.0106267);

// that is what powi currently does at runtime
assert_eq!(x32, x.powi(p));
assert_eq!(x64, (x as f64).powi(p));

// with constant folding, powi is accurate instead
assert_eq!(x64 as f32, 1.0002812_f32.powi(p));

The result from f32::powi changes by 1% depending on constant folding, and there are no subnormals involved in this.

[RalfJung]

Of course, nothing in Rust documentation claims to provide those, and the current note on powi does imply that there will be rounding error.

Yeah I think that is the key point here.

I'll close this issue then.

[quaternic]

I still find it quite problematic, what does f32::powi return?

Raises a number to an integer power.

Using this function is generally faster than using powf. It might have a different sequence of rounding operations than powf, so the results are not guaranteed to agree.

fn main() {
    let p = 1 << 18;
    let x = 1.0002812_f32;
    let xp = x.powi(p);
    for dp in -40..40 {
        println!("x^(p{dp:+}) < x^p == {}", x.powi(p+dp) < xp);
    }
}

That program can print:

x^(p-40) < x^p == true
x^(p-39) < x^p == true
x^(p-38) < x^p == true
x^(p-37) < x^p == false
x^(p-36) < x^p == false
...

Sure, the documentation says that powi may not agree with powf, but it's much more surprising that powi can't agree with itself. To better match reality, the documentation would have to say that the possible relative error is nondeterministic and may grow linearly in the exponent.

This is not the same as e.g. f32::sin (#109118) where the differences are in the last few least significant bits, as far as I know.

[RalfJung]

🤷 We can reopen this as a docs issue if you prefer. It would be great if you could submit a PR since you seem to have the required domain knowledge. :)

[quaternic]

Found the corresponding llvm-issue: llvm/llvm-project#65088