`powd` overflow

[maurolacy]

Hello, and first off, thanks for the library.

We're using rust_decimal to compute a function using fixed point math, involving exp and powd.

Now, what we are noticing is that powd's range is not very large. By example, when raising a Decimal larger than 32_313_447 to a Decimal power smaller than one (0.68 in particular) using powd, we are hitting

thread ... panicked at 'Pow overflowed', .../github.com-1ecc6299db9ec823/rust_decimal-1.16.0/src/maths.rs:283:21

This doesn't make much sense, as a number raised to an exponent smaller than 1 would yield a smaller number (127_755, in this particular case).

Are you aware of this? I can always take a look at the code, but, asking first in case this is a kwown issue / there are known workarounds. What do you suggest?

Thanks,

[maurolacy]

Here's how to reproduce:

$ evcxr
Welcome to evcxr. For help, type :help
>> :dep rust_decimal = { version = "1.16", default-features = false, features = ["maths"] }
>> use rust_decimal::prelude::*;
>> let base = Decimal::new(32_313_447, 0);
>> let exponent = Decimal::new(68, 2);
>> let one = Decimal::new(1, 0);
>> base
32313447
>> exponent
0.68
>> base.powd(exponent)
127750.73491154239696209320793
>> (base + one).powd(exponent)
thread '<unnamed>' panicked at 'Pow overflowed', .../github.com-1ecc6299db9ec823/rust_decimal-1.16.0/src/maths.rs:283:21
stack backtrace:
   0: std::panicking::begin_panic
   1: <rust_decimal::decimal::Decimal as rust_decimal::maths::MathematicalOps>::powd
   2: run_user_code_11
   3: evcxr::runtime::Runtime::run_loop
   4: evcxr::runtime::runtime_hook
   5: evcxr::main
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.
Segmentation fault.
   0: evcxr::runtime::Runtime::install_crash_handlers::segfault_handler
   1: <unknown>
   2: mi_free
   3: alloc::alloc::dealloc
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/alloc/src/alloc.rs:104:14
      <alloc::alloc::Global as core::alloc::Allocator>::deallocate
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/alloc/src/alloc.rs:239:22
      alloc::alloc::box_free
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/alloc/src/alloc.rs:334:9
      panic_unwind::real_imp::cleanup
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/panic_unwind/src/gcc.rs:83:5
      __rust_panic_cleanup
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/panic_unwind/src/lib.rs:99:19
   4: std::panicking::try::cleanup
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/std/src/panicking.rs:360:42
   5: std::panicking::try::do_catch
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/std/src/panicking.rs:404:23
      std::panicking::try
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/std/src/panicking.rs:343:19
      std::panic::catch_unwind
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/std/src/panic.rs:431:14
      std::rt::lang_start_internal
             at /rustc/53cb7b09b00cbea8754ffb78e7e3cb521cb8af4b/library/std/src/rt.rs:34:21
   6: main
   7: __libc_start_main
             at /build/glibc-vjB4T1/glibc-2.28/csu/../csu/libc-start.c:308:16
   8: _start

Child process terminated with status: signal: 6
>> 

[schungx]

This doesn't make much sense, as a number raised to an exponent smaller than 1 would yield a smaller number (127_755, in this particular case).

Decimal tries to maintain decimal digits, so after a few iterations, it is possible for the number of digits to overflow, even though the number gets smaller. If you don't care about rounding errors, then you really should be doing your pow and exp with f32 or f64 instead.

[maurolacy]

Decimal tries to maintain decimal digits, so after a few iterations, it is possible for the number of digits to overflow, even though the number gets smaller. If you don't care about rounding errors, then you really should be doing your pow and exp with f32 or f64 instead.

We cannot. We are forbidden to use floating point at all, because of non-determinism, i.e. different results being obtained in different implementations / architectures. This is in the context of blockchain, so, we need for the results to be fully deterministic.

Isn't there a way to use fixed point arithmetic but truncating the less significant digits, instead of overflowing? I'm now using saturating ops (through unwrap_or), so this is a relatively minor issue. But still, it would be nice to solve it in a more elegant way. It would also probably make this library more useful.

[maurolacy]

Also, I've noticed that sqrt() doesn't suffer from this. Or at least, its range of application is order of magnitudes larger than that of powd. So, powd can for sure be made to behave; at least for exponents smaller than one.

[schungx]

We cannot. We are forbidden to use floating point at all, because of non-determinism, i.e. different results being obtained in different implementations / architectures.

I understand that the IEEE 754 standard defines deterministic calculations for addition, subtraction, multiplication, division plus the square root, meaning that they should generate deterministic output when given the exact same inputs. So your restriction is strange; different implementations / architectures should NOT give different results - they would only be wrong results, meaning that the implementation is wrong

For exponential (and also power), you start running into the Table Maker's Dilemma and sooner or later you get into rounding situations no matter what you use.

It is probably possible to make an exp or pow where, for each round, it then normalizes the number into a fixed number of significant digits before continuing... which may make it accept a wider range of inputs at the risk of having rounding errors starts accumulating.

I'm not really an expert on this, but if you want exp or pow to work over a full range of 32-bit numbers, you may need more than the 96 bits that this crate is keeping. In fact, you may need to use a arbitrary-length number library.

[maurolacy]

if you want exp or pow to work over a full range of 32-bit numbers, you may need more than the 96 bits that this crate is keeping. In fact, you may need to use a arbitrary-length number library.

Yes. In fact it would be nice if this library used i128 as the underlying type, instead of i64. But that is (literally) another issue, planned for v2.

