dnsl48 / fraction

[Rust] Lossless fractions and decimals; drop-in float replacement
Apache License 2.0
83 stars 25 forks source link

Overflow with small number of multiplication #33

Closed ChristopherRabotin closed 4 years ago

ChristopherRabotin commented 4 years ago

In the following example, the computation of g_rad will cause an overflow in either of the functions. This overflow error does not occur when compiled in release mode because rust does not do bound checking in release mode.


extern crate fraction;
use crate::fraction::GenericDecimal;
pub type Decimal = GenericDecimal<u128, u16>;
use std::f64::consts::TAU;

fn as_single_call() -> Decimal {
    let tt = -0.0000000000002962836210611707;

    let circ_ratio = Decimal::from(TAU) / Decimal::from(360.0);

    let multiplication = Decimal::from(35_999.050 * tt);

    let sum = Decimal::from(357.528) + multiplication;

    dbg!(sum);
    dbg!(circ_ratio);

    let g_rad = circ_ratio * sum; // <==== This is the line which fails (line 79 in my editor)
    g_rad
}

fn as_several_calls() -> Decimal {
    let tt = -0.0000000000002962836210611707;

    let g_rad = (Decimal::from(TAU) / Decimal::from(360.0))
        * (Decimal::from(357.528) + Decimal::from(35_999.050 * tt));
    g_rad
}

fn main() {
    let sgl = as_single_call();
    let svl = as_several_calls();
    println!("{}\n{}", sgl, svl);
}

Output in release mode:

RUST_BACKTRACE=1 cargo run --release
   Compiling paraod v0.1.0 (/home/chris/Workspace/rust/paraod)
    Finished release [optimized] target(s) in 0.74s
     Running `target/release/paraod`
[src/main.rs:76] sum = GenericDecimal(357.527999989334071111237863 | prec=24; Rational(Plus, Ratio { numer: 357527999989334071111237863, denom: 1000000000000000000000000 }); 357.527999989334071111237863)
[src/main.rs:77] circ_ratio = GenericDecimal(0.017453292519943 | prec=15; Rational(Plus, Ratio { numer: 3141592653589793, denom: 180000000000000000 }); 0.01745329251994329444444444444444)
0.832634730828580153746246
0.832634730828580153746246

Output in debug mode:

 RUST_BACKTRACE=1 cargo run          
   Compiling paraod v0.1.0 (/home/chris/Workspace/rust/paraod)
    Finished dev [unoptimized + debuginfo] target(s) in 0.40s
     Running `target/debug/paraod`
[src/main.rs:76] sum = GenericDecimal(357.527999989334071111237863 | prec=24; Rational(Plus, Ratio { numer: 357527999989334071111237863, denom: 1000000000000000000000000 }); 357.527999989334071111237863)
[src/main.rs:77] circ_ratio = GenericDecimal(0.017453292519943 | prec=15; Rational(Plus, Ratio { numer: 3141592653589793, denom: 180000000000000000 }); 0.01745329251994329444444444444444)
thread 'main' panicked at 'attempt to multiply with overflow', /home/chris/.rustup/toolchains/stable-x86_64-unknown-linux-gnu/lib/rustlib/src/rust/library/core/src/ops/arith.rs:323:1
stack backtrace:
   0: rust_begin_unwind
             at /rustc/18bf6b4f01a6feaf7259ba7cdae58031af1b7b39/library/std/src/panicking.rs:475
   1: core::panicking::panic_fmt
             at /rustc/18bf6b4f01a6feaf7259ba7cdae58031af1b7b39/library/core/src/panicking.rs:85
   2: core::panicking::panic
             at /rustc/18bf6b4f01a6feaf7259ba7cdae58031af1b7b39/library/core/src/panicking.rs:50
   3: <u128 as core::ops::arith::Mul>::mul
             at /home/chris/.rustup/toolchains/stable-x86_64-unknown-linux-gnu/lib/rustlib/src/rust/library/core/src/ops/arith.rs:316
   4: <num_rational::Ratio<T> as core::ops::arith::Mul>::mul
             at /home/chris/.cargo/registry/src/github.com-1ecc6299db9ec823/num-rational-0.2.4/src/lib.rs:766
   5: <fraction::fraction::GenericFraction<T> as core::ops::arith::Mul>::mul
             at /home/chris/.cargo/registry/src/github.com-1ecc6299db9ec823/fraction-0.6.3/src/fraction/mod.rs:1153
   6: <fraction::decimal::GenericDecimal<T,P> as core::ops::arith::Mul>::mul
             at /home/chris/.cargo/registry/src/github.com-1ecc6299db9ec823/fraction-0.6.3/src/decimal/mod.rs:130
   7: paraod::as_single_call
             at ./src/main.rs:79
   8: paraod::main
             at ./src/main.rs:92
   9: core::ops::function::FnOnce::call_once
             at /home/chris/.rustup/toolchains/stable-x86_64-unknown-linux-gnu/lib/rustlib/src/rust/library/core/src/ops/function.rs:227
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.
dnsl48 commented 4 years ago

