to_degrees/to_radians aren't rounded correctly

[huonw]

#![feature(float_extras)]

fn main() {
    let exact = 57.2957795130823208767981548141051703324054724665643215491602438612028471483215526324409689958511109441862233816328648932;

    assert_eq!(1_f32.to_degrees(), exact);
}

thread '<main>' panicked at 'assertion failed: `(left == right)` (left: `57.295776`, right: `57.29578`)', <anon>:6
playpen: application terminated with error code 101

That computed answer is just one bit wrong.

The above case would be addressed by having the 180/pi value used in the conversion being a correctly rounded constant, rather than computed as 180.0 / consts::PI, although I'm not sure this will get the right answer in all cases. If we decide to care about this, then we'll need an exact multiplication algorithm such as Brisebarre and Muller (2008) Correctly rounded multiplication by arbitrary precision constants. Let C be the exact (i.e. infinite precision) value of 180 / pi, then:

// C as f32
const C_H: f32 = 57.295780181884765625;
// (C - C_H) as f32
const C_L: f32 = -6.6880244276035227812826633453369140625e-7;

fn exact_mul(x: f32) -> f32 {
    let a = C_L * x;
    let b = C_H.mul_add(x, a);
    b
}

fn main() {
    let exact = 57.2957795130823208767981548141051703324054724665643215491602438612028471483215526324409689958511109441862233816328648932;
    assert_eq!(exact_mul(1.0), exact);
}

(There's some conditions on the constant C for this to work for all x, and I don't know/haven't yet checked if 180/pi satisfies them.)

This will is noticably slower than the naive method, especially if there's not hardware support for FMA (mul_add), however, if we do implement this, people who don't care about accuracy can use the naive method trivially. That said, to_degrees is almost always used for output for humans, and is often rounded to many fewer decimal places than full precision, so loss of precision doesn't matter... but also, speed probably doesn't matter so much (the formatting will almost certainly be much slower than the multiplication). On the other hand, to_radians won't be used for human output, typically.

#![feature(test, float_extras)]
extern crate test;

// C as f32
const C_H: f32 = 57.295780181884765625;
// (C - C_H) as f32
const C_L: f32 = -6.6880244276035227812826633453369140625e-7;

fn exact_mul(x: f32) -> f32 {
    let a = C_L * x;
    let b = C_H.mul_add(x, a);
    b
}

const FLOATS: &'static [f32] = &[-360.0, -359.99, -100.0, -1.001, 0.001, -1e-40,
                                 0.0,
                                 1e-40, 0.001, 1.001, 100.0, 359.99, 360.0];

#[bench]
fn exact(b: &mut test::Bencher) {
    b.iter(|| {
        for x in test::black_box(FLOATS) {
            test::black_box(exact_mul(*x));
        }
    })
}

#[bench]
fn inexact(b: &mut test::Bencher) {
    b.iter(|| {
        for x in test::black_box(FLOATS) {
            test::black_box(x.to_degrees());
        }
    })
}

#[bench]
fn format(b: &mut test::Bencher) {
    b.iter(|| {
        let mut buf = [0u8; 100];
        for x in test::black_box(FLOATS) {
            use std::io::prelude::*;
            let _ = write!(&mut buf as &mut [_], "{:.0}", x);
            test::black_box(&buf);
        }
    })
}

test exact   ... bench:         169 ns/iter (+/- 6)
test format  ... bench:       1,238 ns/iter (+/- 41)
test inexact ... bench:          84 ns/iter (+/- 3)

(It's likely that to_radians suffers similarly, since pi / 180 is just as transcendental as 180 / pi, but I haven't found an example with a few tests.)

---

This issue is somewhat of a policy issue: how much do we care about getting correctly rounded floating point answers? This case is pretty simple, but we almost certainly don't want to guarantee 100% precise answers for functions like sin/exp etc., since this is non-trivial/slow to do, and IEEE754 (and typical libm's) don't guarantee exact rounding for all inputs.

[Gankra]

This seems to not be super important based on:

That said, to_degrees is almost always used for output for humans, and is often rounded to many fewer decimal places than full precision, so loss of precision doesn't matter

Unless someone finds some really degenerate precision loss.

While I agree that the perf hit isn't a big deal, I don't see a reason to take the hit unless someone's really hurting for the precision (and then they can of course roll their own in the worst-case).

[hanna-kruppe]

Although the vast majority of code I have contributed to Rust so far is dedicated to correct rounding of floats, I agree. If it's always within 1 ULP, then it's most likely perfectly okay. Theoretically some badly-written code might round-trip between degrees and radians n times, but that is a silly thing to do and dangerous anyway: You're bound to mess up with the units, unless you use a newtype wrapper, but a newtype wrapper can easily be augmented with a custom, correctly-rounded implementation of the conversions.

[huonw]

If it's always within 1 ULP, then it's most likely perfectly okay

Ok, I did an exhaustive test of all f32, and the maximum error (for values for which there's no overflow/underflow) seems to be 2 ulp:

#![feature(float_extras)]
extern crate ieee754; // [dependencies] ieee754 = "0.2"
use ieee754::Ieee754;

use std::f32;

// C as f32
const C_H: f32 = 57.295780181884765625;
// (C - C_H) as f32
const C_L: f32 = -6.6880244276035227812826633453369140625e-7;

fn exact_mul(x: f32) -> f32 {
    let a = C_L * x;
    let b = C_H.mul_add(x, a);
    b
}

fn main() {
    let mut errs = vec![];
    let mut count = 0_usize;
    for x in f32::MIN.upto(f32::MAX) {
        count += 1;
        let exact = exact_mul(x);
        let approx = x.to_degrees();
        if exact != approx {
            if let Some(ulp) = exact.ulp() {
                let ulp_error = (exact - approx).abs() / ulp;
                let bits_wrong = ulp_error.round() as u32;
                errs.push(bits_wrong);
            }
        }
    }
    println!("{}/{} wrong, max ulp err {}",
             errs.len(), count,
             errs.iter().cloned().max().unwrap_or(0))
}

2798869104/4278190079 wrong, max ulp err 2

(Warning, this program requires ~11GB of RAM to finish.)

(There's some conditions on the constant C for this to work for all x, and I don't know/haven't yet checked if 180/pi satisfies them.)

(I've now checked, and it does in single precision, and for all doubles except 8357367214274750 * 2ⁿ.)

[hanna-kruppe]

Thank you for testing this! 1 ULP vs 2 ULP does not make a big difference IMHO. Obviously using the correctly rounded constant in place of PI / 180 is an unconditional win, but halving performance just for a couple of bits of precision is not worth it IMHO. What would be really bad would be catastrophic failure that can result in almost unbounded error (measured in ULP) for sufficiently extreme inputs.

[huonw]

Yeah, I agree.

(I believe there is infinite (absolute and relative) error for two (positive and negative) cases for f64: let x be the float after f64::MAX / C_H, then I believe x * C should round to something slightly less than f64::MAX, but x * C_H (the naive version) overflows to infinity. I don't think this affects f32, and I don't particularly think this matters.)

[brson]

Still repros.

[Gankra]

I suggest closing this as "wont fix", honestly.

[hanna-kruppe]

One thing someone could do with little effort is replacing PI / 180 with the correctly rounded constant, which would hopefully improve accuracy slightly in some cases. Other than that, I agree with WONTFIX.