Avoid losing precision in f64

You're right @quinedot, thanks for your insights.
I've never said I'm not willing to pick some fixed precision as a minimal resolution, perhaps it is a solution indeed.

Well, I've implemented support for all primitive numbers in Rust with a simple "x as f64", but I surely could print-and-parse them as the general impl, nice.

Ah, guess I misinterpreted your objections to "multiply by some power of 10 and round".

Humm yeah @user16251, that was the simplest impl at the time.
Actually, I wanted to support sub-nanoseconds range because of calculations. Imagine one made a for loop handling 1 million items, made an Interval::now() then .elapsed() to get a Duration. Then, this results in e.g. 35ms, which if the user divides by the number of items, he will get a time in the sub-nanos range... I want to be able to print these adequately.

Hi @scottmcm.
Yep, I do exactly that, I round things based on the scale of the number, so it gets more beautiful and useful for our human eyes.

But I'm implementing the v2.0 of human-repr. In this new, unreleased version, I'll support both serde and parsing, which will allow one to make a JSON or YAML config file with human representations like timeout: 10ms or "upload_limit": 10GB.
For this I'm creating an alternate Display format impl that does not do any rounding, it should output a representation that I can parse and retrieve the exact same number back.

Well, if you want to parse a finite base-10 string from a config file, do whole-number arithmetic on it, and get the exact value back, you'd probably want to use some floating-point decimal type. Unfortunately, this is non-trivial to roll yourself, but there are a number of crates that do it, e.g., dec. (There are several more crates under the 'decimal' keyword on crates.io; many of them are for arbitrary-precision decimal numbers.)

You could loop through each of the fractional digits, printing one character at a time, and stopping if n consecutive digit of the same value is encountered. Perhaps there are better heuristics, but this is the simplest one I can think of. In general, fractions will infinitely repeat when turned into decimal, which can be reasonably detected with lookbehind.

many fractions repeat but have a much longer repeat period than just repeating the same digit, e.g. 1/7 repeats the sequence 142857: 1/7 = 0.142857142857142857142857142857142857142857142857142857142857...

1 Like

You might just want to leave the precision up to the user. The documentation for the Display trait includes instructions on how to access the precision (look at the Vector2D example).

The println! and fmt! macros know that floating point numbers have limited precision, and are designed to round the output to exclude digits that could be off solely because of representational error, which is why println!(0.2) displays cleanly even though 0.2 cannot be exactly represented in binary floating point.

However, if you (or the user of your library) are doing some floating point calculations, then this is introducing addition error in your numbers that exceeds representational error, so the macros don't hide it. You'll need to know how much error is being introduced to know how much to round, which will depend entirely on the math being performed. If picking an absolute precision isn't applicable for a general purpose tool like this, it might make sense to use a relative precision. Once you have determined the amount of error, then you can do something like this to round to the desired number of digits when printing:

fn fuzzy_format(number: f64, digits_of_error: usize) -> String {
    let available_digits = f64::DIGITS - digits_of_error;
    let integer_digits = number.log10().ceil() as isize;
    if (integer_digits >= available_digits) {
        // large number with error to the left of the decimal
        // round to desired integer precision before printing, then add zeros back
        let bad_digits = integer_digits - available_digits;
        fmt!("{:.0}", number/10.pow(bad_digits)) + "0".repeat(bad_digits);
    } else {
        // small number with error to the right of the decimal
        // print the number rounding to the maximum number of good digits, then remove the trailing zeros.
        let precision = available_digits - integer_digits;
        fmt!("{:.*}", number, precision).trim_end_matches('0');
    }
}

This allocates on the heap, but it could be modified to avoid this like alice described. I'm also being sloppy mixing signed and unsigned, so this won't compile, but hopefully the idea is useful.

Specifically, we aim to print the decimal string with the fewest significant digits which is closer to the exact float value than any other exact float value, because printing the exact decimal expansion of a binary number representation isn't helpful to anybody.

I believe the easiest and off-the-shelf solution would be to use ryu or lexical's write-to-byte-array APIs for decimalization.

4 Likes

Thank you all for the suggestions!
I'm still filtering and experimenting with them...

My most promising idea now is not re-parsing the f64 without the integer part anymore, but using module to remove it!!

fn main() {
    let n = 23e9 + 1.0f64;
    println!("   original: {}", n);
    
    println!(" first pass: {}", n / 1000.);
    println!("      fract: {} <<<", (n / 1000.).fract());
    
    println!("     module: {}", n % 1000.);
    println!("pass/module: {} <<<", (n % 1000.) / 1000.);
}

//    original: 23000000001
//  first pass: 23000000.001
//       fract: 0.000999998301267624 <<<
//      module: 1
// pass/module: 0.001 <<<

playground

This allows me to divide a much smaller number in every iteration, so it gets pretty much unhindered at the end!
My actual implementation splits the input number into its constituent integer and fractional parts (.trunc() and .fract()), and does this each pass:

fract = (fract + int % DIVISOR) / DIVISOR;
int = (int / DIVISOR).trunc() + fract.trunc();
fract = fract.fract();

In a nutshell, I use fract as an accumulator, dividing itself plus the module of the next step by the divisor, and then splitting it again between the int and fract parts... This made all my tests pass with the precision I wanted! \o/

Well, I haven't tested its performance yet, but I'm excited about it! Should be faster than any re-parsings, and I don't have to fix the precision or estimate errors at all.

P.S. I support three kinds of measurements in human-repr crate, I came up with the above for the simplest one first: counts/bytes. Now I'll have to figure out how to adapt it to both durations and throughputs.