Hi @ChristopherRabotin, thanks for reporting that!

To me this looks correct as we are relying on the num::Ratio implementation underneath, and it does not perform checked math by default (for performance reasons).

In the case you reported the precision of -0.0000000000002962836210611707; appears to be high enough to assume that there won't be enough space on stack for most of math operations with those numbers (as Decimals are using Fractions internally, and those rely on Rational calculations).

To work with those numbers without overflows you may switch to BigDecimal or DynaDecimal. Those types use memory allocation for working with bigger numbers, not fitting into stack.

ChristopherRabotin commented 4 years ago

Thanks for your quick response.

How come that number can't be represented without memory allocation? It fits as an f64, and if PI and TAU can be represented (although they are irrational), I would assume that this other number would be representable. What am I misunderstanding?

I'm using fraction for high precision time computations, but the structure needs to be copy-able so I can't use BigDecimal yet.

dnsl48 commented 4 years ago

I think that comes to the internal representation of our Decimal implementation. Our Decimal is just a wrapper over Fraction type, which in turn is a wrapper around num::Ratio. As such, Decimals use rationals internally for the computations. Given that, sum is internally represented as Ratio { numer: 357527999989334071111237863, denom: 1000000000000000000000000 } circ_ratio is internally represented as Ratio { numer: 3141592653589793, denom: 180000000000000000 }

I didn't manually debug that, but from what I can see, the denominators are different, so to figure out the result of the multiplication we will have to get an intermediate value of these two denoms and numerators multiplied, and those values don't fit on stack.

the structure needs to be copy-able so I can't use BigDecimal yet

Well, that feels like a limitation of this library, but it needs the intermediate values of the computations to be bigger than the computation results. One potential way to work it around might be to use BigDecimal for the computation and then downcast it to the simple Decimal if the final result fits. But obviously that requires some manual juggling in the code. Sorry it's not that simple for your use-case.

Also, there might be other rust libraries out there which are more optimized for decimal calculus (not using Rationals internally).

ChristopherRabotin commented 4 years ago

I like the idea of using BigDecimal only for this calculation and then downcasting it.

I'm using fraction for hifitime and the bug is extremely limited: if you try to convert a time to this odd time system (called Dynamic Barycentric Time, or TDB) which is within +/- 1.5 seconds from 2000-01-01T11:59:27.815064907 then it blows up in debug mode. Every other conversions works fine.

dnsl48 commented 4 years ago

using BigDecimal only for this calculation and then downcasting it

This is what DynaInt type already does. However, the type itself is not implementing Copy trait. However, you could perhaps use that type manually fetching its value with its unpack method. I guess the easiest way to do that would be to use DynaDecimal and then extract its numer/denom values initializing Decimal with unpack'ed` numer and denom values.

ChristopherRabotin commented 4 years ago

I ended up implementing it this way, and it works like a charm since this function returns an f64:

fn inner_g_rad(&self) -> f64 {
        // Let's do this computation on a big decimal and downcast later.
        use crate::fraction::BigDecimal;
        use std::f64::consts::TAU;
        let g_rad = (BigDecimal::from(TAU) / BigDecimal::from(360.0))
            * (BigDecimal::from(357.528)
                + BigDecimal::from(35_999.050 * self.as_tt_centuries_j2k()));
        let g_rad_f64 = g_rad.to_f64();

        let inner = g_rad + BigDecimal::from(0.0167 * g_rad_f64.as_ref().unwrap().sin());
        inner.to_f64().unwrap()
    }
dnsl48 commented 4 years ago

Oh, true that's much simpler if you need f64 as the result :)

ChristopherRabotin commented 4 years ago

Thanks for all your help and for being super responsive.