fast_posit/posit/ops/

mul.rs

use super::*;

impl<
  const N: u32,
  const ES: u32,
  Int: crate::Int,
> Posit<N, ES, Int> {
  #[inline]
  pub(crate) unsafe fn mul_kernel(x: Decoded<N, ES, Int>, y: Decoded<N, ES, Int>) -> (Decoded<N, ES, Int>, Int) {
    // Multiplying two numbers in the form `frac × 2^exp` is much easier than adding them. We have
    //
    //   (x.frac / FRAC_DENOM * 2^x.exp) * (y.frac / FRAC_DENOM * 2^y.exp)
    //   = (x.frac * y.frac) / FRAC_DENOM² * 2^(x.exp + y.exp)
    //   = (x.frac * y.frac / FRAC_DENOM) / FRAC_DENOM * 2^(x.exp + y.exp)
    //
    // In other words: the resulting `exp` is just the sum of the `exp`s, and the `frac` is the
    // product of the `frac`s divided by `FRAC_DENOM`. Since we know `FRAC_DENOM` = `2^FRAC_WIDTH`
    // = `2^(Int::BITS - 2)`, we can re-arrange the expression one more time:
    //
    //   = (x.frac * y.frac / 2^FRAC_WIDTH) / FRAC_DENOM * 2^(x.exp + y.exp)
    //   = ((x.frac * y.frac) >> Int::BITS) / FRAC_DENOM * 2^(x.exp + y.exp + 2)
    //
    // Meaning the result has
    //
    //   frac = (x.frac * y.frac) >> Int::BITS
    //    exp = x.exp + y.exp + 2
    //
    // Only a couple other points to keep in mind:
    //
    //   - The multiplication must use a type with double the precision of `Int`, so that there is
    //     no chance of overflow.
    //   - When we shift the frac right by `Int::BITS`, we must also accumulate the lower
    //     `Int::BITS` to `sticky`.
    //   - The `frac` must start with `0b01` or `0b10`, i.e. it must represent a `frac` in the
    //     range [1., 2.[ or [-2., 1.[, but the result of multiplying the `frac`s may not. When
    //     that happens, we may need to shift 1 or 2 places left. For example: 1. × 1. = 1., but
    //     1.5 × 1.5 = 2.25; the former is good, the latter needs an extra shift by 1 to become
    //     1.125. Of course, if we shift the `frac` left by n places we must subtract n from `exp`.
    //
    // Keeping these points in mind, the final result is
    //
    //   frac = (x.frac * y.frac) << underflow >> Int::BITS
    //    exp = x.exp + y.exp + 2 - underflow

    use crate::underlying::Double;
    let mul = x.frac.doubling_mul(y.frac);
    // SAFETY: `x.frac` and `y.frac` are not 0, so their product cannot be 0; nor can it ever be MIN
    let underflow = unsafe { mul.leading_run_minus_one() };  // Can only be 0,1,2, optimise?
    let (frac, sticky) = (mul << underflow).components_hi_lo();
    let exp = x.exp + y.exp + Int::ONE + Int::ONE - Int::of_u32(underflow);

    (Decoded{frac, exp}, sticky)
  }

  pub(crate) fn mul(self, other: Self) -> Self {
    if self == Self::NAR || other == Self::NAR {
      Self::NAR
    } else if self == Self::ZERO || other == Self::ZERO {
      Self::ZERO
    } else {
      // SAFETY: neither `self` nor `other` are 0 or NaR
      let a = unsafe { self.decode_regular() };
      let b = unsafe { other.decode_regular() };
      // SAFETY: `self` and `other` aren't symmetrical
      let (result, sticky) = unsafe { Self::mul_kernel(a, b) };
      // SAFETY: `result` does not have an underflowing `frac`
      unsafe { result.encode_regular_round(sticky) }
    }
  }
}

use core::ops::{Mul, MulAssign};

super::mk_ops!{Mul, MulAssign, mul, mul_assign}

#[cfg(test)]
mod tests {
  use super::*;
  use malachite::rational::Rational;

  #[allow(dead_code)]
  fn ops() {
    let mut a = crate::p32::ONE;
    let mut b = crate::p32::MINUS_ONE;
    let _ = a * b;
    let _ = &a * b;
    let _ = a * &b;
    let _ = &a * &b;
    a *= b;
    b *= &a;
  }

  /// Aux function: check that `a * b` is rounded correctly.
  fn is_correct_rounded<const N: u32, const ES: u32, Int: crate::Int>(
    a: Posit<N, ES, Int>,
    b: Posit<N, ES, Int>,
  ) -> bool
  where
    Rational: TryFrom<Posit<N, ES, Int>>,
    <Rational as TryFrom<Posit<N, ES, Int>>>::Error: core::fmt::Debug
  {
    let mul_posit = a * b;
    if let (Ok(a), Ok(b)) = (Rational::try_from(a), Rational::try_from(b)) {
      let mul_exact = a * b;
      super::rational::is_correct_rounded(mul_exact, mul_posit)
    } else {
      mul_posit == Posit::<N, ES, Int>::NAR
    }
  }

  // TODO Factor all these into a macro

  #[test]
  fn posit_10_0_exhaustive() {
    for a in Posit::<10, 0, i16>::cases_exhaustive_all() {
      for b in Posit::<10, 0, i16>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }

  #[test]
  fn posit_10_1_exhaustive() {
    for a in Posit::<10, 1, i16>::cases_exhaustive_all() {
      for b in Posit::<10, 1, i16>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }

  #[test]
  fn posit_10_2_exhaustive() {
    for a in Posit::<10, 2, i16>::cases_exhaustive_all() {
      for b in Posit::<10, 2, i16>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }

  #[test]
  fn posit_10_3_exhaustive() {
    for a in Posit::<10, 3, i16>::cases_exhaustive_all() {
      for b in Posit::<10, 3, i16>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }

  #[test]
  fn posit_8_0_exhaustive() {
    for a in Posit::<8, 0, i8>::cases_exhaustive_all() {
      for b in Posit::<8, 0, i8>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }

  #[test]
  fn p8_exhaustive() {
    for a in crate::p8::cases_exhaustive_all() {
      for b in crate::p8::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }

  use proptest::prelude::*;
  const PROPTEST_CASES: u32 = if cfg!(debug_assertions) {0x1_0000} else {0x80_0000};
  proptest!{
    #![proptest_config(ProptestConfig::with_cases(PROPTEST_CASES))]

    #[test]
    fn p16_proptest(
      a in crate::p16::cases_proptest(),
      b in crate::p16::cases_proptest(),
    ) {
      assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
    }

    #[test]
    fn p32_proptest(
      a in crate::p32::cases_proptest(),
      b in crate::p32::cases_proptest(),
    ) {
      assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
    }

    #[test]
    fn p64_proptest(
      a in crate::p64::cases_proptest(),
      b in crate::p64::cases_proptest(),
    ) {
      assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
    }
  }

  #[test]
  fn posit_3_0_exhaustive() {
    for a in Posit::<3, 0, i8>::cases_exhaustive_all() {
      for b in Posit::<3, 0, i8>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }
  #[test]
  fn posit_4_0_exhaustive() {
    for a in Posit::<4, 0, i8>::cases_exhaustive_all() {
      for b in Posit::<4, 0, i8>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }
  #[test]
  fn posit_4_1_exhaustive() {
    for a in Posit::<4, 1, i8>::cases_exhaustive_all() {
      for b in Posit::<4, 1, i8>::cases_exhaustive_all() {
        assert!(is_correct_rounded(a, b), "{:?} * {:?}", a, b)
      }
    }
  }
}