Rounding errors for decimals without a finite representation

[GM6VSxh91jVuHdt]

Summary

rust-decimal has rounding errors for calcluations involving numbers that do not have finite decimal representations. Other decimal libraries do not have similar issues.

rust-decimal

I am having an issue with rounding errors. As a simple example, calculating 1 / (1/3) gives an incorrect result. Code to reproduce the issue:

fn main() {
    use rust_decimal::Decimal;
    assert_eq!(
        Decimal::new(3,0),
        (Decimal::from(1) / (Decimal::from(1) / Decimal::from(3)))
    );
}

results in:

thread 'main' panicked at 'assertion failed: `(left == right)`
  left: `3`,
 right: `3.0000000000000000000000000003`', src/main.rs:3:5
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Replacing 1/3 with any repeating decimal will display the same issue.

Similarly, calculating 1 / (1/e) displays a rounding error. This would suggest that the issue is for any decimal without a finite representation.

fn main() {
    use rust_decimal::Decimal;
    assert_eq!(
        Decimal::E,
        (Decimal::from(1) / (Decimal::from(1) / Decimal::E))
    );
}

thread 'main' panicked at 'assertion failed: `(left == right)`
  left: `2.7182818284590452353602874714`,
 right: `2.7182818284590452353602874711`', src/main.rs:3:5
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Comparison with other decimal implementations

Other implementations of decimals give correct results:

BigDecimal crate

BigDecimal code that passes:

fn main() {
    use bigdecimal::BigDecimal;
    assert_eq!(
        BigDecimal::from(3),
        BigDecimal::from(1) / (BigDecimal::from(1) / BigDecimal::from(3))
    )
}

Decimal crate

Decimal code that passes:

fn main() {
    use decimal::d128;
    assert_eq!(d128!(3), d128!(1) / (d128!(1) / d128!(3)))
}

Python standard library

Python standard library:

Python 3.8.10 (default, Jun  2 2021, 10:49:15) 
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from decimal import * 
>>> Decimal(1) / (Decimal(1) / Decimal(3))
Decimal('3.000000000000000000000000000')

Fix

I have not compared the code of the libraries above, so I am not sure where rust-decimal's behaviour differs. However, the fact that they all work correctly suggest that this isn't an intractable problem. Happy to help however I usefully can.

[schungx]

I suppose if we always lose the last digit then we'd be avoiding all the rounding errors...

I suspect other libraries might simply restrict the max. number of significant digits to one fewer than the max?

[paupino]

Just to add: rust-decimal is consistent with how .NET handles this: https://dotnetfiddle.net/jOdEyU

To be quite honest, handling rational numbers is out of scope of the library as it stands. i.e. we handle rational numbers to the best of our ability to a precision of 28. For example, the same issue would exist with (1/3) * 3.

I wouldn't be opposed to another crate extending features for rational support though (e.g. rust-decimal-rational) but isn't on any short term roadmaps.

[paupino]

It looks like some of the other libraries handle the inverse logic differently indicating a significant figure rounding. e.g. if we were to use the multiplication example:

Python:

>>> from decimal import *
>>> Decimal(3) * (Decimal(1) / Decimal(3))
Decimal('0.9999999999999999999999999999')

Decimal:

fn main() {
    use decimal::d128;
    println!("{}", d128!(3) * (d128!(1) / d128!(3)))
}

outputs: 0.9999999999999999999999999999999999

Big Decimal

fn main() {
    use bigdecimal::BigDecimal;
    println!("{}", BigDecimal::from(3) * (BigDecimal::from(1) / BigDecimal::from(3)));
}

outputs: 0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999.

[paupino]

Hit "comment" too early... 😂

What I was going to suggest is that perhaps by implementing #389 we could cater for similar behavior using an opt in approach. i.e. It seems like they're handling the base decimal similarly - however are rounding differently. Of course, the only precise answer is to implement rational numbers but in lieu of that, significant figure rounding may be suitable?

[schungx]

I may say, if the library only guarantees 28 digits of precision, we should always round the result to 28 digits.

The 3.0000000000000000000000000003 in the example above has 29 digits, and that last digit is always going to have errors crept in due to rounding. In fact, all counter examples above have 29 digits of precision, so for a 28-precision library the last one is always going to be imprecise.

A 96-bit mantissa can only hold 28 full decimal digits (the 29th one incomplete), meaning that numeric results should always normalize to 10^28, throwing out the remainder (if <= 4) or adding one (if >= 5). Essentially rounding the mantissa to 28 digits.

If we always normalize all results, then comparisons and display no longer will have this problem.

[paupino]

I think there is some confusion about the max precision being 28. It doesn't limit the number of digits being used but rather constrains the equation:

m / 10^e

Where m (the mantissa) is constrained by 96 bits (i.e. 0 <= m <= 2^96) and e has the range 0 <= e <= 28.

The example in question is still valid since it is within these constraints. Ultimately it has a mantissa of 30000000000000000000000000003 and scale of 28. The mantissa of course fits within 96 bits as per:

01100000 11101111 01101011 00011010 10111010 01101111 00000111 00100011 00110000 00000000 00000000 00000011

Ultimately, because the mantissa is 96 bits, it can have more than 28 digits (base 10) representing it. i.e. The maximum represented value is 79228162514264337593543950335 which of course is 29 digits. It's also why our string buffer size "optimum" is 31 - to fit the 29 digits + possible 2 extra for a leading 0..