Be careful, this can still lose stability in some cases. The only way to guarantee stability is to decimalize the original value. [playground]

fn main() {
    let n = 123456789012345.0f64;

    dbg!(n);
    for i in 1..=f64::DIGITS as i32 {
        let d = 10.0f64.powi(i);
        eprintln!();
        dbg!(d);
        dbg!(n % d);
        dbg!((n % d) / d);
        let (mut n, mut iter) = (n, 0.0);
        for _ in 1..=i {
            let m = n % 10.0;
            iter += m;
            n -= m;
            n /= 10.0;
            iter /= 10.0;
        }
        dbg!(iter);
    }
}
[src/main.rs:4] n = 123456789012345.0

[src/main.rs:8] d = 10.0
[src/main.rs:9] n % d = 5.0
[src/main.rs:10] (n % d) / d = 0.5
[src/main.rs:19] iter = 0.5

[src/main.rs:8] d = 100.0
[src/main.rs:9] n % d = 45.0
[src/main.rs:10] (n % d) / d = 0.45
[src/main.rs:19] iter = 0.45

[src/main.rs:8] d = 1000.0
[src/main.rs:9] n % d = 345.0
[src/main.rs:10] (n % d) / d = 0.345
[src/main.rs:19] iter = 0.34500000000000003

[src/main.rs:8] d = 10000.0
[src/main.rs:9] n % d = 2345.0
[src/main.rs:10] (n % d) / d = 0.2345
[src/main.rs:19] iter = 0.23450000000000001

[src/main.rs:8] d = 100000.0
[src/main.rs:9] n % d = 12345.0
[src/main.rs:10] (n % d) / d = 0.12345
[src/main.rs:19] iter = 0.12344999999999999

[src/main.rs:8] d = 1000000.0
[src/main.rs:9] n % d = 12345.0
[src/main.rs:10] (n % d) / d = 0.012345
[src/main.rs:19] iter = 0.012344999999999998

[src/main.rs:8] d = 10000000.0
[src/main.rs:9] n % d = 9012345.0
[src/main.rs:10] (n % d) / d = 0.9012345
[src/main.rs:19] iter = 0.9012344999999999

[src/main.rs:8] d = 100000000.0
[src/main.rs:9] n % d = 89012345.0
[src/main.rs:10] (n % d) / d = 0.89012345
[src/main.rs:19] iter = 0.8901234499999999

[src/main.rs:8] d = 1000000000.0
[src/main.rs:9] n % d = 789012345.0
[src/main.rs:10] (n % d) / d = 0.789012345
[src/main.rs:19] iter = 0.789012345

[src/main.rs:8] d = 10000000000.0
[src/main.rs:9] n % d = 6789012345.0
[src/main.rs:10] (n % d) / d = 0.6789012345
[src/main.rs:19] iter = 0.6789012345

[src/main.rs:8] d = 100000000000.0
[src/main.rs:9] n % d = 56789012345.0
[src/main.rs:10] (n % d) / d = 0.56789012345
[src/main.rs:19] iter = 0.56789012345

[src/main.rs:8] d = 1000000000000.0
[src/main.rs:9] n % d = 456789012345.0
[src/main.rs:10] (n % d) / d = 0.456789012345
[src/main.rs:19] iter = 0.45678901234499997

[src/main.rs:8] d = 10000000000000.0
[src/main.rs:9] n % d = 3456789012345.0
[src/main.rs:10] (n % d) / d = 0.3456789012345
[src/main.rs:19] iter = 0.34567890123449996

[src/main.rs:8] d = 100000000000000.0
[src/main.rs:9] n % d = 23456789012345.0
[src/main.rs:10] (n % d) / d = 0.23456789012345
[src/main.rs:19] iter = 0.23456789012344997

[src/main.rs:8] d = 1000000000000000.0
[src/main.rs:9] n % d = 123456789012345.0
[src/main.rs:10] (n % d) / d = 0.123456789012345
[src/main.rs:19] iter = 0.123456789012345

Alternatively, (finding and) extracting the largest unit of measure first, instead of incrementally extracting the smallest uom first, will result in better stability, but will still make up data beyond about 15 significant decimal figures. As a bonus, it'll also avoid the failure case where large - small == large. [playground]

fn main() {
    const K: f64 = 1e3;
    const M: f64 = 1e6;
    const G: f64 = 1e9;
    const T: f64 = 1e12;
    const P: f64 = 1e15;

    let mut b = dbg!(123456789123456789.0f64);
    let pb = extract(&mut b, P);
    let tb = extract(&mut b, T);
    let gb = extract(&mut b, G);
    let mb = extract(&mut b, M);
    let kb = extract(&mut b, K);
    println!("{pb}PB {tb}TB {gb}GB {mb}MB {kb}KB {b}B");
}

fn extract(n: &mut f64, d: f64) -> f64 {
    let rem = *n % d;
    let out = *n - rem;
    *n -= out;
    out / d
}
[src/main.rs:8] 123456789123456789.0f64 = 1.2345678912345678e17
123PB 456TB 789GB 123MB 456KB 784B
2 Likes

Awesome @CAD97, thank you very much!
I was tinkering with your first example here, and it is almost exactly what I did, you just subtracted the remainder instead of truncating the division as I did. It was very informative to get to see the precision being "progressively lost"......

Interesting alternative, I've never thought about extracting from the largest onwards like that... I'll surely play with it, thank you! Perhaps that method can also work for the durations and throughputs I'll work on next!