dragon.rs - source
Description: Source of the Rust file `library/core/src/num/flt2dec/strategy/dragon.rs`.

#
core/num/flt2dec/strategy/
dragon.rs

    1//! Almost direct (but slightly optimized) Rust translation of Figure 3 of "Printing
    2//! Floating-Point Numbers Quickly and Accurately"[^1].
    3//!
    4//! [^1]: Burger, R. G. and Dybvig, R. K. 1996. Printing floating-point numbers
    5//!   quickly and accurately. SIGPLAN Not. 31, 5 (May. 1996), 108-116.
    6
    7use crate::cmp::Ordering;
    8use crate::mem::MaybeUninit;
    9use crate::num::bignum::{Big32x40 as Big, Digit32 as Digit};
    10use crate::num::flt2dec::estimator::estimate_scaling_factor;
    11use crate::num::flt2dec::{Decoded, MAX_SIG_DIGITS, round_up};
    12
    13static POW10: [Digit; 10] =
    14    [1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000];
    15// precalculated arrays of `Digit`s for 5^(2^n).
    16static POW5TO16: [Digit; 2] = [0x86f26fc1, 0x23];
    17static POW5TO32: [Digit; 3] = [0x85acef81, 0x2d6d415b, 0x4ee];
    18static POW5TO64: [Digit; 5] = [0xbf6a1f01, 0x6e38ed64, 0xdaa797ed, 0xe93ff9f4, 0x184f03];
    19static POW5TO128: [Digit; 10] = [
    20    0x2e953e01, 0x3df9909, 0xf1538fd, 0x2374e42f, 0xd3cff5ec, 0xc404dc08, 0xbccdb0da, 0xa6337f19,
    21    0xe91f2603, 0x24e,
    22];
    23static POW5TO256: [Digit; 19] = [
    24    0x982e7c01, 0xbed3875b, 0xd8d99f72, 0x12152f87, 0x6bde50c6, 0xcf4a6e70, 0xd595d80f, 0x26b2716e,
    25    0xadc666b0, 0x1d153624, 0x3c42d35a, 0x63ff540e, 0xcc5573c0, 0x65f9ef17, 0x55bc28f2, 0x80dcc7f7,
    26    0xf46eeddc, 0x5fdcefce, 0x553f7,
    27];
    28
    29#[doc(hidden)]
    30pub fn mul_pow10(x: &mut Big, n: usize) -> &mut Big {
    31    debug_assert!(n < 512);
    32    // Save ourself the left shift for the smallest cases.
    33    if n < 8 {
    34        return x.mul_small(POW10[n & 7]);
    35    }
    36    // Multiply by the powers of 5 and shift the 2s in at the end.
    37    // This keeps the intermediate products smaller and faster.
    38    if n & 7 != 0 {
    39        x.mul_small(POW10[n & 7] >> (n & 7));
    40    }
    41    if n & 8 != 0 {
    42        x.mul_small(POW10[8] >> 8);
    43    }
    44    if n & 16 != 0 {
    45        x.mul_digits(&POW5TO16);
    46    }
    47    if n & 32 != 0 {
    48        x.mul_digits(&POW5TO32);
    49    }
    50    if n & 64 != 0 {
    51        x.mul_digits(&POW5TO64);
    52    }
    53    if n & 128 != 0 {
    54        x.mul_digits(&POW5TO128);
    55    }
    56    if n & 256 != 0 {
    57        x.mul_digits(&POW5TO256);
    58    }
    59    x.mul_pow2(n)
    60}
    61
    62fn div_2pow10(x: &mut Big, mut n: usize) -> &mut Big {
    63    let largest = POW10.len() - 1;
    64    while n > largest {
    65        x.div_rem_small(POW10[largest]);
    66        n -= largest;
    67    }
    68    x.div_rem_small(POW10[n] << 1);
    69    x
    70}
    71
    72// only usable when `x < 16 * scale`; `scaleN` should be `scale.mul_small(N)`
    73fn div_rem_upto_16<'a>(
    74    x: &'a mut Big,
    75    scale: &Big,
    76    scale2: &Big,
    77    scale4: &Big,
    78    scale8: &Big,
    79) -> (u8, &'a mut Big) {
    80    let mut d = 0;
    81    if *x >= *scale8 {
    82        x.sub(scale8);
    83        d += 8;
    84    }
    85    if *x >= *scale4 {
    86        x.sub(scale4);
    87        d += 4;
    88    }
    89    if *x >= *scale2 {
    90        x.sub(scale2);
    91        d += 2;
    92    }
    93    if *x >= *scale {
    94        x.sub(scale);
    95        d += 1;
    96    }
    97    debug_assert!(*x < *scale);
    98    (d, x)
    99}
    100
    101/// The shortest mode implementation for Dragon.
    102pub fn format_shortest<'a>(
    103    d: &Decoded,
    104    buf: &'a mut [MaybeUninit<u8>],
    105) -> (/*digits*/ &'a [u8], /*exp*/ i16) {
    106    // the number `v` to format is known to be:
    107    // - equal to `mant * 2^exp`;
    108    // - preceded by `(mant - 2 * minus) * 2^exp` in the original type; and
    109    // - followed by `(mant + 2 * plus) * 2^exp` in the original type.
    110    //
    111    // obviously, `minus` and `plus` cannot be zero. (for infinities, we use out-of-range values.)
    112    // also we assume that at least one digit is generated, i.e., `mant` cannot be zero too.
    113    //
    114    // this also means that any number between `low = (mant - minus) * 2^exp` and
    115    // `high = (mant + plus) * 2^exp` will map to this exact floating point number,
    116    // with bounds included when the original mantissa was even (i.e., `!mant_was_odd`).
    117
    118    assert!(d.mant > 0);
    119    assert!(d.minus > 0);
    120    assert!(d.plus > 0);
    121    assert!(d.mant.checked_add(d.plus).is_some());
    122    assert!(d.mant.checked_sub(d.minus).is_some());
    123    assert!(buf.len() >= MAX_SIG_DIGITS);
    124
    125    // `a.cmp(&b) < rounding` is `if d.inclusive {a <= b} else {a < b}`
    126    let rounding = if d.inclusive { Ordering::Greater } else { Ordering::Equal };
    127
    128    // estimate `k_0` from original inputs satisfying `10^(k_0-1) < high <= 10^(k_0+1)`.
    129    // the tight bound `k` satisfying `10^(k-1) < high <= 10^k` is calculated later.
    130    let mut k = estimate_scaling_factor(d.mant + d.plus, d.exp);
    131
    132    // convert `{mant, plus, minus} * 2^exp` into the fractional form so that:
    133    // - `v = mant / scale`
    134    // - `low = (mant - minus) / scale`
    135    // - `high = (mant + plus) / scale`
    136    let mut mant = Big::from_u64(d.mant);
    137    let mut minus = Big::from_u64(d.minus);
    138    let mut plus = Big::from_u64(d.plus);
    139    let mut scale = Big::from_small(1);
    140    if d.exp < 0 {
    141        scale.mul_pow2(-d.exp as usize);
    142    } else {
    143        mant.mul_pow2(d.exp as usize);
    144        minus.mul_pow2(d.exp as usize);
    145        plus.mul_pow2(d.exp as usize);
    146    }
    147
    148    // divide `mant` by `10^k`. now `scale / 10 < mant + plus <= scale * 10`.
    149    if k >= 0 {
    150        mul_pow10(&mut scale, k as usize);
    151    } else {
    152        mul_pow10(&mut mant, -k as usize);
    153        mul_pow10(&mut minus, -k as usize);
    154        mul_pow10(&mut plus, -k as usize);
    155    }
    156
    157    // fixup when `mant + plus > scale` (or `>=`).
    158    // we are not actually modifying `scale`, since we can skip the initial multiplication instead.
    159    // now `scale < mant + plus <= scale * 10` and we are ready to generate digits.
    160    //
    161    // note that `d[0]` *can* be zero, when `scale - plus < mant < scale`.
    162    // in this case rounding-up condition (`up` below) will be triggered immediately.
    163    if scale.cmp(mant.clone().add(&plus)) < rounding {
    164        // equivalent to scaling `scale` by 10
    165        k += 1;
    166    } else {
    167        mant.mul_small(10);
    168        minus.mul_small(10);
    169        plus.mul_small(10);
    170    }
    171
    172    // cache `(2, 4, 8) * scale` for digit generation.
    173    let mut scale2 = scale.clone();
    174    scale2.mul_pow2(1);
    175    let mut scale4 = scale.clone();
    176    scale4.mul_pow2(2);
    177    let mut scale8 = scale.clone();
    178    scale8.mul_pow2(3);
    179
    180    let mut down;
    181    let mut up;
    182    let mut i = 0;
    183    loop {
    184        // invariants, where `d[0..n-1]` are digits generated so far:
    185        // - `v = mant / scale * 10^(k-n-1) + d[0..n-1] * 10^(k-n)`
    186        // - `v - low = minus / scale * 10^(k-n-1)`
    187        // - `high - v = plus / scale * 10^(k-n-1)`
    188        // - `(mant + plus) / scale <= 10` (thus `mant / scale < 10`)
    189        // where `d[i..j]` is a shorthand for `d[i] * 10^(j-i) + ... + d[j-1] * 10 + d[j]`.
    190
    191        // generate one digit: `d[n] = floor(mant / scale) < 10`.
    192        let (d, _) = div_rem_upto_16(&mut mant, &scale, &scale2, &scale4, &scale8);
    193        debug_assert!(d < 10);
    194        buf[i] = MaybeUninit::new(b'0' + d);
    195        i += 1;
    196
    197        // this is a simplified description of the modified Dragon algorithm.
    198        // many intermediate derivations and completeness arguments are omitted for convenience.
    199        //
    200        // start with modified invariants, as we've updated `n`:
    201        // - `v = mant / scale * 10^(k-n) + d[0..n-1] * 10^(k-n)`
    202        // - `v - low = minus / scale * 10^(k-n)`
    203        // - `high - v = plus / scale * 10^(k-n)`
    204        //
    205        // assume that `d[0..n-1]` is the shortest representation between `low` and `high`,
    206        // i.e., `d[0..n-1]` satisfies both of the following but `d[0..n-2]` doesn't:
    207        // - `low < d[0..n-1] * 10^(k-n) < high` (bijectivity: digits round to `v`); and
    208        // - `abs(v / 10^(k-n) - d[0..n-1]) <= 1/2` (the last digit is correct).
    209        //
    210        // the second condition simplifies to `2 * mant <= scale`.
    211        // solving invariants in terms of `mant`, `low` and `high` yields
    212        // a simpler version of the first condition: `-plus < mant < minus`.
    213        // since `-plus < 0 <= mant`, we have the correct shortest representation
    214        // when `mant < minus` and `2 * mant <= scale`.
    215        // (the former becomes `mant <= minus` when the original mantissa is even.)
    216        //
    217        // when the second doesn't hold (`2 * mant > scale`), we need to increase the last digit.
    218        // this is enough for restoring that condition: we already know that
    219        // the digit generation guarantees `0 <= v / 10^(k-n) - d[0..n-1] < 1`.
    220        // in this case, the first condition becomes `-plus < mant - scale < minus`.
    221        // since `mant < scale` after the generation, we have `scale < mant + plus`.
    222        // (again, this becomes `scale <= mant + plus` when the original mantissa is even.)
    223        //
    224        // in short:
    225        // - stop and round `down` (keep digits as is) when `mant < minus` (or `<=`).
    226        // - stop and round `up` (increase the last digit) when `scale < mant + plus` (or `<=`).
    227        // - keep generating otherwise.
    228        down = mant.cmp(&minus) < rounding;
    229        up = scale.cmp(mant.clone().add(&plus)) < rounding;
    230        if down || up {
    231            break;
    232        } // we have the shortest representation, proceed to the rounding
    233
    234        // restore the invariants.
    235        // this makes the algorithm always terminating: `minus` and `plus` always increases,
    236        // but `mant` is clipped modulo `scale` and `scale` is fixed.
    237        mant.mul_small(10);
    238        minus.mul_small(10);
    239        plus.mul_small(10);
    240    }
    241
    242    // rounding up happens when
    243    // i) only the rounding-up condition was triggered, or
    244    // ii) both conditions were triggered and tie breaking prefers rounding up.
    245    if up && (!down || *mant.mul_pow2(1) >= scale) {
    246        // if rounding up changes the length, the exponent should also change.
    247        // it seems that this condition is very hard to satisfy (possibly impossible),
    248        // but we are just being safe and consistent here.
    249        // SAFETY: we initialized that memory above.
    250        if let Some(c) = round_up(unsafe { buf[..i].assume_init_mut() }) {
    251            buf[i] = MaybeUninit::new(c);
    252            i += 1;
    253            k += 1;
    254        }
    255    }
    256
    257    // SAFETY: we initialized that memory above.
    258    (unsafe { buf[..i].assume_init_ref() }, k)
    259}
    260
    261/// The exact and fixed mode implementation for Dragon.
    262pub fn format_exact<'a>(
    263    d: &Decoded,
    264    buf: &'a mut [MaybeUninit<u8>],
    265    limit: i16,
    266) -> (/*digits*/ &'a [u8], /*exp*/ i16) {
    267    assert!(d.mant > 0);
    268    assert!(d.minus > 0);
    269    assert!(d.plus > 0);
    270    assert!(d.mant.checked_add(d.plus).is_some());
    271    assert!(d.mant.checked_sub(d.minus).is_some());
    272
    273    // estimate `k_0` from original inputs satisfying `10^(k_0-1) < v <= 10^(k_0+1)`.
    274    let mut k = estimate_scaling_factor(d.mant, d.exp);
    275
    276    // `v = mant / scale`.
    277    let mut mant = Big::from_u64(d.mant);
    278    let mut scale = Big::from_small(1);
    279    if d.exp < 0 {
    280        scale.mul_pow2(-d.exp as usize);
    281    } else {
    282        mant.mul_pow2(d.exp as usize);
    283    }
    284
    285    // divide `mant` by `10^k`. now `scale / 10 < mant <= scale * 10`.
    286    if k >= 0 {
    287        mul_pow10(&mut scale, k as usize);
    288    } else {
    289        mul_pow10(&mut mant, -k as usize);
    290    }
    291
    292    // fixup when `mant + plus >= scale`, where `plus / scale = 10^-buf.len() / 2`.
    293    // in order to keep the fixed-size bignum, we actually use `mant + floor(plus) >= scale`.
    294    // we are not actually modifying `scale`, since we can skip the initial multiplication instead.
    295    // again with the shortest algorithm, `d[0]` can be zero but will be eventually rounded up.
    296    if *div_2pow10(&mut scale.clone(), buf.len()).add(&mant) >= scale {
    297        // equivalent to scaling `scale` by 10
    298        k += 1;
    299    } else {
    300        mant.mul_small(10);
    301    }
    302
    303    // if we are working with the last-digit limitation, we need to shorten the buffer
    304    // before the actual rendering in order to avoid double rounding.
    305    // note that we have to enlarge the buffer again when rounding up happens!
    306    let mut len = if k < limit {
    307        // oops, we cannot even produce *one* digit.
    308        // this is possible when, say, we've got something like 9.5 and it's being rounded to 10.
    309        // we return an empty buffer, with an exception of the later rounding-up case
    310        // which occurs when `k == limit` and has to produce exactly one digit.
    311        0
    312    } else if ((k as i32 - limit as i32) as usize) < buf.len() {
    313        (k - limit) as usize
    314    } else {
    315        buf.len()
    316    };
    317
    318    if len > 0 {
    319        // cache `(2, 4, 8) * scale` for digit generation.
    320        // (this can be expensive, so do not calculate them when the buffer is empty.)
    321        let mut scale2 = scale.clone();
    322        scale2.mul_pow2(1);
    323        let mut scale4 = scale.clone();
    324        scale4.mul_pow2(2);
    325        let mut scale8 = scale.clone();
    326        scale8.mul_pow2(3);
    327
    328        for i in 0..len {
    329            if mant.is_zero() {
    330                // following digits are all zeroes, we stop here
    331                // do *not* try to perform rounding! rather, fill remaining digits.
    332                for c in &mut buf[i..len] {
    333                    *c = MaybeUninit::new(b'0');
    334                }
    335                // SAFETY: we initialized that memory above.
    336                return (unsafe { buf[..len].assume_init_ref() }, k);
    337            }
    338
    339            let mut d = 0;
    340            if mant >= scale8 {
    341                mant.sub(&scale8);
    342                d += 8;
    343            }
    344            if mant >= scale4 {
    345                mant.sub(&scale4);
    346                d += 4;
    347            }
    348            if mant >= scale2 {
    349                mant.sub(&scale2);
    350                d += 2;
    351            }
    352            if mant >= scale {
    353                mant.sub(&scale);
    354                d += 1;
    355            }
    356            debug_assert!(mant < scale);
    357            debug_assert!(d < 10);
    358            buf[i] = MaybeUninit::new(b'0' + d);
    359            mant.mul_small(10);
    360        }
    361    }
    362
    363    // rounding up if we stop in the middle of digits
    364    // if the following digits are exactly 5000..., check the prior digit and try to
    365    // round to even (i.e., avoid rounding up when the prior digit is even).
    366    let order = mant.cmp(scale.mul_small(5));
    367    if order == Ordering::Greater
    368        || (order == Ordering::Equal
    369            // SAFETY: `buf[len-1]` is initialized.
    370            && len > 0 && unsafe { buf[len - 1].assume_init() } & 1 == 1)
    371    {
    372        // if rounding up changes the length, the exponent should also change.
    373        // but we've been requested a fixed number of digits, so do not alter the buffer...
    374        // SAFETY: we initialized that memory above.
    375        if let Some(c) = round_up(unsafe { buf[..len].assume_init_mut() }) {
    376            // ...unless we've been requested the fixed precision instead.
    377            // we also need to check that, if the original buffer was empty,
    378            // the additional digit can only be added when `k == limit` (edge case).
    379            k += 1;
    380            if k > limit && len < buf.len() {
    381                buf[len] = MaybeUninit::new(c);
    382                len += 1;
    383            }
    384        }
    385    }
    386
    387    // SAFETY: we initialized that memory above.
    388    (unsafe { buf[..len].assume_init_ref() }, k)
    389}