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}