fast_posit/posit/ops/
div.rs

1use super::*;
2
3impl<
4  const N: u32,
5  const ES: u32,
6  Int: crate::Int,
7> Posit<N, ES, Int> {
8  /// Return a [normalised](Decoded::is_normalised) `Decoded` that's the result of dividing `x` by
9  /// `y`, plus the sticky bit.
10  ///
11  /// # Safety
12  ///
13  /// `x` and `y` have to be [normalised](Decoded::is_normalised), or calling this function
14  /// is *undefined behaviour*.
15  #[inline]
16  pub(crate) unsafe fn div_kernel(x: Decoded<N, ES, Int>, y: Decoded<N, ES, Int>) -> (Decoded<N, ES, Int>, Int) {
17    // Let's use ÷ to denote true mathematical division, and / denote integer division *that rounds
18    // down* (i.e. towards negative infinity, not towards zero). To divide two numbers in the form
19    // `frac × 2^exp`, we have:
20    //
21    //   (x.frac ÷ FRAC_DENOM * 2^x.exp) ÷ (y.frac ÷ FRAC_DENOM * 2^y.exp)
22    //   = (x.frac ÷ y.frac) * 2^(x.exp - y.exp)
23    //   = (x.frac ÷ y.frac * FRAC_DENOM) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
24    //
25    // Since we know `FRAC_DENOM` = `2^FRAC_WIDTH`, we can re-arrange the expression one more
26    // time:
27    //
28    //   = (x.frac ÷ y.frac * 2^FRAC_WIDTH) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
29    //   = ((x.frac ÷ y.frac) << FRAC_WIDTH) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
30    //   = ((x.frac << FRAC_WIDTH) / y.frac) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
31    //
32    // Meaning the result has
33    //
34    //   frac = (x.frac << FRAC_WIDTH) / y.frac
35    //    exp = x.exp - y.exp
36    //
37    // But this is not quite correct. This is because `(x.frac << FRAC_WIDTH) / y.frac` may
38    // underflow, which means some bits will be lost at the end. To avoid this, we compute the
39    // `underflow` first, then adjust the shift amount and the exponent accordingly.
40    //
41    //   frac = (x.frac << (FRAC_WIDTH + underflow)) / y.frac
42    //    exp = x.exp - y.exp - underflow
43
44    // TODO: The current implementation does two divisions, which is expensive. But the `underflow`
45    // can really only be [0,1,2]. Maybe we can determine this by just looking at the signs and
46    // relative magnitudes of the `frac`s, without dividing. Then we only need to do the second
47    // division.
48    // SAFETY: `y` is normalised, so `y.frac` cannot be 0 nor -1.
49    let (div, _) = unsafe { x.frac.shift_div_rem(y.frac, Decoded::<N, ES, Int>::FRAC_WIDTH) };
50    // SAFETY: `x.frac` and `y.frac` are not 0, so `div` cannot be 0; nor can it ever be MIN.
51    let underflow = unsafe { div.leading_run_minus_one() };
52
53    // SAFETY: `y` is normalised, so `y.frac` cannot be 0 nor -1.
54    let (frac, sticky) = unsafe { x.frac.shift_div_rem(y.frac, Decoded::<N, ES, Int>::FRAC_WIDTH + underflow) };
55    let exp = x.exp - y.exp - Int::of_u32(underflow);
56
57    (Decoded{frac, exp}, sticky)
58  }
59
60  pub(crate) fn div(self, other: Self) -> Self {
61    if self == Self::NAR || other == Self::NAR || other == Self::ZERO {
62      Self::NAR
63    } else if self == Self::ZERO {
64      Self::ZERO
65    } else {
66      // SAFETY: neither `self` nor `other` are 0 or NaR
67      let a = unsafe { self.decode_regular() };
68      let b = unsafe { other.decode_regular() };
69      let (result, sticky) = unsafe { Self::div_kernel(a, b) };
70      // SAFETY: `result.is_normalised()` holds
71      unsafe { result.encode_regular_round(sticky) }
72    }
73  }
74}
75
76use core::ops::{Div, DivAssign};
77super::mk_ops!{Div, DivAssign, div, div_assign}
78
79#[cfg(test)]
80mod tests {
81  super::mk_tests!{/, /=}
82}
fast_posit/posit/ops/div.rs

fast_posit/posit/ops/
div.rs
