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  #[inline]
9  pub(crate) unsafe fn div_kernel(x: Decoded<N, ES, Int>, y: Decoded<N, ES, Int>) -> (Decoded<N, ES, Int>, Int) {
10    // Let's use ÷ to denote true mathematical division, and / denote integer division *that rounds
11    // down* (i.e. towards negative infinity, not towards zero). To divide two numbers in the form
12    // `frac × 2^exp`, we have:
13    //
14    //   (x.frac ÷ FRAC_DENOM * 2^x.exp) ÷ (y.frac ÷ FRAC_DENOM * 2^y.exp)
15    //   = (x.frac ÷ y.frac) * 2^(x.exp - y.exp)
16    //   = (x.frac ÷ y.frac * FRAC_DENOM) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
17    //
18    // Since we know `FRAC_DENOM` = `2^FRAC_WIDTH`, we can re-arrange the expression one more
19    // time:
20    //
21    //   = (x.frac ÷ y.frac * 2^FRAC_WIDTH) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
22    //   = ((x.frac ÷ y.frac) << FRAC_WIDTH) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
23    //   = ((x.frac << FRAC_WIDTH) / y.frac) ÷ FRAC_DENOM * 2^(x.exp - y.exp)
24    //
25    // Meaning the result has
26    //
27    //   frac = (x.frac << FRAC_WIDTH) / y.frac
28    //    exp = x.exp - y.exp
29    //
30    // But this is not quite correct. This is because `(x.frac << FRAC_WIDTH) / y.frac` may
31    // underflow, which means some bits will be lost at the end. To avoid this, we compute the
32    // `underflow` first, then adjust the shift amount and the exponent accordingly.
33    //
34    //   frac = (x.frac << (FRAC_WIDTH + underflow)) / y.frac
35    //    exp = x.exp - y.exp - underflow
36
37    // TODO: The current implementation does two divisions, which is expensive. But the `underflow`
38    // can really only be [0,1,2]. Maybe we can determine this by just looking at the signs and
39    // relative magnitudes of the `frac`s, without dividing. Then we only need to do the second
40    // division.
41    let (div, _) = x.frac.shift_div_rem(y.frac, Decoded::<N, ES, Int>::FRAC_WIDTH);
42    // SAFETY: `x.frac` and `y.frac` are not 0, so `div` cannot be 0; nor can it ever be MIN
43    let underflow = unsafe { div.leading_run_minus_one() };
44
45    let (frac, sticky) = x.frac.shift_div_rem(y.frac, Decoded::<N, ES, Int>::FRAC_WIDTH + underflow);
46    let exp = x.exp - y.exp - Int::of_u32(underflow);
47
48    (Decoded{frac, exp}, sticky)
49  }
50
51  pub(crate) fn div(self, other: Self) -> Self {
52    if self == Self::NAR || other == Self::NAR || other == Self::ZERO {
53      Self::NAR
54    } else if self == Self::ZERO {
55      Self::ZERO
56    } else {
57      // SAFETY: neither `self` nor `other` are 0 or NaR
58      let a = unsafe { self.decode_regular() };
59      let b = unsafe { other.decode_regular() };
60      let (result, sticky) = unsafe { Self::div_kernel(a, b) };
61      // SAFETY: `result` does not have an underflowing `frac`
62      unsafe { result.encode_regular_round(sticky) }
63    }
64  }
65}
66
67use core::ops::{Div, DivAssign};
68
69super::mk_ops!{Div, DivAssign, div, div_assign}
70
71#[cfg(test)]
72mod tests {
73  use super::*;
74  use malachite::rational::Rational;
75
76  #[allow(dead_code)]
77  fn ops() {
78    let mut a = crate::p32::ONE;
79    let mut b = crate::p32::MINUS_ONE;
80    let _ = a / b;
81    let _ = &a / b;
82    let _ = a / &b;
83    let _ = &a / &b;
84    a /= b;
85    b /= &a;
86  }
87
88  /// Aux function: check that `a / b` is rounded correctly.
89  fn is_correct_rounded<const N: u32, const ES: u32, Int: crate::Int>(
90    a: Posit<N, ES, Int>,
91    b: Posit<N, ES, Int>,
92  ) -> bool
93  where
94    Rational: TryFrom<Posit<N, ES, Int>>,
95    <Rational as TryFrom<Posit<N, ES, Int>>>::Error: core::fmt::Debug
96  {
97    let mul_posit = a / b;
98    if let (Ok(a), Ok(b)) = (Rational::try_from(a), Rational::try_from(b))
99    && b != Rational::from(0) {
100      let mul_exact = a / b;
101      super::rational::is_correct_rounded(mul_exact, mul_posit)
102    } else {
103      mul_posit == Posit::<N, ES, Int>::NAR
104    }
105  }
106
107  // TODO Factor all these into a macro
108
109  #[test]
110  fn posit_10_0_exhaustive() {
111    for a in Posit::<10, 0, i16>::cases_exhaustive_all() {
112      for b in Posit::<10, 0, i16>::cases_exhaustive_all() {
113        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
114      }
115    }
116  }
117
118  #[test]
119  fn posit_10_1_exhaustive() {
120    for a in Posit::<10, 1, i16>::cases_exhaustive_all() {
121      for b in Posit::<10, 1, i16>::cases_exhaustive_all() {
122        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
123      }
124    }
125  }
126
127  #[test]
128  fn posit_10_2_exhaustive() {
129    for a in Posit::<10, 2, i16>::cases_exhaustive_all() {
130      for b in Posit::<10, 2, i16>::cases_exhaustive_all() {
131        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
132      }
133    }
134  }
135
136  #[test]
137  fn posit_10_3_exhaustive() {
138    for a in Posit::<10, 3, i16>::cases_exhaustive_all() {
139      for b in Posit::<10, 3, i16>::cases_exhaustive_all() {
140        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
141      }
142    }
143  }
144
145  #[test]
146  fn posit_8_0_exhaustive() {
147    for a in Posit::<8, 0, i8>::cases_exhaustive_all() {
148      for b in Posit::<8, 0, i8>::cases_exhaustive_all() {
149        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
150      }
151    }
152  }
153
154  #[test]
155  fn p8_exhaustive() {
156    for a in crate::p8::cases_exhaustive_all() {
157      for b in crate::p8::cases_exhaustive_all() {
158        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
159      }
160    }
161  }
162
163  use proptest::prelude::*;
164  const PROPTEST_CASES: u32 = if cfg!(debug_assertions) {0x1_0000} else {0x80_0000};
165  proptest!{
166    #![proptest_config(ProptestConfig::with_cases(PROPTEST_CASES))]
167
168    #[test]
169    fn p16_proptest(
170      a in crate::p16::cases_proptest(),
171      b in crate::p16::cases_proptest(),
172    ) {
173      assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
174    }
175
176    #[test]
177    fn p32_proptest(
178      a in crate::p32::cases_proptest(),
179      b in crate::p32::cases_proptest(),
180    ) {
181      assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
182    }
183
184    #[test]
185    fn p64_proptest(
186      a in crate::p64::cases_proptest(),
187      b in crate::p64::cases_proptest(),
188    ) {
189      assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
190    }
191  }
192
193  #[test]
194  fn posit_3_0_exhaustive() {
195    for a in Posit::<3, 0, i8>::cases_exhaustive_all() {
196      for b in Posit::<3, 0, i8>::cases_exhaustive_all() {
197        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
198      }
199    }
200  }
201  #[test]
202  fn posit_4_0_exhaustive() {
203    for a in Posit::<4, 0, i8>::cases_exhaustive_all() {
204      for b in Posit::<4, 0, i8>::cases_exhaustive_all() {
205        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
206      }
207    }
208  }
209  #[test]
210  fn posit_4_1_exhaustive() {
211    for a in Posit::<4, 1, i8>::cases_exhaustive_all() {
212      for b in Posit::<4, 1, i8>::cases_exhaustive_all() {
213        assert!(is_correct_rounded(a, b), "{:?} / {:?}", a, b)
214      }
215    }
216  }
217}
fast_posit/posit/ops/div.rs

fast_posit/posit/ops/
div.rs