For now, it would be good for us if powd could work smoothly up to 5 x 10^12 or so (with exponents smaller than one). This is more or less sqrt()'s current range.

[schungx]

https://github.com/paupino/rust-decimal/blob/master/src/maths.rs#L316-L323

This is where it calculates the pow by converting it into exp, and this may be the place where it overflows.

Since you know the exponent will be <1, I suggest you break down the steps and display intermediate values to find out which step it is overflowing. Then you may simply add normalize calls?

[maurolacy]

The problem is a checked_mul overflow, in https://github.com/paupino/rust-decimal/blob/master/src/maths.rs#L187.

I guess there's not much to do, except for increasing the number of bits, or improving / optimizing the impl somehow (if possible).

[schungx]

Mul overflows because you have too many significant digits. You can avoid it by trimming down the number of digits.

[paupino]

As @schungx correctly points out, it's due to the current implementation attempting to maintain precision. This "overflowing" precision is not uncommon in simple multiplication and division scenarios - i.e. typically we'd call this an "underflow" and automatically round the number to make it fit into the allocated bits. Of course, when we have more complex scenarios whereby a multiplication is just a small part of the formula you'll start losing precision later down the track each time you round - especially as smaller numbers trend to zero.

The current pow implementation was purposely taking a lazy approach by way of leveraging checked_mul. The underlying implementation of checked_mul actually uses 192 bits (e.g. 96 + 96) to calculate the multiplication product before "shrinking" it back to the required 96 bits to help handle the underflow scenarios. I remember as I was implementing powd I started to expose the raw product however the size of the change quickly spiraled out of control - consequently I used the checked_ functions as a shortcut.

Anyway, one possible solution would be to keep the 196 bit precision product in tact until the calculation is complete - i.e. only rounding to 96 bits at the last moment. Unfortunately, this requires a fair bit of refactoring to get going since other operations would also need to be able to support Buf24 data structures. The other solution (which @schungx mentions above) is to "normalize" the product as you go effectively ignoring the overflow and saturating the result (by rounding). The risk of course is that you lose precision - especially since it can trend to zero.

The better solution is to of course keep as much precision, for as long as possible. This is definitely a large chunk of work which makes it tempting to wait until v2 for (since variable precision could make this much easier). The easier solution (for now) would be to provide saturating logic in the power function - perhaps as a separate function. It'd still require a bit of fiddling around to get going but could perhaps solve your requirement.

[maurolacy]

Thanks for the detailed explanation.

I tried the "normalizing" approach mentioned above (basically, calling term.normalize_assign() before / after checked_mul) to no avail.

[maurolacy]

OK, best solution for me is to dump the last term (27) of the power series when computing the exponential. That introduces an error of roughly one decimal digit, but extends the range of powd by about 40 times.

That said, I don't think this is the right thing to do. Will just introduce an upper cut-off / saturating value in my function, and that's it.

Thanks for your time and your work.

[maurolacy]

Of course, I can always do:

        for n in 2..=27 {
            term = term.checked_mul(self.div(Decimal::new(n, 0)))?;
            let next = result + term;
            let diff = (next - result).abs();
            result = next;
            if diff <= tolerance {
                break;
            }
        }

basically distributing the powers over the factorial.

But I guess this opens up its own can of worms, in terms of precision.

Do you think this method has merit? Problems as I see it are:

• More expensive / less performant.
• Gives wildly inaccurate results for larger exponents (e^100 or so); instead of just overflowing.
• powd is less precise for large bases, too.

Nice thing is that the number of terms can be easily adjusted. And, results are pretty good for small exponents (which is my case in particular).

[schungx]

We can also use different algorithms for different ranges of inputs, picking the best one depending on the actual input values.

[schungx]

I tried the "normalizing" approach mentioned above (basically, calling term.normalize_assign() before / after checked_mul) to no avail.

You'd need to normalize to a precision that checked_mul no longer overflows/underflows.

[dovahcrow]

We got a similar issue with the following code:

>> let v = Decimal::from_f64(13.815517970299580976037625513).unwrap();
>> v.exp()

We are also using Decimal in the blockchain.

[arrtchiu]

Am I correct in thinking a division also hits a similar error because of maintaining precision? In my case I'm just simply dividing two numbers and somehow must have hit the exact right inputs to trigger this:

thread 'main' panicked at 'Division overflowed', /home/ubuntu/.cargo/registry/src/github.com-1ecc6299db9ec823/rust_decimal-1.27.0/src/arithmetic_impls.rs:218:44
stack backtrace:
   0: rust_begin_unwind
             at /rustc/2c8cc343237b8f7d5a3c3703e3a87f2eb2c54a74/library/std/src/panicking.rs:575:5
   1: core::panicking::panic_fmt
             at /rustc/2c8cc343237b8f7d5a3c3703e3a87f2eb2c54a74/library/core/src/panicking.rs:64:14

[paupino]

@arrtchiu For division to overflow it would typically mean the quotient was being divided by a very small number which is a slightly different problem. Behind the scenes, division actually does what powd should be doing - it expands the result to 192 bits and then tries to strip the underflow so that it can fit back into the 96 bit mantissa. In other words, it's unlikely to be the same problem - namely because it would be a true overflow if we couldn't fit it back into a 96 bits.

That said, can you provide the inputs that you were using? That'll help define if there is an issue or if it is expected behavior.

[arrtchiu]

Thanks @paupino, makes sense. I don't have the inputs unfortunately but if I see it again I'll instrument and try to get a repro case.