Skip to main content

malachite_nz/integer_vector/arithmetic/
dot_sign.rs

1// Copyright © 2026 Mikhail Hogrefe
2//
3// This file is part of Malachite.
4//
5// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
6// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
7// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
8
9use crate::integer::Integer;
10use crate::integer_vector::IntegerVector;
11use crate::integer_vector::arithmetic::dot_general::vec_dot_general;
12use crate::natural::Natural;
13use core::cmp::Ordering::{self, *};
14use malachite_base::num::arithmetic::traits::{DotSign, Sign};
15use malachite_base::num::comparison::traits::{OrdProductNeg, OrdProducts, OrdProductsNeg};
16use malachite_base::num::logic::traits::SignificantBits;
17
18// What one pass over the terms finds out, for one sign: how many nonzero terms have it, and the
19// index and bit length (the sum of its factors' bit lengths) of the largest of them.
20#[derive(Clone, Copy, Default)]
21struct Side {
22    count: u64,
23    max_bits: u64,
24    max_index: usize,
25}
26
27impl Side {
28    const fn add(&mut self, index: usize, bits: u64) {
29        if self.count == 0 || bits > self.max_bits {
30            self.max_bits = bits;
31            self.max_index = index;
32        }
33        self.count += 1;
34    }
35
36    // An upper bound on the bit length of the side's sum: each of the `count` terms is below
37    // `2^max_bits`.
38    fn sum_bits_upper_bound(self) -> u64 {
39        self.max_bits + self.count.significant_bits()
40    }
41}
42
43// The sign of the dot product of `xs` and `ys`, which have the same length. Two-dimensional vectors
44// go straight to `Integer::cmp_products_neg`. Otherwise the terms are split by sign in one pass; an
45// empty side, a single term on each side, or bit lengths far enough apart decide without any
46// product, as does a largest term that outweighs all the others. Then, if the largest product is at
47// least twice as wide as any other, the dot product of the rest is computed and the largest term is
48// weighed against it by `Integer::cmp_product_neg`; otherwise the whole dot product is computed.
49fn dot_sign(xs: &[Integer], ys: &[Integer]) -> Ordering {
50    assert_eq!(
51        xs.len(),
52        ys.len(),
53        "the vectors must have the same dimension"
54    );
55    // In two dimensions the sign of `x_0 y_0 + x_1 y_1` is the comparison of `x_0 y_0` with `-x_1
56    // y_1`, which `cmp_products_neg` makes with all of its screens.
57    if let ([x_0, x_1], [y_0, y_1]) = (xs, ys) {
58        return Integer::cmp_products_neg(x_0, y_0, x_1, y_1);
59    }
60    let mut positive = Side::default();
61    let mut negative = Side::default();
62    // The bit lengths of the largest and second-largest terms of either sign, and the index of the
63    // largest.
64    let mut top_bits = 0;
65    let mut top_index = 0;
66    let mut second_bits = 0;
67    for (i, (x, y)) in xs.iter().zip(ys).enumerate() {
68        if x.abs == 0u32 || y.abs == 0u32 {
69            continue;
70        }
71        let bits = x.abs.significant_bits() + y.abs.significant_bits();
72        if x.sign == y.sign {
73            positive.add(i, bits);
74        } else {
75            negative.add(i, bits);
76        }
77        if bits > top_bits {
78            second_bits = top_bits;
79            top_bits = bits;
80            top_index = i;
81        } else if bits > second_bits {
82            second_bits = bits;
83        }
84    }
85    match (positive.count, negative.count) {
86        (0, 0) => return Equal,
87        (_, 0) => return Greater,
88        (0, _) => return Less,
89        (1, 1) => {
90            let (p, n) = (positive.max_index, negative.max_index);
91            return Natural::cmp_products(&xs[p].abs, &ys[p].abs, &xs[n].abs, &ys[n].abs);
92        }
93        _ => {}
94    }
95    // A side's largest term is at least `2^(max_bits - 2)`, so it alone outweighs a whole other
96    // side whose sum is below `2^(max_bits - 2)`.
97    if positive.sum_bits_upper_bound() + 2 <= negative.max_bits {
98        return Less;
99    } else if negative.sum_bits_upper_bound() + 2 <= positive.max_bits {
100        return Greater;
101    }
102    // The largest term is at least `2^(top_bits - 2)`, and the other terms, each below
103    // `2^second_bits`, sum to less than `2^(second_bits + bits(count - 1))`; so if the largest term
104    // is far enough above the others, its sign is the dot product's.
105    let k = top_index;
106    let others_bits_upper_bound =
107        second_bits + (positive.count + negative.count - 1).significant_bits();
108    if others_bits_upper_bound + 2 <= top_bits {
109        return if xs[k].sign == ys[k].sign {
110            Greater
111        } else {
112            Less
113        };
114    }
115    // Setting the largest term aside pays only when its product is a large share of the work, at
116    // least twice as wide as any other; then it is weighed against the rest by `cmp_product_neg`,
117    // which forms it only if its screens cannot decide. Otherwise the whole dot product is computed
118    // in one call.
119    if top_bits >= second_bits << 1 {
120        let left = vec_dot_general(None, false, &xs[..k], &ys[..k], false);
121        let rest = vec_dot_general(Some(&left), false, &xs[k + 1..], &ys[k + 1..], false);
122        Integer::cmp_product_neg(&xs[k], &ys[k], &rest)
123    } else {
124        vec_dot_general(None, false, xs, ys, false).sign()
125    }
126}
127
128impl DotSign<Self> for IntegerVector {
129    /// Determines the sign of the dot product of two [`IntegerVector`]s, taking the first by value
130    /// and the second by value.
131    ///
132    /// The result is the ordering of the dot product against 0. It is found with as few
133    /// multiplications as possible. For two-dimensional vectors it is
134    /// [`Integer::cmp_products_neg`](malachite_base::num::comparison::traits::OrdProductsNeg) on
135    /// the two pairs of elements. Otherwise the terms are split by sign, and when one side is
136    /// empty, when there is one term on each side (compared through
137    /// [`Natural::cmp_products`](malachite_base::num::comparison::traits::OrdProducts)), or when
138    /// the bit lengths of the largest terms on the two sides are far enough apart, no product is
139    /// formed at all; nor is one when the largest term alone outweighs all the others. Otherwise,
140    /// if the largest term's product is at least twice as wide as any other, the dot product of the
141    /// other terms is computed and the largest term is weighed against it by
142    /// [`Integer::cmp_product_neg`](malachite_base::num::comparison::traits::OrdProductNeg), which
143    /// forms that product only if it must; if not, the whole dot product is computed.
144    ///
145    /// $$
146    /// f(v, w) = \operatorname{sgn}\left ( \sum_{i=0}^{n-1} v_iw_i \right ),
147    /// $$
148    ///
149    /// where $n$ is the dimension of $v$ and $w$.
150    ///
151    /// # Worst-case complexity
152    /// $T(n, m) = O(nm \log m \log\log m)$
153    ///
154    /// $M(m) = O(m \log m)$
155    ///
156    /// where $T$ is time, $M$ is additional memory, $n$ is the dimension of the vectors, and $m$ is
157    /// the largest number of significant bits of any element.
158    ///
159    /// # Panics
160    /// Panics if the vectors have different dimensions.
161    ///
162    /// # Examples
163    /// ```
164    /// use core::cmp::Ordering::*;
165    /// use core::str::FromStr;
166    /// use malachite_base::num::arithmetic::traits::DotSign;
167    /// use malachite_nz::integer_vector::IntegerVector;
168    ///
169    /// let v = IntegerVector::from_str("(3, 7)").unwrap();
170    /// let w = IntegerVector::from_str("(8, -3)").unwrap();
171    /// assert_eq!(v.clone().dot_sign(w.clone()), Greater);
172    /// let w = IntegerVector::from_str("(7, -3)").unwrap();
173    /// assert_eq!(v.clone().dot_sign(w.clone()), Equal);
174    /// ```
175    #[inline]
176    fn dot_sign(self, other: Self) -> Ordering {
177        dot_sign(&self.elements, &other.elements)
178    }
179}
180
181impl DotSign<&Self> for IntegerVector {
182    /// Determines the sign of the dot product of two [`IntegerVector`]s, taking the first by value
183    /// and the second by reference.
184    ///
185    /// The result is the ordering of the dot product against 0. It is found with as few
186    /// multiplications as possible. For two-dimensional vectors it is
187    /// [`Integer::cmp_products_neg`](malachite_base::num::comparison::traits::OrdProductsNeg) on
188    /// the two pairs of elements. Otherwise the terms are split by sign, and when one side is
189    /// empty, when there is one term on each side (compared through
190    /// [`Natural::cmp_products`](malachite_base::num::comparison::traits::OrdProducts)), or when
191    /// the bit lengths of the largest terms on the two sides are far enough apart, no product is
192    /// formed at all; nor is one when the largest term alone outweighs all the others. Otherwise,
193    /// if the largest term's product is at least twice as wide as any other, the dot product of the
194    /// other terms is computed and the largest term is weighed against it by
195    /// [`Integer::cmp_product_neg`](malachite_base::num::comparison::traits::OrdProductNeg), which
196    /// forms that product only if it must; if not, the whole dot product is computed.
197    ///
198    /// $$
199    /// f(v, w) = \operatorname{sgn}\left ( \sum_{i=0}^{n-1} v_iw_i \right ),
200    /// $$
201    ///
202    /// where $n$ is the dimension of $v$ and $w$.
203    ///
204    /// # Worst-case complexity
205    /// $T(n, m) = O(nm \log m \log\log m)$
206    ///
207    /// $M(m) = O(m \log m)$
208    ///
209    /// where $T$ is time, $M$ is additional memory, $n$ is the dimension of the vectors, and $m$ is
210    /// the largest number of significant bits of any element.
211    ///
212    /// # Panics
213    /// Panics if the vectors have different dimensions.
214    ///
215    /// # Examples
216    /// ```
217    /// use core::cmp::Ordering::*;
218    /// use core::str::FromStr;
219    /// use malachite_base::num::arithmetic::traits::DotSign;
220    /// use malachite_nz::integer_vector::IntegerVector;
221    ///
222    /// let v = IntegerVector::from_str("(3, 7)").unwrap();
223    /// let w = IntegerVector::from_str("(8, -3)").unwrap();
224    /// assert_eq!(v.clone().dot_sign(&w), Greater);
225    /// let w = IntegerVector::from_str("(7, -3)").unwrap();
226    /// assert_eq!(v.clone().dot_sign(&w), Equal);
227    /// ```
228    #[inline]
229    fn dot_sign(self, other: &Self) -> Ordering {
230        dot_sign(&self.elements, &other.elements)
231    }
232}
233
234impl DotSign<IntegerVector> for &IntegerVector {
235    /// Determines the sign of the dot product of two [`IntegerVector`]s, taking the first by
236    /// reference and the second by value.
237    ///
238    /// The result is the ordering of the dot product against 0. It is found with as few
239    /// multiplications as possible. For two-dimensional vectors it is
240    /// [`Integer::cmp_products_neg`](malachite_base::num::comparison::traits::OrdProductsNeg) on
241    /// the two pairs of elements. Otherwise the terms are split by sign, and when one side is
242    /// empty, when there is one term on each side (compared through
243    /// [`Natural::cmp_products`](malachite_base::num::comparison::traits::OrdProducts)), or when
244    /// the bit lengths of the largest terms on the two sides are far enough apart, no product is
245    /// formed at all; nor is one when the largest term alone outweighs all the others. Otherwise,
246    /// if the largest term's product is at least twice as wide as any other, the dot product of the
247    /// other terms is computed and the largest term is weighed against it by
248    /// [`Integer::cmp_product_neg`](malachite_base::num::comparison::traits::OrdProductNeg), which
249    /// forms that product only if it must; if not, the whole dot product is computed.
250    ///
251    /// $$
252    /// f(v, w) = \operatorname{sgn}\left ( \sum_{i=0}^{n-1} v_iw_i \right ),
253    /// $$
254    ///
255    /// where $n$ is the dimension of $v$ and $w$.
256    ///
257    /// # Worst-case complexity
258    /// $T(n, m) = O(nm \log m \log\log m)$
259    ///
260    /// $M(m) = O(m \log m)$
261    ///
262    /// where $T$ is time, $M$ is additional memory, $n$ is the dimension of the vectors, and $m$ is
263    /// the largest number of significant bits of any element.
264    ///
265    /// # Panics
266    /// Panics if the vectors have different dimensions.
267    ///
268    /// # Examples
269    /// ```
270    /// use core::cmp::Ordering::*;
271    /// use core::str::FromStr;
272    /// use malachite_base::num::arithmetic::traits::DotSign;
273    /// use malachite_nz::integer_vector::IntegerVector;
274    ///
275    /// let v = IntegerVector::from_str("(3, 7)").unwrap();
276    /// let w = IntegerVector::from_str("(8, -3)").unwrap();
277    /// assert_eq!((&v).dot_sign(w.clone()), Greater);
278    /// let w = IntegerVector::from_str("(7, -3)").unwrap();
279    /// assert_eq!((&v).dot_sign(w.clone()), Equal);
280    /// ```
281    #[inline]
282    fn dot_sign(self, other: IntegerVector) -> Ordering {
283        dot_sign(&self.elements, &other.elements)
284    }
285}
286
287impl DotSign<&IntegerVector> for &IntegerVector {
288    /// Determines the sign of the dot product of two [`IntegerVector`]s, taking the first by
289    /// reference and the second by reference.
290    ///
291    /// The result is the ordering of the dot product against 0. It is found with as few
292    /// multiplications as possible. For two-dimensional vectors it is
293    /// [`Integer::cmp_products_neg`](malachite_base::num::comparison::traits::OrdProductsNeg) on
294    /// the two pairs of elements. Otherwise the terms are split by sign, and when one side is
295    /// empty, when there is one term on each side (compared through
296    /// [`Natural::cmp_products`](malachite_base::num::comparison::traits::OrdProducts)), or when
297    /// the bit lengths of the largest terms on the two sides are far enough apart, no product is
298    /// formed at all; nor is one when the largest term alone outweighs all the others. Otherwise,
299    /// if the largest term's product is at least twice as wide as any other, the dot product of the
300    /// other terms is computed and the largest term is weighed against it by
301    /// [`Integer::cmp_product_neg`](malachite_base::num::comparison::traits::OrdProductNeg), which
302    /// forms that product only if it must; if not, the whole dot product is computed.
303    ///
304    /// $$
305    /// f(v, w) = \operatorname{sgn}\left ( \sum_{i=0}^{n-1} v_iw_i \right ),
306    /// $$
307    ///
308    /// where $n$ is the dimension of $v$ and $w$.
309    ///
310    /// # Worst-case complexity
311    /// $T(n, m) = O(nm \log m \log\log m)$
312    ///
313    /// $M(m) = O(m \log m)$
314    ///
315    /// where $T$ is time, $M$ is additional memory, $n$ is the dimension of the vectors, and $m$ is
316    /// the largest number of significant bits of any element.
317    ///
318    /// # Panics
319    /// Panics if the vectors have different dimensions.
320    ///
321    /// # Examples
322    /// ```
323    /// use core::cmp::Ordering::*;
324    /// use core::str::FromStr;
325    /// use malachite_base::num::arithmetic::traits::DotSign;
326    /// use malachite_nz::integer_vector::IntegerVector;
327    ///
328    /// let v = IntegerVector::from_str("(3, 7)").unwrap();
329    /// let w = IntegerVector::from_str("(8, -3)").unwrap();
330    /// assert_eq!((&v).dot_sign(&w), Greater);
331    /// let w = IntegerVector::from_str("(7, -3)").unwrap();
332    /// assert_eq!((&v).dot_sign(&w), Equal);
333    /// ```
334    #[inline]
335    fn dot_sign(self, other: &IntegerVector) -> Ordering {
336        dot_sign(&self.elements, &other.elements)
337    }
338}