malachite_float/float/arithmetic/dot.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::InnerFloat::{Finite, Infinity, NaN, Zero};
10use crate::emulate_float_slice_float_slice_to_float_fn;
11use crate::float::arithmetic::sum::{complete_sum_result, max_prec, update_zero_sign};
12use crate::{
13 Float, float_infinity, float_nan, float_negative_infinity, float_negative_zero, float_zero,
14};
15use alloc::vec::Vec;
16use core::cmp::Ordering::{self, *};
17use malachite_base::num::basic::floats::PrimitiveFloat;
18use malachite_base::num::conversion::traits::ExactFrom;
19use malachite_base::num::logic::traits::SignificantBits;
20use malachite_base::rounding_modes::RoundingMode::{self, *};
21use malachite_nz::natural::Natural;
22use malachite_nz::natural::arithmetic::float::sum::{FloatSumInput, sum_float_significands};
23
24/// Computes the dot product of two equal-length slices of primitive floats, with a single rounding.
25///
26/// The result is correctly rounded to the nearest value: the products are exact, the sum is
27/// computed as if in infinite precision, and only a single rounding is performed, at the end. This
28/// includes gradual underflow: results in the subnormal range are correctly rounded to their
29/// reduced precisions. Intermediate overflow and underflow cannot occur.
30///
31/// $$
32/// f((x_i)_ {i=0}^{n-1}, (y_i)_ {i=0}^{n-1}) = \sum_ {i=0}^{n-1} x_i y_i + \varepsilon.
33/// $$
34/// - If $\sum_ {i=0}^{n-1} x_i y_i$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored or
35/// assumed to be 0.
36/// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, then $|\varepsilon| \leq
37/// 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p}$, where $p$ is the precision of the
38/// output (typically 24 if `T` is a [`f32`] and 53 if `T` is a [`f64`], but less if the output is
39/// subnormal).
40///
41/// See [`Float::dot_prec_round`] for a description of the special cases, which follow the rules of
42/// multiplication for each term and the rules of addition for their combination.
43///
44/// If the result overflows, $\pm\infty$ is returned, and if it underflows, $\pm0.0$ is returned.
45///
46/// # Worst-case complexity
47/// $T(n) = O(n)$
48///
49/// $M(n) = O(n)$
50///
51/// where $T$ is time, $M$ is additional memory, and $n$ is `xs.len()`: the products are
52/// constant-size, and a primitive float's exponent range is bounded, so the summation window is
53/// repositioned only a constant number of times.
54///
55/// # Panics
56/// Panics if `xs` and `ys` have different lengths.
57///
58/// # Examples
59/// ```
60/// use malachite_base::num::float::NiceFloat;
61/// use malachite_float::float::arithmetic::dot::primitive_float_dot;
62///
63/// // A naive fold overflows on the first product; the correctly-rounded dot product does not.
64/// let xs = [1.0e300f64, 1.0e300];
65/// let ys = [1.0e300f64, -1.0e300];
66/// assert_eq!(NiceFloat(primitive_float_dot(&xs, &ys)), NiceFloat(0.0));
67/// ```
68#[allow(clippy::type_repetition_in_bounds)]
69#[inline]
70pub fn primitive_float_dot<T: PrimitiveFloat>(xs: &[T], ys: &[T]) -> T
71where
72 Float: From<T> + PartialOrd<T>,
73 for<'a> T: ExactFrom<&'a Float>,
74{
75 emulate_float_slice_float_slice_to_float_fn(Float::dot_prec, xs, ys)
76}
77
78// A Float's finite fields: the sign, the exponent, the precision, and a reference to the
79// significand.
80fn parts(f: &Float) -> (bool, i32, u64, &Natural) {
81 let Float(Finite {
82 sign,
83 exponent,
84 precision,
85 significand,
86 }) = f
87 else {
88 unreachable!()
89 };
90 (*sign, *exponent, *precision, significand)
91}
92
93impl Float {
94 /// Computes the dot product of two equal-length slices of [`Float`]s, rounding the result to
95 /// the specified precision and with the specified rounding mode. An [`Ordering`] is also
96 /// returned, indicating whether the rounded dot product is less than, equal to, or greater than
97 /// the exact dot product. Although `NaN`s are not comparable to any [`Float`], whenever this
98 /// function returns a `NaN` it also returns `Equal`.
99 ///
100 /// The products are never rounded, and only a single rounding is performed, at the end: the
101 /// result is the correctly-rounded exact dot product. Intermediate overflow and underflow
102 /// cannot occur: each product is computed exactly at the significand level, with its exponent
103 /// tracked over a range twice as wide as a [`Float`]'s, and only the final sum is subject to
104 /// the exponent range check. (MPFR's `mpfr_dot`, which is documented as experimental, computes
105 /// each product at full precision and requires those multiplications to be exact, so it does
106 /// not handle inputs whose products leave the exponent range.)
107 ///
108 /// See [`RoundingMode`] for a description of the possible rounding modes.
109 ///
110 /// $$
111 /// f((x_i)_ {i=0}^{n-1}, (y_i)_ {i=0}^{n-1}, p, m) = \sum_ {i=0}^{n-1} x_i y_i +
112 /// \varepsilon.
113 /// $$
114 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored
115 /// or assumed to be 0.
116 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, and $m$ is not `Nearest`, then
117 /// $|\varepsilon| < 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p+1}$.
118 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, and $m$ is `Nearest`, then
119 /// $|\varepsilon| \leq 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p}$.
120 ///
121 /// If the output has a precision, it is `prec`.
122 ///
123 /// Each term $x_iy_i$ follows the rules of [`Float`] multiplication, and the terms are then
124 /// combined following the rules of [`Float`] addition:
125 /// - The dot product of empty slices is $0.0$.
126 /// - If any term is a `NaN` — because an input is `NaN`, or because a zero is paired with an
127 /// infinity — the dot product is `NaN`.
128 /// - If two infinite terms have different signs, the dot product is `NaN`. Otherwise, if any
129 /// term is infinite, the dot product is an infinity of that sign.
130 /// - If every term is a zero and all the terms have the same sign, the dot product is a zero of
131 /// that sign. If they do not all have the same sign, the dot product is $0.0$, unless $m$ is
132 /// `Floor`, in which case it is $-0.0$.
133 /// - If some terms are nonzero but the exact dot product is zero, the dot product is $0.0$,
134 /// unless $m$ is `Floor`, in which case it is $-0.0$.
135 ///
136 /// Overflow and underflow:
137 /// - If $f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)\geq 2^{2^{30}-1}$ and $m$ is `Ceiling`,
138 /// `Up`, or `Nearest`, $\infty$ is returned instead.
139 /// - If $f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)\geq 2^{2^{30}-1}$ and $m$ is `Floor` or
140 /// `Down`, $(1-(1/2)^p)2^{2^{30}-1}$ is returned instead.
141 /// - If $f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)\leq -2^{2^{30}-1}$ and $m$ is `Floor`,
142 /// `Up`, or `Nearest`, $-\infty$ is returned instead.
143 /// - If $f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)\leq -2^{2^{30}-1}$ and $m$ is `Ceiling`
144 /// or `Down`, $-(1-(1/2)^p)2^{2^{30}-1}$ is returned instead.
145 /// - If $0<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<2^{-2^{30}}$, and $m$ is `Floor` or
146 /// `Down`, $0.0$ is returned instead.
147 /// - If $0<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<2^{-2^{30}}$, and $m$ is `Ceiling` or
148 /// `Up`, $2^{-2^{30}}$ is returned instead.
149 /// - If $0<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)\leq2^{-2^{30}-1}$, and $m$ is
150 /// `Nearest`, $0.0$ is returned instead.
151 /// - If $2^{-2^{30}-1}<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<2^{-2^{30}}$, and $m$ is
152 /// `Nearest`, $2^{-2^{30}}$ is returned instead.
153 /// - If $-2^{-2^{30}}<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<0$, and $m$ is `Ceiling` or
154 /// `Down`, $-0.0$ is returned instead.
155 /// - If $-2^{-2^{30}}<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<0$, and $m$ is `Floor` or
156 /// `Up`, $-2^{-2^{30}}$ is returned instead.
157 /// - If $-2^{-2^{30}-1}\leq f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<0$, and $m$ is
158 /// `Nearest`, $-0.0$ is returned instead.
159 /// - If $-2^{-2^{30}}<f((x_i)_ {i=0}^{n-1},(y_i)_ {i=0}^{n-1},p,m)<-2^{-2^{30}-1}$, and $m$ is
160 /// `Nearest`, $-2^{-2^{30}}$ is returned instead.
161 ///
162 /// If you know you'll be using `Nearest`, consider using [`Float::dot_prec`] instead. If you
163 /// know that your target precision is the maximum of the precisions of the inputs, consider
164 /// using [`Float::dot_round`] instead. If both of these things are true, consider using
165 /// [`Float::dot`] instead.
166 ///
167 /// # Worst-case complexity
168 /// $T(n, m, p) = O(n + m (n + p) + m \log m \log\log m)$
169 ///
170 /// $M(n, m, p) = O(n + p + m \log m)$
171 ///
172 /// where $T$ is time, $M$ is additional memory, $n$ is `xs.len()`, $m$ is the sum of the
173 /// significant bits of the elements of `xs` and `ys`, and $p$ is `prec`: each term is an exact
174 /// significand product (mul-class in the pair's bits), and the terms then feed the summation
175 /// kernel, which inherits the summation bound.
176 ///
177 /// # Panics
178 /// Panics if `prec` is zero, if `xs` and `ys` have different lengths, or if `rm` is `Exact` and
179 /// the exact dot product is not exactly representable with `prec` bits.
180 ///
181 /// # Examples
182 /// ```
183 /// use malachite_base::num::arithmetic::traits::PowerOf2;
184 /// use malachite_base::num::basic::traits::{One, Two};
185 /// use malachite_base::rounding_modes::RoundingMode::*;
186 /// use malachite_float::Float;
187 /// use std::cmp::Ordering::*;
188 ///
189 /// let xs = [Float::ONE, Float::TWO, Float::from(3)];
190 /// let ys = [Float::from(4), Float::from(5), Float::from(6)];
191 /// let (dot, o) = Float::dot_prec_round(&xs, &ys, 10, Floor);
192 /// assert_eq!(dot.to_string(), "32.000");
193 /// assert_eq!(o, Equal);
194 ///
195 /// // 0.25 * 0.25 + 2 * 5 = 10.0625
196 /// let xs = [Float::power_of_2(-2i64), Float::TWO];
197 /// let ys = [Float::power_of_2(-2i64), Float::from(5)];
198 /// let (dot, o) = Float::dot_prec_round(&xs, &ys, 3, Floor);
199 /// assert_eq!(dot.to_string(), "10.0");
200 /// assert_eq!(o, Less);
201 ///
202 /// let (dot, o) = Float::dot_prec_round(&xs, &ys, 3, Ceiling);
203 /// assert_eq!(dot.to_string(), "12.0");
204 /// assert_eq!(o, Greater);
205 /// ```
206 pub fn dot_prec_round(
207 xs: &[Self],
208 ys: &[Self],
209 prec: u64,
210 rm: RoundingMode,
211 ) -> (Self, Ordering) {
212 // The dot product is correctly rounded: the products are never rounded, and only a single
213 // rounding is performed, at the end. Unlike MPFR's experimental `mpfr_dot`, intermediate
214 // overflow and underflow cannot occur: each product is computed exactly at the significand
215 // level, with its exponent tracked in an `i64` (twice the `Float` exponent range fits
216 // comfortably), and only the final sum is subject to the exponent range check.
217 //
218 // The `Exact` rounding mode is handled by computing with `Nearest` and panicking if the
219 // result is inexact.
220 assert_ne!(prec, 0);
221 assert_eq!(
222 xs.len(),
223 ys.len(),
224 "dot product requires slices of equal length"
225 );
226 let n = xs.len();
227 if n == 0 {
228 return (float_zero!(), Equal);
229 } else if n == 1 {
230 return xs[0].mul_prec_round_ref_ref(&ys[0], prec, rm);
231 }
232 // Classify each term x * y according to the multiplication rules, then combine the terms
233 // according to the addition rules, determining the sign of an infinite result, the sign of
234 // an all-zero result, and the regular terms.
235 let mut sign_inf = 0i8;
236 let mut sign_zero = 0i8;
237 let mut regulars: Vec<(&Self, &Self)> = Vec::new();
238 for (x, y) in xs.iter().zip(ys.iter()) {
239 if x.is_nan() || y.is_nan() {
240 return (float_nan!(), Equal);
241 }
242 let term_sign = if x.is_sign_negative() == y.is_sign_negative() {
243 1
244 } else {
245 -1
246 };
247 if x.is_infinite() || y.is_infinite() {
248 if *x == 0u32 || *y == 0u32 {
249 // A zero times an infinity is NaN.
250 return (float_nan!(), Equal);
251 }
252 if sign_inf == 0 {
253 sign_inf = term_sign;
254 } else if sign_inf != term_sign {
255 // Infinite terms of opposite signs add to NaN.
256 return (float_nan!(), Equal);
257 }
258 } else if *x == 0u32 || *y == 0u32 {
259 if regulars.is_empty() {
260 // This choice is sticky when new zeros are considered.
261 update_zero_sign(&mut sign_zero, term_sign, rm);
262 }
263 } else {
264 regulars.push((x, y));
265 }
266 }
267 // At this point the result cannot be NaN.
268 if sign_inf != 0 {
269 return if sign_inf > 0 {
270 (float_infinity!(), Equal)
271 } else {
272 (float_negative_infinity!(), Equal)
273 };
274 }
275 // At this point every term is finite.
276 if regulars.is_empty() {
277 // All the terms were zeros (and there is at least one). The dot product is zero with
278 // sign sign_zero.
279 assert_ne!(sign_zero, 0);
280 return if sign_zero > 0 {
281 (float_zero!(), Equal)
282 } else {
283 (float_negative_zero!(), Equal)
284 };
285 }
286 // Optimize the case where there are only one or two regular terms, delegating to the
287 // correctly-rounded multiplication and fused multiply-add-multiply.
288 if regulars.len() == 1 {
289 return regulars[0]
290 .0
291 .mul_prec_round_ref_ref(regulars[0].1, prec, rm);
292 } else if regulars.len() == 2 {
293 return regulars[0].0.mul_add_mul_prec_round_ref_ref_ref_ref(
294 regulars[0].1,
295 regulars[1].0,
296 regulars[1].1,
297 prec,
298 rm,
299 );
300 }
301 let (kernel_rm, exact) = if rm == Exact {
302 (Nearest, true)
303 } else {
304 (rm, false)
305 };
306 // Compute each product exactly at the significand level. A Float's significand is stored
307 // limb-aligned with its top bit set, and its value is significand * 2^(exponent - 64 *
308 // len); the product of two such significands has its top bit either exactly at the combined
309 // width (in which case the product is already aligned) or one position below it (in which
310 // case a shift by 1 restores the alignment and the exponent decreases by 1).
311 let terms: Vec<(bool, i64, u64, Natural)> = regulars
312 .iter()
313 .map(|&(x, y)| {
314 let (sx, ex, px, sig_x) = parts(x);
315 let (sy, ey, py, sig_y) = parts(y);
316 let mut s = sig_x * sig_y;
317 let mut exp = i64::from(ex) + i64::from(ey);
318 let full = sig_x.significant_bits() + sig_y.significant_bits();
319 let mut term_prec = px + py;
320 if s.significant_bits() < full {
321 s <<= 1;
322 exp -= 1;
323 term_prec -= 1;
324 }
325 (sx == sy, exp, term_prec, s)
326 })
327 .collect();
328 let inputs: Vec<FloatSumInput> = terms
329 .iter()
330 .map(|(sign, exp, term_prec, s)| FloatSumInput {
331 sign: *sign,
332 exp: *exp,
333 prec: *term_prec,
334 significand: s,
335 })
336 .collect();
337 complete_sum_result(
338 sum_float_significands(&inputs, prec, kernel_rm),
339 prec,
340 rm,
341 exact,
342 "Inexact Float dot product",
343 )
344 }
345
346 /// Computes the dot product of two equal-length slices of [`Float`]s, rounding the result to
347 /// the nearest value of the specified precision. An [`Ordering`] is also returned, indicating
348 /// whether the rounded dot product is less than, equal to, or greater than the exact dot
349 /// product. Although `NaN`s are not comparable to any [`Float`], whenever this function returns
350 /// a `NaN` it also returns `Equal`.
351 ///
352 /// The products are never rounded, and only a single rounding is performed, at the end: the
353 /// result is the correctly-rounded exact dot product, and intermediate overflow and underflow
354 /// cannot occur.
355 ///
356 /// If the dot product is equidistant from two [`Float`]s with the specified precision, the
357 /// [`Float`] with fewer 1s in its binary expansion is chosen. See [`RoundingMode`] for a
358 /// description of the `Nearest` rounding mode.
359 ///
360 /// $$
361 /// f((x_i)_ {i=0}^{n-1}, (y_i)_ {i=0}^{n-1}, p) = \sum_ {i=0}^{n-1} x_i y_i + \varepsilon.
362 /// $$
363 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored
364 /// or assumed to be 0.
365 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, then $|\varepsilon| \leq
366 /// 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p}$.
367 ///
368 /// If the output has a precision, it is `prec`.
369 ///
370 /// See [`Float::dot_prec_round`] for a description of the special cases and of overflow and
371 /// underflow behavior.
372 ///
373 /// If you know that your target precision is the maximum of the precisions of the inputs,
374 /// consider using [`Float::dot`] instead.
375 ///
376 /// # Worst-case complexity
377 /// $T(n, m, p) = O(n + m (n + p) + m \log m \log\log m)$
378 ///
379 /// $M(n, m, p) = O(n + p + m \log m)$
380 ///
381 /// where $T$ is time, $M$ is additional memory, $n$ is `xs.len()`, $m$ is the sum of the
382 /// significant bits of the elements of `xs` and `ys`, and $p$ is `prec`.
383 ///
384 /// # Panics
385 /// Panics if `prec` is zero or if `xs` and `ys` have different lengths.
386 ///
387 /// # Examples
388 /// ```
389 /// use malachite_base::num::basic::traits::{One, Two};
390 /// use malachite_float::Float;
391 /// use std::cmp::Ordering::*;
392 ///
393 /// let xs = [Float::ONE, Float::TWO, Float::from(3)];
394 /// let ys = [Float::from(4), Float::from(5), Float::from(6)];
395 /// let (dot, o) = Float::dot_prec(&xs, &ys, 10);
396 /// assert_eq!(dot.to_string(), "32.000");
397 /// assert_eq!(o, Equal);
398 ///
399 /// let (dot, o) = Float::dot_prec(&xs, &ys, 3);
400 /// assert_eq!(dot.to_string(), "32.0");
401 /// assert_eq!(o, Equal);
402 /// ```
403 #[inline]
404 pub fn dot_prec(xs: &[Self], ys: &[Self], prec: u64) -> (Self, Ordering) {
405 Self::dot_prec_round(xs, ys, prec, Nearest)
406 }
407
408 /// Computes the dot product of two equal-length slices of [`Float`]s, rounding the result with
409 /// the specified rounding mode. The precision of the result is the maximum of the precisions of
410 /// the inputs (or 1 if there are no inputs). An [`Ordering`] is also returned, indicating
411 /// whether the rounded dot product is less than, equal to, or greater than the exact dot
412 /// product. Although `NaN`s are not comparable to any [`Float`], whenever this function returns
413 /// a `NaN` it also returns `Equal`.
414 ///
415 /// The products are never rounded, and only a single rounding is performed, at the end: the
416 /// result is the correctly-rounded exact dot product, and intermediate overflow and underflow
417 /// cannot occur.
418 ///
419 /// See [`RoundingMode`] for a description of the possible rounding modes.
420 ///
421 /// $$
422 /// f((x_i)_ {i=0}^{n-1}, (y_i)_ {i=0}^{n-1}, m) = \sum_ {i=0}^{n-1} x_i y_i + \varepsilon.
423 /// $$
424 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored
425 /// or assumed to be 0.
426 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, and $m$ is not `Nearest`, then
427 /// $|\varepsilon| < 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p+1}$, where $p$ is
428 /// the maximum precision of the inputs.
429 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, and $m$ is `Nearest`, then
430 /// $|\varepsilon| \leq 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p}$, where $p$ is
431 /// the maximum precision of the inputs.
432 ///
433 /// See [`Float::dot_prec_round`] for a description of the special cases and of overflow and
434 /// underflow behavior.
435 ///
436 /// If you know you'll be using `Nearest`, consider using [`Float::dot`] instead.
437 ///
438 /// # Worst-case complexity
439 /// $T(n, m) = O(m (n + m))$
440 ///
441 /// $M(n, m) = O(n + m \log m)$
442 ///
443 /// where $T$ is time, $M$ is additional memory, $n$ is `xs.len()`, and $m$ is the sum of the
444 /// significant bits of the elements of `xs` and `ys`.
445 ///
446 /// # Panics
447 /// Panics if `xs` and `ys` have different lengths, or if `rm` is `Exact` and the exact dot
448 /// product is not exactly representable with the maximum of the precisions of the inputs.
449 ///
450 /// # Examples
451 /// ```
452 /// use malachite_base::num::arithmetic::traits::PowerOf2;
453 /// use malachite_base::num::basic::traits::Two;
454 /// use malachite_base::rounding_modes::RoundingMode::*;
455 /// use malachite_float::Float;
456 /// use std::cmp::Ordering::*;
457 ///
458 /// // 0.25 * 0.25 + 2 * 5 = 10.0625, whose inputs have maximum precision 3
459 /// let xs = [Float::power_of_2(-2i64), Float::TWO];
460 /// let ys = [Float::power_of_2(-2i64), Float::from(5)];
461 /// let (dot, o) = Float::dot_round(&xs, &ys, Floor);
462 /// assert_eq!(dot.to_string(), "10.0");
463 /// assert_eq!(o, Less);
464 ///
465 /// let (dot, o) = Float::dot_round(&xs, &ys, Ceiling);
466 /// assert_eq!(dot.to_string(), "12.0");
467 /// assert_eq!(o, Greater);
468 /// ```
469 #[inline]
470 pub fn dot_round(xs: &[Self], ys: &[Self], rm: RoundingMode) -> (Self, Ordering) {
471 Self::dot_prec_round(xs, ys, max_prec(xs.iter().chain(ys.iter())), rm)
472 }
473
474 /// Computes the dot product of two equal-length slices of [`Float`]s. The precision of the
475 /// result is the maximum of the precisions of the inputs (or 1 if there are no inputs), and the
476 /// dot product is rounded to nearest.
477 ///
478 /// The products are never rounded, and only a single rounding is performed, at the end: the
479 /// result is the correctly-rounded exact dot product, and intermediate overflow and underflow
480 /// cannot occur.
481 ///
482 /// $$
483 /// f((x_i)_ {i=0}^{n-1}, (y_i)_ {i=0}^{n-1}) = \sum_ {i=0}^{n-1} x_i y_i + \varepsilon.
484 /// $$
485 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is infinite, zero, or `NaN`, $\varepsilon$ may be ignored
486 /// or assumed to be 0.
487 /// - If $\sum_ {i=0}^{n-1} x_i y_i$ is finite and nonzero, then $|\varepsilon| \leq
488 /// 2^{\lfloor\log_2 |\sum_ {i=0}^{n-1} x_i y_i|\rfloor-p}$, where $p$ is the maximum
489 /// precision of the inputs.
490 ///
491 /// See [`Float::dot_prec_round`] for a description of the special cases and of overflow and
492 /// underflow behavior.
493 ///
494 /// # Worst-case complexity
495 /// $T(n, m) = O(m (n + m))$
496 ///
497 /// $M(n, m) = O(n + m \log m)$
498 ///
499 /// where $T$ is time, $M$ is additional memory, $n$ is `xs.len()`, and $m$ is the sum of the
500 /// significant bits of the elements of `xs` and `ys`.
501 ///
502 /// # Panics
503 /// Panics if `xs` and `ys` have different lengths.
504 ///
505 /// # Examples
506 /// ```
507 /// use malachite_base::num::basic::traits::{One, Two};
508 /// use malachite_float::Float;
509 ///
510 /// let xs = [Float::ONE, Float::TWO, Float::from(3)];
511 /// let ys = [Float::from(4), Float::from(5), Float::from(6)];
512 /// assert_eq!(Float::dot(&xs, &ys).to_string(), "32.0");
513 /// ```
514 #[inline]
515 pub fn dot(xs: &[Self], ys: &[Self]) -> Self {
516 Self::dot_round(xs, ys, Nearest).0
517 }
518}