Regarding normalization: the library already does try to round up/down to ensure it fits within those constraints (e.g. see

rust-decimal/src/ops/div.rs

Line 593 in 3c0f76b

fn unscale_from_overflow(num: &mut Buf12, scale: i32, sticky: bool) -> Result<i32, DivError> {

). I'd be a bit nervous about doing further rounding at the library level as it stands.

Just to elaborate: the reason the result ends up the way it is is because of:

1. 1/3 is evaluated to 0.3333333333333333333333333333
2. 1 is "rescaled" to the same scale as max (i.e. 28).
3. We now have two mantissa's with identical scales: 10000000000000000000000000000 / 3333333333333333333333333333
4. The "full" answer to this of course is: 3.0000000000000000000000000003000000000000000000000000000300000000 etc
5. The library attempts to represent as much of this as possible - hence: 30000000000000000000000000003 and scale 28 or: 3.0000000000000000000000000003

It's not that the answer is incorrect, it's simply that information is lost on the first calculation. Having the library make the assumption that the number needs to be further rounded feels a little bit risky. My gut is to recommend explicit rounding so that it is clear what is happening within the constraints of the library. An accurate approach of course is to handle rational numbers explicitly. That particular approach could be achieved by num-rational if solution exactness is required.

[schungx]

Ultimately, because the mantissa is 96 bits, it can have more than 28 digits (base 10) representing it. i.e. The maximum represented value is 79228162514264337593543950335 which of course is 29 digits. It's also why our string buffer size "optimum" is 31 - to fit the 29 digits + possible 2 extra for a leading 0..

I think that's the main problem, because 2^96 as a mantissa cannot represent 29 full digits (one of the digits is at most 7, missing 8 and 9). Therefore, you technically cannot represent exact precision of all 29-digit numbers. Which then leads to the fact that, in any calculation, the 29th digit will always be suspect due to imprecise precision.

Thus using a 96-bit mantissa gives you either: a 28-digit precise number, or... a 29-digit number which has only 28 digits if it starts with 8 or 9...

It's not that the answer is incorrect, it's simply that information is lost on the first calculation.

I agree the answer is not incorrect. It is correct to the extent of the last digit of precision, as we can see.

My gut is to recommend explicit rounding so that it is clear what is happening within the constraints of the library.

This sounds like be a good way to handle this. Explicitly rounding a number to 28 digits will normalize any number, so if this step is done to the calculation result (and before any result comparison), then precise equality should be guaranteed.

So the test validations should be revised to:

fn main() {
    use rust_decimal::Decimal;

    let expected = Decimal::new(3,0);
    let actual = Decimal::from(1) / (Decimal::from(1) / Decimal::from(3));
    let normalized = actual.round_dp_with_strategy(27, MidpointAwayFromZero);

    assert_ne!(expected, actual);
    assert_eq!(expected, normalized);
}

I would suggest adding a normalize function to Decimal that simply rounds the last digit of the mantissa, instead of having the user guess how many decimal digits there are in the number (not that it would be terribly difficult)...

[paupino]

I think that's the main problem, because 2^96 as a mantissa cannot represent 29 full digits

Got it, I see what you mean - that makes sense.

I would suggest adding a normalize function to Decimal that simply rounds the last digit of the mantissa, instead of having the user guess how many decimal digits there are in the number (not that it would be terribly difficult)...

There is a normalize function already however unfortunately it'll only strip trailing zero's. In this case there is none due to that pesky 3. That being said, it's possible that the function could be adjusted to also capture the partial bit scenario as you've mentioned and normalize to 3. I'd need to look into it a bit more and would ultimately want to have a few extra tests to really battle test that behavior. But it would be useful in scenarios like this...

[schungx]

I agree that normally this would not be a problem to most actual use cases... it would just look annoying to people who dislikes seeing 3.0000000003 in print-outs.

On the other hand, however, if a user program actually depends on the fact that x = 1 / (1 / x) in exact terms, which seems to work for floats anyway... and they simply do == on these numbers, then potentially the small rounding will cause logic to fail.

So, in other words, it is probably prudent to do this just in case programmers think Decimal is always precise and uses == which doesn't get caught by clippy.

[GM6VSxh91jVuHdt]

Thank you both for your detailed responses, I am grateful for your time!

I realise now that I have misunderstood decimal floating point arithmetic. I had thought the whole reason to use crates like these was to avoid rounding errors Indeed, from the readme of this repo:

A Decimal implementation written in pure Rust suitable for financial calculations that require significant integral and fractional digits with no round-off errors.

However, reading further from the python stdlib and Knuth have made it clear that we can't ever escape rounding errors, just representation errors; you can only be precise up to a finite number of digits.

That misunderstanding cleared up, it is worth noting that Java's BigDecimal will actually throw an exception if you ask for an exact result but the result has an inifinite decimal representation. This prevents any surprises like this!

My gut is to recommend explicit rounding so that it is clear what is happening within the constraints of the library. An accurate approach of course is to handle rational numbers explicitly. That particular approach could be achieved by num-rational if solution exactness is required.

I agree - it might be worth a disclaimer/warning in the readme/docs that this crate/all decimals are precise only to x many digits, and therefore unexpected results like this may occur.

---

I'm also not sure that rounding to 28 digits will always work; consider 1 / (1 / 11000), which unrounded is 10999.9999999999999999999989. Rounded to 28 digits this doesn't give 11000.

[schungx]

Actually speaking, Decimal is exact unless one of the numbers cannot be represented by a finite number of digits in decimal. For example, 1/3 is infinite in decimal representation.

Once you have that, it is impossible not to have rounding errors.

Since infinite representations can only be created via divison, so if you don't have divison, you won't have errors.

Of course, you only have max 28 exact digits. However you can get more digits by moving to more bits. But once you have infinite digit numbers, more bits will not help.