implicit_predicates_core/scratch/
product.rs

1//! The product of two expansions in an arena.
2
3use super::{Scratch, ScratchExhausted, Span};
4use crate::expansion;
5
6impl Scratch<'_> {
7    /// Writes the expansion of `a * b` and returns its span.
8    ///
9    /// The product is the sum, over the components `c` of the shorter
10    /// operand, of the longer operand scaled by `c` ([`expansion::scale`]),
11    /// accumulated with [`expansion::sum`]. It is exact, and the result is a
12    /// valid expansion of at most `2 * a.len() * b.len()` components.
13    ///
14    /// Requirement, as for [`scale`](Self::scale): every component of `a`
15    /// and of `b` lies in the product domain (see
16    /// [`eft::in_product_domain`](crate::eft::in_product_domain)), since
17    /// every product of a component of `a` by one of `b` is computed exactly
18    /// with [`eft::two_product`](crate::eft::two_product); debug builds
19    /// panic if one does not. The components of the result can lie outside
20    /// that domain, so a product of products needs the caller to bound the
21    /// magnitudes involved.
22    ///
23    /// Needs `4 * a.len() * b.len() + 2 * max(a.len(), b.len())` free
24    /// components while it works, or none if either operand is empty. The
25    /// result ends up just above the previous top, like that of every other
26    /// operation, and the working space beyond it is released.
27    ///
28    /// # Errors
29    ///
30    /// [`ScratchExhausted`] if that space is not free; the arena is
31    /// unchanged.
32    ///
33    /// # Examples
34    ///
35    /// ```
36    /// use implicit_predicates_core::Scratch;
37    ///
38    /// let tiny = 2f64.powi(-60);
39    /// let mut buf = [0.0; 32];
40    /// let mut scratch = Scratch::borrowed(&mut buf);
41    /// let a = scratch.push(&[tiny, 1.0])?; // 1 + 2^-60
42    /// let b = scratch.push(&[-tiny, 1.0])?; // 1 - 2^-60
43    /// let p = scratch.product(a, b)?;
44    /// assert_eq!(scratch.get(p), &[-tiny * tiny, 1.0]); // 1 - 2^-120
45    /// assert_eq!(scratch.used(), 4 + p.len());
46    /// assert_eq!(scratch.high_water(), 4 + 4 * 2 * 2 + 2 * 2);
47    /// # Ok::<(), implicit_predicates_core::ScratchExhausted>(())
48    /// ```
49    pub fn product(&mut self, a: Span, b: Span) -> Result<Span, ScratchExhausted> {
50        let needed = working_space(a.len, b.len).ok_or(ScratchExhausted)?;
51        self.write_result([a, b], needed, |[a, b], work| {
52            let (short, long) = if a.len() <= b.len() { (a, b) } else { (b, a) };
53            write_product(short, long, work)
54        })
55    }
56}
57
58/// The components [`Scratch::product`] needs above the top for operands of
59/// `la` and `lb` components, or `None` if that does not fit in a `u64`.
60fn working_space(la: u32, lb: u32) -> Option<u64> {
61    if la == 0 || lb == 0 {
62        return Some(0);
63    }
64    let longer = u64::from(la.max(lb));
65    (u64::from(la) * u64::from(lb))
66        .checked_mul(4)?
67        .checked_add(2 * longer)
68}
69
70/// Writes the expansion of `short * long` at the start of `work` and
71/// returns its length. `work` holds exactly [`working_space`] components,
72/// split into three zones: two of `2 * short.len() * long.len()` that take
73/// turns holding the running sum, and one of `2 * long.len()` for each
74/// partial product.
75fn write_product(short: &[f64], long: &[f64], work: &mut [f64]) -> usize {
76    if short.is_empty() {
77        return 0;
78    }
79    let zone = 2 * short.len() * long.len();
80    let (sums, partial) = work.split_at_mut(2 * zone);
81    let mut in_first_zone = true;
82    let mut len = 0;
83    for &c in short {
84        let partial_len = expansion::scale(long, c, partial);
85        let (first, second) = sums.split_at_mut(zone);
86        let (from, to) = if in_first_zone {
87            (first, second)
88        } else {
89            (second, first)
90        };
91        len = expansion::sum(&from[..len], &partial[..partial_len], to);
92        in_first_zone = !in_first_zone;
93    }
94    if !in_first_zone {
95        sums.copy_within(zone..zone + len, 0);
96    }
97    len
98}
99
100#[cfg(test)]
101mod tests {
102    use super::working_space;
103
104    #[test]
105    fn working_space_is_exact_and_never_overflows() {
106        assert_eq!(working_space(0, 7), Some(0));
107        assert_eq!(working_space(7, 0), Some(0));
108        assert_eq!(working_space(2, 3), Some(4 * 6 + 2 * 3));
109        assert_eq!(working_space(u32::MAX, 1), Some(6 * u64::from(u32::MAX)));
110        assert_eq!(working_space(u32::MAX, u32::MAX), None);
111    }
112}
implicit_predicates_core/scratch/product.rs

implicit_predicates_core/scratch/
product.rs
