fpdec_core/
lib.rs

1// ---------------------------------------------------------------------------
2// Copyright:   (c) 2021 ff. Michael Amrhein (michael@adrhinum.de)
3// License:     This program is part of a larger application. For license
4//              details please read the file LICENSE.TXT provided together
5//              with the application.
6// ---------------------------------------------------------------------------
7// $Source: fpdec-core/src/lib.rs $
8// $Revision: 2023-05-12T15:46:50+02:00 $
9
10#![doc = include_str ! ("../README.md")]
11#![cfg_attr(not(feature = "std"), no_std)]
12// activate some rustc lints
13#![deny(non_ascii_idents)]
14#![warn(missing_debug_implementations)]
15#![warn(missing_docs)]
16#![warn(trivial_casts, trivial_numeric_casts)]
17#![warn(unused)]
18#![allow(dead_code)]
19// activate some clippy lints
20#![warn(clippy::cast_possible_truncation)]
21// #![warn(clippy::cast_possible_wrap)]
22#![warn(clippy::cast_precision_loss)]
23// #![warn(clippy::cast_sign_loss)]
24#![warn(clippy::cognitive_complexity)]
25#![warn(clippy::decimal_literal_representation)]
26#![warn(clippy::enum_glob_use)]
27#![warn(clippy::equatable_if_let)]
28#![warn(clippy::fallible_impl_from)]
29#![warn(clippy::if_not_else)]
30#![warn(clippy::if_then_some_else_none)]
31#![warn(clippy::implicit_clone)]
32#![warn(clippy::integer_division)]
33#![warn(clippy::manual_assert)]
34#![warn(clippy::match_same_arms)]
35// #![warn(clippy::mismatching_type_param_order)] TODO: enable when 1.62
36// stable
37#![warn(clippy::missing_const_for_fn)]
38#![warn(clippy::missing_errors_doc)]
39#![warn(clippy::missing_panics_doc)]
40#![warn(clippy::multiple_crate_versions)]
41// #![warn(clippy::multiple_inherent_impl)]
42#![warn(clippy::must_use_candidate)]
43#![warn(clippy::needless_pass_by_value)]
44#![warn(clippy::print_stderr)]
45#![warn(clippy::print_stdout)]
46#![warn(clippy::semicolon_if_nothing_returned)]
47#![warn(clippy::str_to_string)]
48#![warn(clippy::string_to_string)]
49#![warn(clippy::undocumented_unsafe_blocks)]
50#![warn(clippy::unicode_not_nfc)]
51#![warn(clippy::unimplemented)]
52#![warn(clippy::unseparated_literal_suffix)]
53#![warn(clippy::unused_self)]
54#![warn(clippy::unwrap_in_result)]
55#![warn(clippy::use_self)]
56#![warn(clippy::used_underscore_binding)]
57#![warn(clippy::wildcard_imports)]
58
59use core::{cmp::Ordering, ops::Neg};
60
61pub use parser::{str_to_dec, ParseDecimalError};
62pub use powers_of_ten::{checked_mul_pow_ten, mul_pow_ten, ten_pow};
63pub use rounding::{
64    i128_div_rounded, i128_mul_div_ten_pow_rounded, i128_shifted_div_rounded,
65    Round, RoundingMode,
66};
67
68mod parser;
69mod powers_of_ten;
70mod rounding;
71
72/// The maximum number of fractional decimal digits supported by `Decimal`.
73pub const MAX_N_FRAC_DIGITS: u8 = 18;
74
75#[doc(hidden)]
76#[inline]
77#[must_use]
78pub fn adjust_coeffs(x: i128, p: u8, y: i128, q: u8) -> (i128, i128) {
79    match p.cmp(&q) {
80        Ordering::Equal => (x, y),
81        Ordering::Greater => (x, mul_pow_ten(y, p - q)),
82        Ordering::Less => (mul_pow_ten(x, q - p), y),
83    }
84}
85
86#[doc(hidden)]
87#[inline]
88#[must_use]
89pub fn checked_adjust_coeffs(
90    x: i128,
91    p: u8,
92    y: i128,
93    q: u8,
94) -> (Option<i128>, Option<i128>) {
95    match p.cmp(&q) {
96        Ordering::Equal => (Some(x), Some(y)),
97        Ordering::Greater => (Some(x), checked_mul_pow_ten(y, p - q)),
98        Ordering::Less => (checked_mul_pow_ten(x, q - p), Some(y)),
99    }
100}
101
102/// Return `(q, r)` with `q = x / y` and `r = x % y`, so that `x = q * y + r`,
103/// where q is rounded against floor so that r, if non-zero, has the same sign
104/// as y and `0 <= abs(r) < abs(y)`.
105#[doc(hidden)]
106#[inline]
107#[must_use]
108#[allow(clippy::integer_division)]
109pub const fn i128_div_mod_floor(x: i128, y: i128) -> (i128, i128) {
110    let (q, r) = (x / y, x % y);
111    if (r > 0 && y < 0) || (r < 0 && y > 0) {
112        (q - 1, r + y)
113    } else {
114        (q, r)
115    }
116}
117
118// The following code is a partial copy of rust core int_log10.rs
119// TODO: remove after feature(int_log) got stable
120
121// 0 < val <= u8::MAX
122#[allow(clippy::unusual_byte_groupings)]
123#[doc(hidden)]
124#[inline]
125#[must_use]
126pub const fn u8(val: u8) -> u32 {
127    let val = val as u32;
128
129    // For better performance, avoid branches by assembling the solution
130    // in the bits above the low 8 bits.
131
132    // Adding c1 to val gives 10 in the top bits for val < 10, 11 for val >=
133    // 10
134    const C1: u32 = 0b11_00000000 - 10; // 758
135                                        // Adding c2 to val gives 01 in the top bits for val < 100, 10 for val >=
136                                        // 100
137    const C2: u32 = 0b10_00000000 - 100; // 412
138
139    // Value of top bits:
140    //            +c1  +c2  1&2
141    //     0..=9   10   01   00 = 0
142    //   10..=99   11   01   01 = 1
143    // 100..=255   11   10   10 = 2
144    ((val + C1) & (val + C2)) >> 8
145}
146
147// 0 < val < 100_000
148#[allow(clippy::unusual_byte_groupings)]
149#[doc(hidden)]
150#[inline]
151const fn less_than_5(val: u32) -> u32 {
152    // Similar to u8, when adding one of these constants to val,
153    // we get two possible bit patterns above the low 17 bits,
154    // depending on whether val is below or above the threshold.
155    const C1: u32 = 0b011_00000000000000000 - 10; // 393206
156    const C2: u32 = 0b100_00000000000000000 - 100; // 524188
157    const C3: u32 = 0b111_00000000000000000 - 1000; // 916504
158    const C4: u32 = 0b100_00000000000000000 - 10000; // 514288
159
160    // Value of top bits:
161    //                +c1  +c2  1&2  +c3  +c4  3&4   ^
162    //         0..=9  010  011  010  110  011  010  000 = 0
163    //       10..=99  011  011  011  110  011  010  001 = 1
164    //     100..=999  011  100  000  110  011  010  010 = 2
165    //   1000..=9999  011  100  000  111  011  011  011 = 3
166    // 10000..=99999  011  100  000  111  100  100  100 = 4
167    (((val + C1) & (val + C2)) ^ ((val + C3) & (val + C4))) >> 17
168}
169
170// 0 < val <= u16::MAX
171#[doc(hidden)]
172#[inline]
173#[must_use]
174pub const fn u16(val: u16) -> u32 {
175    less_than_5(val as u32)
176}
177
178// 0 < val <= u32::MAX
179#[doc(hidden)]
180#[inline]
181#[must_use]
182pub const fn u32(mut val: u32) -> u32 {
183    let mut log = 0;
184    if val >= 100_000 {
185        val /= 100_000;
186        log += 5;
187    }
188    log + less_than_5(val)
189}
190
191// 0 < val <= u64::MAX
192#[doc(hidden)]
193#[inline]
194#[must_use]
195#[allow(clippy::cast_possible_truncation)]
196pub const fn u64(mut val: u64) -> u32 {
197    let mut log = 0;
198    if val >= 10_000_000_000 {
199        val /= 10_000_000_000;
200        log += 10;
201    }
202    if val >= 100_000 {
203        val /= 100_000;
204        log += 5;
205    }
206    log + less_than_5(val as u32)
207}
208
209// 0 < val <= u128::MAX
210#[doc(hidden)]
211#[inline]
212#[must_use]
213#[allow(clippy::cast_possible_truncation)]
214pub const fn u128(mut val: u128) -> u32 {
215    let mut log = 0;
216    if val >= 100_000_000_000_000_000_000_000_000_000_000 {
217        val /= 100_000_000_000_000_000_000_000_000_000_000;
218        log += 32;
219        return log + u32(val as u32);
220    }
221    if val >= 10_000_000_000_000_000 {
222        val /= 10_000_000_000_000_000;
223        log += 16;
224    }
225    log + u64(val as u64)
226}
227
228// --- end of copied code ---
229
230#[doc(hidden)]
231#[inline]
232#[must_use]
233#[allow(clippy::cast_possible_truncation)]
234pub const fn i128_magnitude(i: i128) -> u8 {
235    // TODO: change after feature(int_log) got stable:
236    // i.log10() as u8
237    u128(i.unsigned_abs()) as u8
238}
239
240/// Return the index of the most significant bit of an u128.
241/// The given u128 must not be zero!
242fn u128_msb(mut i: u128) -> u8 {
243    debug_assert_ne!(i, 0);
244    const IDX_MAP: [u8; 16] =
245        [0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4];
246    let mut n: u8 = 0;
247    if i & 0xffffffffffffffff0000000000000000_u128 != 0 {
248        n = 64;
249        i >>= 64;
250    };
251    if i & 0x0000000000000000ffffffff00000000_u128 != 0 {
252        n += 32;
253        i >>= 32;
254    };
255    if i & 0x000000000000000000000000ffff0000_u128 != 0 {
256        n += 16;
257        i >>= 16;
258    };
259    if i & 0x0000000000000000000000000000ff00_u128 != 0 {
260        n += 8;
261        i >>= 8;
262    };
263    if i & 0x000000000000000000000000000000f0_u128 != 0 {
264        n += 4;
265        i >>= 4;
266    };
267    n + IDX_MAP[i as usize] - 1
268}
269
270#[inline(always)]
271const fn u128_hi(u: u128) -> u128 {
272    u >> 64
273}
274
275#[inline(always)]
276const fn u128_lo(u: u128) -> u128 {
277    u & 0xffffffffffffffff
278}
279
280#[inline(always)]
281const fn u128_mul_u128(x: u128, y: u128) -> (u128, u128) {
282    let xh = u128_hi(x);
283    let xl = u128_lo(x);
284    let yh = u128_hi(y);
285    let yl = u128_lo(y);
286    let mut t = xl * yl;
287    let mut rl = u128_lo(t);
288    t = xl * yh + u128_hi(t);
289    let mut rh = u128_hi(t);
290    t = xh * yl + u128_lo(t);
291    rl += u128_lo(t) << 64;
292    rh += xh * yh + u128_hi(t);
293    (rh, rl)
294}
295
296// Calculate x = x / y in place, where x = xh * 2^128 + xl, and return x % y.
297// Adapted from
298// D. E. Knuth, The Art of Computer Programming, Vol. 2, Ch. 4.3.1,
299// Exercise 16
300#[inline(always)]
301#[allow(clippy::integer_division)]
302fn u256_idiv_u64(xh: &mut u128, xl: &mut u128, y: u64) -> u128 {
303    if y == 1 {
304        return 0;
305    }
306    let y = y as u128;
307    let mut th = u128_hi(*xh);
308    let mut r = th % y;
309    let mut tl = (r << 64) + u128_lo(*xh);
310    *xh = ((th / y) << 64) + tl / y;
311    r = tl % y;
312    th = (r << 64) + u128_hi(*xl);
313    r = th % y;
314    tl = (r << 64) + u128_lo(*xl);
315    *xl = ((th / y) << 64) + tl / y;
316    tl % y
317}
318
319// Calculate x = x / y in place, where x = xh * 2^128 + xl, and return x % y.
320// Specialized version adapted from
321// Henry S. Warren, Hacker’s Delight,
322// originally found at http://www.hackersdelight.org/HDcode/divlu.c.txt.
323// That code is in turn based on Algorithm D from
324// D. E. Knuth, The Art of Computer Programming, Vol. 2, Ch. 4.3.1,
325// adapted to the special case m = 4 and n = 2 and xh < y (!).
326// The link given above does not exist anymore, but the code can still be
327// found at https://github.com/hcs0/Hackers-Delight/blob/master/divlu.c.txt.
328#[inline(always)]
329#[allow(clippy::integer_division)]
330fn u256_idiv_u128_special(xh: &mut u128, xl: &mut u128, mut y: u128) -> u128 {
331    debug_assert!(*xh < y);
332    const B: u128 = 1 << 64;
333    // Normalize dividend and divisor, so that y > 2^127 (i.e. highest bit
334    // set)
335    let n_bits = 127 - u128_msb(y);
336    y <<= n_bits;
337    let yn1 = u128_hi(y);
338    let yn0 = u128_lo(y);
339    // bits to be shifted from xl to xh:
340    let sh = if n_bits == 0 {
341        0
342    } else {
343        *xl >> (128 - n_bits)
344    };
345    let xn32 = *xh << n_bits | sh;
346    let xn10 = *xl << n_bits;
347    let xn1 = u128_hi(xn10);
348    let xn0 = u128_lo(xn10);
349    let mut q1 = xn32 / yn1;
350    let mut rhat = xn32 % yn1;
351    // Now we have
352    // q1 * yn1 + rhat = xn32
353    // so that
354    // q1 * yn1 * 2^64 + rhat * 2^64 + xn1 = xn32 * 2^64 + xn1
355    while q1 >= B || q1 * yn0 > rhat * B + xn1 {
356        q1 -= 1;
357        rhat += yn1;
358        if rhat >= B {
359            break;
360        }
361    }
362    // The loop did not change the equation given above. It was terminated if
363    // either q1 < 2^64 or rhat >= 2^64 or q1 * yn0 > rhat * 2^64 + xn1.
364    // In these cases follows:
365    // q1 * yn0 <= rhat * 2^64 + xn1, therefor
366    // q1 * yn1 * 2^64 + q1 * yn0 <= xn32 * 2^64 + xn1, and
367    // q1 * y <= xn32 * 2^64 + xn1, and
368    // xn32 * 2^64 + xn1 - q1 * y >= 0.
369    // That means that the add-back step in Knuth's algorithm is not required.
370
371    // Since the final quotient is < 2^128, this must also be true for
372    // xn32 * 2^64 + xn1 - q1 * y. Thus, in the following we can safely
373    // ignore any possible overflow in xn32 * 2^64 or q1 * y.
374    let t = xn32
375        .wrapping_mul(B)
376        .wrapping_add(xn1)
377        .wrapping_sub(q1.wrapping_mul(y));
378    let mut q0 = t / yn1;
379    rhat = t % yn1;
380    while q0 >= B || q0 * yn0 > rhat * B + xn0 {
381        q0 -= 1;
382        rhat += yn1;
383        if rhat >= B {
384            break;
385        }
386    }
387    // Write back result
388    *xh = 0;
389    *xl = q1 * B + q0;
390    // Denormalize remainder
391    (t.wrapping_mul(B)
392        .wrapping_add(xn0)
393        .wrapping_sub(q0.wrapping_mul(y)))
394        >> n_bits
395}
396
397// Calculate x = x / y in place, where x = xh * 2^128 + xl, and return x % y.
398#[inline(always)]
399#[allow(clippy::cast_possible_truncation)]
400fn u256_idiv_u128(xh: &mut u128, xl: &mut u128, y: u128) -> u128 {
401    if u128_hi(y) == 0 {
402        return u256_idiv_u64(xh, xl, u128_lo(y) as u64);
403    }
404    if *xh < y {
405        return u256_idiv_u128_special(xh, xl, y);
406    }
407    let mut t = *xh % y;
408    let r = u256_idiv_u128_special(&mut t, xl, y);
409    *xh /= y;
410    r
411}
412
413/// Return `Some<(q, r)>` with `q = (x * 10^p) / y` and `r = (x * 10^p) % y`,
414/// so that `(x * 10^p) = q * y + r`, where q is rounded against floor so that
415/// r, if non-zero, has the same sign as y and `0 <= abs(r) < abs(y)`, or
416/// return `None` if |q| > i128::MAX.
417#[doc(hidden)]
418#[must_use]
419pub fn i128_shifted_div_mod_floor(
420    x: i128,
421    p: u8,
422    y: i128,
423) -> Option<(i128, i128)> {
424    let (mut xh, mut xl) =
425        u128_mul_u128(x.unsigned_abs(), ten_pow(p) as u128);
426    let r = u256_idiv_u128(&mut xh, &mut xl, y.unsigned_abs());
427    if xh != 0 || xl > i128::MAX as u128 {
428        return None;
429    }
430    // xl <= i128::MAX, so xl as i128 is safe.
431    let mut q = xl as i128;
432    // r < y, so r as i128 is safe.
433    let mut r = r as i128;
434    if x.is_negative() {
435        if y.is_negative() {
436            r = r.neg();
437        } else {
438            q = q.neg() - 1;
439            r = y - r;
440        }
441    } else if y.is_negative() {
442        q = q.neg() - 1;
443        r -= y;
444    }
445    Some((q, r))
446}
447
448/// Return `Some<(q, r)>` with `q = (x1 * x2) / y` and `r = (x1 * x2) % y`,
449/// so that `(x1 * x2) = q * y + r`, where q is rounded against floor so that
450/// r, if non-zero, has the same sign as y and `0 <= abs(r) < abs(y)`, or
451/// return `None` if |q| > i128::MAX.
452#[doc(hidden)]
453#[must_use]
454pub fn i256_div_mod_floor(
455    x1: i128,
456    x2: i128,
457    y: i128,
458) -> Option<(i128, i128)> {
459    debug_assert!(y > 0);
460    let (mut xh, mut xl) =
461        u128_mul_u128(x1.unsigned_abs(), x2.unsigned_abs());
462    let r = u256_idiv_u128(&mut xh, &mut xl, y.unsigned_abs());
463    if xh != 0 || xl > i128::MAX as u128 {
464        return None;
465    }
466    // xl <= i128::MAX, so xl as i128 is safe.
467    let mut q = xl as i128;
468    // r < y, so r as i128 is safe.
469    let mut r = r as i128;
470    if x1.is_negative() != x2.is_negative() {
471        q = q.neg() - 1;
472        r = y - r;
473    }
474    Some((q, r))
475}