fast_posit/posit/ops/
div.rs

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

fast_posit/posit/ops/
div.rs
