fast_posit/posit/ops/
mul.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 multiplying `x`
9  /// and `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 mul_kernel(x: Decoded<N, ES, Int>, y: Decoded<N, ES, Int>) -> (Decoded<N, ES, Int>, Int) {
17    // Multiplying two numbers in the form `frac × 2^exp` is much easier than adding them. We have
18    //
19    //   (x.frac / FRAC_DENOM * 2^x.exp) * (y.frac / FRAC_DENOM * 2^y.exp)
20    //   = (x.frac * y.frac) / FRAC_DENOM² * 2^(x.exp + y.exp)
21    //   = (x.frac * y.frac / FRAC_DENOM) / FRAC_DENOM * 2^(x.exp + y.exp)
22    //
23    // In other words: the resulting `exp` is just the sum of the `exp`s, and the `frac` is the
24    // product of the `frac`s divided by `FRAC_DENOM`. Since we know `FRAC_DENOM` = `2^FRAC_WIDTH`
25    // = `2^(Int::BITS - 2)`, we can re-arrange the expression one more time:
26    //
27    //   = (x.frac * y.frac / 2^FRAC_WIDTH) / FRAC_DENOM * 2^(x.exp + y.exp)
28    //   = ((x.frac * y.frac) >> Int::BITS) / FRAC_DENOM * 2^(x.exp + y.exp + 2)
29    //
30    // Meaning the result has
31    //
32    //   frac = (x.frac * y.frac) >> Int::BITS
33    //    exp = x.exp + y.exp + 2
34    //
35    // Only a couple other points to keep in mind:
36    //
37    //   - The multiplication must use a type with double the precision of `Int`, so that there is
38    //     no chance of overflow.
39    //   - When we shift the frac right by `Int::BITS`, we must also accumulate the lower
40    //     `Int::BITS` to `sticky`.
41    //   - The `frac` must start with `0b01` or `0b10`, i.e. it must represent a `frac` in the
42    //     range [1., 2.[ or [-2., 1.[, but the result of multiplying the `frac`s may not. When
43    //     that happens, we may need to shift 1 or 2 places left. For example: 1. × 1. = 1., but
44    //     1.5 × 1.5 = 2.25; the former is good, the latter needs an extra shift by 1 to become
45    //     1.125. Of course, if we shift the `frac` left by n places we must subtract n from `exp`.
46    //
47    // Keeping these points in mind, the final result is
48    //
49    //   frac = (x.frac * y.frac) << underflow >> Int::BITS
50    //    exp = x.exp + y.exp + 2 - underflow
51
52    use crate::underlying::Double;
53    let mul = x.frac.doubling_mul(y.frac);
54    // SAFETY: `x.frac` and `y.frac` are not 0, so their product cannot be 0; nor can it ever be MIN
55    let underflow = unsafe { mul.leading_run_minus_one() };  // Can only be 0,1,2, optimise?
56    let (frac, sticky) = (mul << underflow).components_hi_lo();
57    let exp = x.exp + y.exp + Int::ONE + Int::ONE - Int::of_u32(underflow);
58
59    (Decoded{frac, exp}, sticky)
60  }
61
62  pub(crate) fn mul(self, other: Self) -> Self {
63    if self == Self::NAR || other == Self::NAR {
64      Self::NAR
65    } else if self == Self::ZERO || other == 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      // SAFETY: `self` and `other` aren't symmetrical
72      let (result, sticky) = unsafe { Self::mul_kernel(a, b) };
73      // SAFETY: `result.is_normalised()` holds
74      unsafe { result.encode_regular_round(sticky) }
75    }
76  }
77}
78
79use core::ops::{Mul, MulAssign};
80super::mk_ops!{Mul, MulAssign, mul, mul_assign}
81
82#[cfg(test)]
83mod tests {
84  super::mk_tests!{*, *=}
85}
fast_posit/posit/ops/mul.rs

fast_posit/posit/ops/
mul.rs
