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

fast_posit/posit/ops/
mul.rs
