multidouble/
precision.rs

1//! The types and operations of one multi-double precision behind one trait, so that code is written once over
2//! double-double ([`Double`]), quad-double ([`Quad`]) and octo-double ([`Octo`]).
3//!
4//! The precisions have the same kinds of types (scalar, 8-lane vector, forward-mode dual, power series, dense LU
5//! and sum-of-products accumulator) with methods that differ in name or arity: `QD::new(a, b, c, d)` against
6//! `OD::new([..; 8])`, `Fl::load_parts` against `Fl::load_parts_od`, `coef_qd` against `coef_od`, `QDLu::factor` of
7//! four parts against `ODLu::factor` of eight. [`Precision`] names each kind once and delegates to those methods, so a
8//! generic function gives bitwise the results of the code written for one precision.
9//!
10//! ```
11//! use multidouble::precision::{Octo, Precision, Quad};
12//! use multidouble::simd::Native;
13//!
14//! // the same generic code in both precisions: a sum of products, accumulated with one rounding at the end
15//! fn norm2<P: Precision<Native>>(x: &[f64]) -> f64 {
16//!     let v: Vec<P::Scalar> = x.iter().map(|&a| P::new(&[a])).collect();
17//!     P::to_dd(P::dot(v.iter().map(|&a| (a, a)))).hi
18//! }
19//! assert_eq!(norm2::<Quad>(&[3.0, 4.0]), 25.0);
20//! assert_eq!(norm2::<Octo>(&[3.0, 4.0]), 25.0);
21//! ```
22use crate::dd::DD;
23use crate::dual::{Du, Lane, Ou, Qu};
24use crate::dv::DV;
25use crate::field::{Fl, MatLane};
26use crate::jet::{Jet, OJet, QJet};
27use crate::linalg::{DDLu, ODLu, QDLu};
28use crate::od::{OdAcc, OD};
29use crate::ov::OV;
30use crate::qd::{QdAcc, QD};
31use crate::qv::QV;
32use crate::simd::Simd;
33use std::fmt::Debug;
34use std::ops::{Add, Mul, Neg, Sub};
35
36/// One multi-double precision on the backend `S`: its types, and the operations whose names differ between
37/// precisions. Every method delegates to the precision's own method and returns its result bit for bit.
38pub trait Precision<S: Simd>: Copy + Send + Sync + 'static {
39    /// the number of f64 parts of a scalar: 2 for double-double, 4 for quad-double, 8 for octo-double
40    const PARTS: usize;
41    /// the scalar: `DD`, `QD`, `OD`
42    type Scalar: Copy
43        + Send
44        + Sync
45        + Debug
46        + PartialEq
47        + Add<Output = Self::Scalar>
48        + Sub<Output = Self::Scalar>
49        + Mul<Output = Self::Scalar>
50        + Neg<Output = Self::Scalar>;
51    /// 8 scalars in lanes: `DV`, `QV`, `OV`
52    type V: MatLane<S = S>;
53    /// forward-mode dual, value and one tangent in this precision: `Du`, `Qu`, `Ou`
54    type U: MatLane<S = S>;
55    /// truncated power series with N coefficients in this precision: `Jet<N>`, `QJet<N>`, `OJet<N>`
56    type J<const N: usize>: MatLane<S = S>;
57    /// dense LU with partial pivoting in this precision: `DDLu`, `QDLu`, `ODLu`
58    type Lu: Send + Sync;
59
60    /// zero and one
61    const ZERO: Self::Scalar;
62    const ONE: Self::Scalar;
63
64    /// a scalar from up to PARTS parts, the missing ones zero, renormalized (`QD::new`, `OD::new`)
65    ///
66    /// # Panics
67    /// `p` has more than PARTS parts.
68    fn new(p: &[f64]) -> Self::Scalar;
69    /// a scalar from exactly PARTS parts as they are, not renormalized (the struct literal `QD { x }`, `OD { x }`)
70    ///
71    /// # Panics
72    /// `p` does not have exactly PARTS parts.
73    fn from_raw(p: &[f64]) -> Self::Scalar;
74    /// part q of a scalar
75    ///
76    /// # Panics
77    /// `q` is not below PARTS.
78    fn part(x: Self::Scalar, q: usize) -> f64;
79    /// from double-double, exactly
80    fn from_dd(d: DD) -> Self::Scalar;
81    /// rounded to double-double
82    fn to_dd(x: Self::Scalar) -> DD;
83    /// sum of the products a b over the pairs, in this precision's accumulator (`QdAcc`, `OdAcc`): `add_prod` per
84    /// pair, `finish` at the end
85    fn dot(pairs: impl Iterator<Item = (Self::Scalar, Self::Scalar)>) -> Self::Scalar;
86
87    /// the scalar in every lane (`QV::splat`, `OV::splat`)
88    fn splat(x: Self::Scalar) -> Self::V;
89    /// a dual from its value and tangent lanes
90    fn dual(v: Self::V, d: Self::V) -> Self::U;
91    /// coefficient k of a lane in this precision (`Lane::coef_qd`, `Lane::coef_od`)
92    fn coef<L: Lane<S = S>>(l: L, k: usize) -> Self::V;
93    /// set coefficient k of a lane from this precision (`Lane::set_coef_qd`, `Lane::set_coef_od`)
94    fn set_coef<L: Lane<S = S>>(l: &mut L, k: usize, c: Self::V);
95
96    /// a field of 8 B values from PARTS part arrays of at least 8 B values each (`Fl::load_parts`,
97    /// `Fl::load_parts_od`)
98    ///
99    /// # Panics
100    /// `parts` does not hold exactly PARTS arrays, or one is shorter than 8 B.
101    fn load_field<const B: usize>(parts: &[&[f64]]) -> Fl<Self::V, B>;
102    /// value j of a field (`Fl::get`)
103    fn field_get<const B: usize>(f: &Fl<Self::V, B>, j: usize) -> Self::Scalar;
104    /// coefficient k of a field in this precision (`Fl::coef_qd`, `Fl::coef_od`)
105    fn field_coef<L: Lane<S = S>, const B: usize>(f: &Fl<L, B>, k: usize) -> Fl<Self::V, B>;
106    /// set coefficient k of a field from this precision (`Fl::set_coef_qd`, `Fl::set_coef_od`)
107    fn field_set_coef<L: Lane<S = S>, const B: usize>(f: &mut Fl<L, B>, k: usize, c: &Fl<Self::V, B>);
108    /// a dual field from its values and tangents, exactly (`Fl::from_vt_qd`, `Fl::from_vt_od`)
109    fn field_dual<const B: usize>(v: &Fl<Self::V, B>, t: &Fl<Self::V, B>) -> Fl<Self::U, B>;
110
111    /// factor the n x n matrix of PARTS column-major part arrays with leading dimension ld (`QDLu::factor`,
112    /// `ODLu::factor`). Each entry's parts must be normalized, as `new` returns them: the pivot search reads the
113    /// leading parts only.
114    ///
115    /// # Errors
116    /// `c` does not hold PARTS arrays, or a pivot is zero or not finite.
117    ///
118    /// # Panics
119    /// `ld` is not a multiple of 8, `ld < n`, or a part array is shorter than `ld * n`.
120    fn lu_factor(n: usize, ld: usize, c: Vec<Vec<f64>>, par: bool) -> Result<Self::Lu, String>;
121    /// solve in place with the factored matrix (`QDLu::solve`, `ODLu::solve`)
122    fn lu_solve(lu: &Self::Lu, b: &mut [Self::Scalar]);
123}
124
125/// Double-double: `DD`, `DV`, `Du`, `Jet<N>`, `DDLu`; its dot product is the plain sum of double-double products
126#[derive(Clone, Copy, Debug)]
127pub struct Double;
128
129/// Quad-double: `QD`, `QV`, `Qu`, `QJet<N>`, `QDLu`, `QdAcc`
130#[derive(Clone, Copy, Debug)]
131pub struct Quad;
132
133/// Octo-double: `OD`, `OV`, `Ou`, `OJet<N>`, `ODLu`, `OdAcc`
134#[derive(Clone, Copy, Debug)]
135pub struct Octo;
136
137fn parts_of<const P: usize>(p: &[f64], what: &str) -> [f64; P] {
138    assert!(p.len() <= P, "{what}: at most {P} parts, not {}", p.len());
139    std::array::from_fn(|q| if q < p.len() { p[q] } else { 0.0 })
140}
141
142impl<S: Simd> Precision<S> for Double {
143    const PARTS: usize = 2;
144    type Scalar = DD;
145    type V = DV<S>;
146    type U = Du<S>;
147    type J<const N: usize> = Jet<S, N>;
148    type Lu = DDLu;
149    const ZERO: DD = DD::ZERO;
150    const ONE: DD = DD::ONE;
151
152    #[inline]
153    fn new(p: &[f64]) -> DD {
154        let x: [f64; 2] = parts_of(p, "Double::new");
155        DD::new(x[0], x[1])
156    }
157    #[inline]
158    fn from_raw(p: &[f64]) -> DD {
159        assert_eq!(p.len(), 2, "Double::from_raw: 2 parts");
160        DD { hi: p[0], lo: p[1] }
161    }
162    #[inline]
163    fn part(x: DD, q: usize) -> f64 {
164        [x.hi, x.lo][q]
165    }
166    #[inline]
167    fn from_dd(d: DD) -> DD {
168        d
169    }
170    #[inline]
171    fn to_dd(x: DD) -> DD {
172        x
173    }
174    fn dot(pairs: impl Iterator<Item = (DD, DD)>) -> DD {
175        let mut acc = DD::ZERO;
176        for (a, b) in pairs {
177            acc = acc + a * b;
178        }
179        acc
180    }
181    #[inline]
182    fn splat(x: DD) -> DV<S> {
183        DV::splat(x)
184    }
185    #[inline]
186    fn dual(v: DV<S>, d: DV<S>) -> Du<S> {
187        Du { v, d }
188    }
189    #[inline]
190    fn coef<L: Lane<S = S>>(l: L, k: usize) -> DV<S> {
191        l.coef(k)
192    }
193    #[inline]
194    fn set_coef<L: Lane<S = S>>(l: &mut L, k: usize, c: DV<S>) {
195        l.set_coef(k, c)
196    }
197    fn load_field<const B: usize>(parts: &[&[f64]]) -> Fl<DV<S>, B> {
198        assert_eq!(parts.len(), 2, "Double::load_field: 2 parts");
199        Fl::load(parts[0], parts[1])
200    }
201    #[inline]
202    fn field_get<const B: usize>(f: &Fl<DV<S>, B>, j: usize) -> DD {
203        f.get(j)
204    }
205    #[inline]
206    fn field_coef<L: Lane<S = S>, const B: usize>(f: &Fl<L, B>, k: usize) -> Fl<DV<S>, B> {
207        f.coef(k)
208    }
209    #[inline]
210    fn field_set_coef<L: Lane<S = S>, const B: usize>(f: &mut Fl<L, B>, k: usize, c: &Fl<DV<S>, B>) {
211        f.set_coef(k, c)
212    }
213    #[inline]
214    fn field_dual<const B: usize>(v: &Fl<DV<S>, B>, t: &Fl<DV<S>, B>) -> Fl<Du<S>, B> {
215        Fl::<Du<S>, B>::from_vt(v, t)
216    }
217    fn lu_factor(n: usize, ld: usize, c: Vec<Vec<f64>>, par: bool) -> Result<DDLu, String> {
218        let [hi, lo]: [Vec<f64>; 2] = c
219            .try_into()
220            .map_err(|c: Vec<Vec<f64>>| format!("Double::lu_factor: 2 parts, not {}", c.len()))?;
221        DDLu::factor(n, ld, hi, lo, par)
222    }
223    #[inline]
224    fn lu_solve(lu: &DDLu, b: &mut [DD]) {
225        lu.solve(b)
226    }
227}
228
229impl<S: Simd> Precision<S> for Quad {
230    const PARTS: usize = 4;
231    type Scalar = QD;
232    type V = QV<S>;
233    type U = Qu<S>;
234    type J<const N: usize> = QJet<S, N>;
235    type Lu = QDLu;
236    const ZERO: QD = QD::ZERO;
237    const ONE: QD = QD::ONE;
238
239    #[inline]
240    fn new(p: &[f64]) -> QD {
241        let x: [f64; 4] = parts_of(p, "Quad::new");
242        QD::new(x[0], x[1], x[2], x[3])
243    }
244    #[inline]
245    fn from_raw(p: &[f64]) -> QD {
246        assert_eq!(p.len(), 4, "Quad::from_raw: 4 parts");
247        QD {
248            x: [p[0], p[1], p[2], p[3]],
249        }
250    }
251    #[inline]
252    fn part(x: QD, q: usize) -> f64 {
253        x.x[q]
254    }
255    #[inline]
256    fn from_dd(d: DD) -> QD {
257        QD::from_dd(d)
258    }
259    #[inline]
260    fn to_dd(x: QD) -> DD {
261        x.to_dd()
262    }
263    fn dot(pairs: impl Iterator<Item = (QD, QD)>) -> QD {
264        let mut acc = QdAcc::zero();
265        for (a, b) in pairs {
266            acc = acc.add_prod(a, b);
267        }
268        acc.finish()
269    }
270    #[inline]
271    fn splat(x: QD) -> QV<S> {
272        QV::splat(x)
273    }
274    #[inline]
275    fn dual(v: QV<S>, d: QV<S>) -> Qu<S> {
276        Qu { v, d }
277    }
278    #[inline]
279    fn coef<L: Lane<S = S>>(l: L, k: usize) -> QV<S> {
280        l.coef_qd(k)
281    }
282    #[inline]
283    fn set_coef<L: Lane<S = S>>(l: &mut L, k: usize, c: QV<S>) {
284        l.set_coef_qd(k, c)
285    }
286    fn load_field<const B: usize>(parts: &[&[f64]]) -> Fl<QV<S>, B> {
287        assert_eq!(parts.len(), 4, "Quad::load_field: 4 parts");
288        Fl::load_parts(parts[0], parts[1], parts[2], parts[3])
289    }
290    #[inline]
291    fn field_get<const B: usize>(f: &Fl<QV<S>, B>, j: usize) -> QD {
292        f.get(j)
293    }
294    #[inline]
295    fn field_coef<L: Lane<S = S>, const B: usize>(f: &Fl<L, B>, k: usize) -> Fl<QV<S>, B> {
296        f.coef_qd(k)
297    }
298    #[inline]
299    fn field_set_coef<L: Lane<S = S>, const B: usize>(f: &mut Fl<L, B>, k: usize, c: &Fl<QV<S>, B>) {
300        f.set_coef_qd(k, c)
301    }
302    #[inline]
303    fn field_dual<const B: usize>(v: &Fl<QV<S>, B>, t: &Fl<QV<S>, B>) -> Fl<Qu<S>, B> {
304        Fl::<Qu<S>, B>::from_vt_qd(v, t)
305    }
306    fn lu_factor(n: usize, ld: usize, c: Vec<Vec<f64>>, par: bool) -> Result<QDLu, String> {
307        let c: [Vec<f64>; 4] = c
308            .try_into()
309            .map_err(|c: Vec<Vec<f64>>| format!("Quad::lu_factor: 4 parts, not {}", c.len()))?;
310        QDLu::factor(n, ld, c, par)
311    }
312    #[inline]
313    fn lu_solve(lu: &QDLu, b: &mut [QD]) {
314        lu.solve(b)
315    }
316}
317
318impl<S: Simd> Precision<S> for Octo {
319    const PARTS: usize = 8;
320    type Scalar = OD;
321    type V = OV<S>;
322    type U = Ou<S>;
323    type J<const N: usize> = OJet<S, N>;
324    type Lu = ODLu;
325    const ZERO: OD = OD::ZERO;
326    const ONE: OD = OD::ONE;
327
328    #[inline]
329    fn new(p: &[f64]) -> OD {
330        OD::new(parts_of(p, "Octo::new"))
331    }
332    #[inline]
333    fn from_raw(p: &[f64]) -> OD {
334        assert_eq!(p.len(), 8, "Octo::from_raw: 8 parts");
335        OD {
336            x: std::array::from_fn(|q| p[q]),
337        }
338    }
339    #[inline]
340    fn part(x: OD, q: usize) -> f64 {
341        x.x[q]
342    }
343    #[inline]
344    fn from_dd(d: DD) -> OD {
345        OD::from_dd(d)
346    }
347    #[inline]
348    fn to_dd(x: OD) -> DD {
349        x.to_dd()
350    }
351    fn dot(pairs: impl Iterator<Item = (OD, OD)>) -> OD {
352        let mut acc = OdAcc::zero();
353        for (a, b) in pairs {
354            acc = acc.add_prod(a, b);
355        }
356        acc.finish()
357    }
358    #[inline]
359    fn splat(x: OD) -> OV<S> {
360        OV::splat(x)
361    }
362    #[inline]
363    fn dual(v: OV<S>, d: OV<S>) -> Ou<S> {
364        Ou { v, d }
365    }
366    #[inline]
367    fn coef<L: Lane<S = S>>(l: L, k: usize) -> OV<S> {
368        l.coef_od(k)
369    }
370    #[inline]
371    fn set_coef<L: Lane<S = S>>(l: &mut L, k: usize, c: OV<S>) {
372        l.set_coef_od(k, c)
373    }
374    fn load_field<const B: usize>(parts: &[&[f64]]) -> Fl<OV<S>, B> {
375        assert_eq!(parts.len(), 8, "Octo::load_field: 8 parts");
376        Fl::load_parts_od(std::array::from_fn(|q| parts[q]))
377    }
378    #[inline]
379    fn field_get<const B: usize>(f: &Fl<OV<S>, B>, j: usize) -> OD {
380        f.get(j)
381    }
382    #[inline]
383    fn field_coef<L: Lane<S = S>, const B: usize>(f: &Fl<L, B>, k: usize) -> Fl<OV<S>, B> {
384        f.coef_od(k)
385    }
386    #[inline]
387    fn field_set_coef<L: Lane<S = S>, const B: usize>(f: &mut Fl<L, B>, k: usize, c: &Fl<OV<S>, B>) {
388        f.set_coef_od(k, c)
389    }
390    #[inline]
391    fn field_dual<const B: usize>(v: &Fl<OV<S>, B>, t: &Fl<OV<S>, B>) -> Fl<Ou<S>, B> {
392        Fl::<Ou<S>, B>::from_vt_od(v, t)
393    }
394    fn lu_factor(n: usize, ld: usize, c: Vec<Vec<f64>>, par: bool) -> Result<ODLu, String> {
395        let c: [Vec<f64>; 8] = c
396            .try_into()
397            .map_err(|c: Vec<Vec<f64>>| format!("Octo::lu_factor: 8 parts, not {}", c.len()))?;
398        ODLu::factor(n, ld, c, par)
399    }
400    #[inline]
401    fn lu_solve(lu: &ODLu, b: &mut [OD]) {
402        lu.solve(b)
403    }
404}
405
406#[cfg(test)]
407mod tests {
408    //! Each method of `Quad` and `Octo` against the precision's own method on the same inputs, bit for bit. The
409    //! control `mutated` swaps two parts on input and must differ.
410    use super::*;
411    use crate::simd::Native;
412
413    type N = Native;
414
415    fn qd_parts(seed: u64) -> [f64; 4] {
416        let a = 1.0 + (seed as f64) / 7.0;
417        let x = QD::from_dd(DD::from_f64(a)) / QD::from_dd(DD::from_f64(3.0 + seed as f64));
418        x.x
419    }
420
421    fn od_parts(seed: u64) -> [f64; 8] {
422        let a = 1.0 + (seed as f64) / 7.0;
423        let x = OD::from_dd(DD::from_f64(a)) / OD::from_dd(DD::from_f64(3.0 + seed as f64));
424        x.x
425    }
426
427    #[test]
428    fn precision_scalars_match_their_types() {
429        for s in 0..20 {
430            let p = qd_parts(s);
431            assert_eq!(<Quad as Precision<N>>::new(&p), QD::new(p[0], p[1], p[2], p[3]));
432            assert_eq!(<Quad as Precision<N>>::new(&p[..2]), QD::new(p[0], p[1], 0.0, 0.0));
433            assert_eq!(<Quad as Precision<N>>::from_raw(&p).x, p);
434            let o = od_parts(s);
435            assert_eq!(<Octo as Precision<N>>::new(&o), OD::new(o));
436            let mut half = o;
437            half[4..].fill(0.0);
438            assert_eq!(<Octo as Precision<N>>::new(&o[..4]), OD::new(half));
439            assert_eq!(<Octo as Precision<N>>::from_raw(&o).x, o);
440            for q in 0..8 {
441                assert_eq!(<Octo as Precision<N>>::part(OD { x: o }, q), o[q]);
442            }
443        }
444    }
445
446    #[test]
447    fn precision_dot_matches_accumulators() {
448        let xs: Vec<QD> = (0..13).map(|s| QD { x: qd_parts(s) }).collect();
449        let mut acc = QdAcc::zero();
450        for w in xs.windows(2) {
451            acc = acc.add_prod(w[0], w[1]);
452        }
453        assert_eq!(
454            <Quad as Precision<N>>::dot(xs.windows(2).map(|w| (w[0], w[1]))),
455            acc.finish()
456        );
457        let os: Vec<OD> = (0..13).map(|s| OD { x: od_parts(s) }).collect();
458        let mut acc = OdAcc::zero();
459        for w in os.windows(2) {
460            acc = acc.add_prod(w[0], w[1]);
461        }
462        assert_eq!(
463            <Octo as Precision<N>>::dot(os.windows(2).map(|w| (w[0], w[1]))),
464            acc.finish()
465        );
466    }
467
468    /// the field loads, reads, coefficient access and dual seeds against the precision's own methods
469    #[test]
470    fn precision_fields_match_their_types() {
471        const B: usize = 2;
472        let cols: Vec<Vec<f64>> = (0..8)
473            .map(|q| (0..8 * B).map(|i| od_parts(i as u64)[q]).collect())
474            .collect();
475        let refs: Vec<&[f64]> = cols.iter().map(|c| c.as_slice()).collect();
476        let fq = <Quad as Precision<N>>::load_field::<B>(&refs[..4]);
477        let gq = Fl::<QV<N>, B>::load_parts(refs[0], refs[1], refs[2], refs[3]);
478        let fo = <Octo as Precision<N>>::load_field::<B>(&refs);
479        let go = Fl::<OV<N>, B>::load_parts_od(std::array::from_fn(|q| refs[q]));
480        let unit_q = Fl::<QV<N>, B>::unit(5);
481        let unit_o = Fl::<OV<N>, B>::unit(5);
482        let dq = <Quad as Precision<N>>::field_dual(&fq, &unit_q);
483        let doo = <Octo as Precision<N>>::field_dual(&fo, &unit_o);
484        for j in 0..8 * B {
485            assert_eq!(<Quad as Precision<N>>::field_get(&fq, j), gq.get(j));
486            assert_eq!(<Octo as Precision<N>>::field_get(&fo, j), go.get(j));
487            assert_eq!(
488                <Quad as Precision<N>>::field_coef(&dq, 1).get(j),
489                Fl::<Qu<N>, B>::from_vt_qd(&gq, &unit_q).coef_qd(1).get(j)
490            );
491            assert_eq!(<Octo as Precision<N>>::field_coef(&doo, 0).get(j), go.get(j));
492            assert_eq!(
493                <Octo as Precision<N>>::field_coef(&doo, 1).get(j),
494                if j == 5 { OD::ONE } else { OD::ZERO }
495            );
496        }
497        let mut jq = Fl::<QJet<N, 4>, B>::zero();
498        <Quad as Precision<N>>::field_set_coef(&mut jq, 2, &gq);
499        let mut jo = Fl::<OJet<N, 4>, B>::zero();
500        <Octo as Precision<N>>::field_set_coef(&mut jo, 2, &go);
501        for j in 0..8 * B {
502            assert_eq!(<Quad as Precision<N>>::field_coef(&jq, 2).get(j), gq.get(j));
503            assert_eq!(<Octo as Precision<N>>::field_coef(&jo, 2).get(j), go.get(j));
504        }
505    }
506
507    #[test]
508    fn precision_lu_matches_its_type() {
509        let n = 5;
510        let ld = 8; // a multiple of 8, as the LUs need
511        let mat = |p: usize, s: u64| -> Vec<Vec<f64>> {
512            (0..p)
513                .map(|q| {
514                    (0..ld * n)
515                        .map(|i| {
516                            if q == 0 {
517                                ((i * 7 + 3) % 11) as f64 - 5.0 + if i % ld == i / ld { 9.0 } else { 0.0 }
518                            } else {
519                                od_parts(s + i as u64)[q] * 1e-20
520                            }
521                        })
522                        .collect()
523                })
524                .collect()
525        };
526        let cq = mat(4, 1);
527        let mut bq: Vec<QD> = (0..n).map(|i| QD { x: qd_parts(i as u64) }).collect();
528        let mut bq2 = bq.clone();
529        let lu = <Quad as Precision<N>>::lu_factor(n, ld, cq.clone(), false).unwrap();
530        <Quad as Precision<N>>::lu_solve(&lu, &mut bq);
531        QDLu::factor(n, ld, cq.try_into().unwrap(), false)
532            .unwrap()
533            .solve(&mut bq2);
534        assert_eq!(bq, bq2);
535        let co = mat(8, 2);
536        let mut bo: Vec<OD> = (0..n).map(|i| OD { x: od_parts(i as u64) }).collect();
537        let mut bo2 = bo.clone();
538        let lu = <Octo as Precision<N>>::lu_factor(n, ld, co.clone(), false).unwrap();
539        <Octo as Precision<N>>::lu_solve(&lu, &mut bo);
540        ODLu::factor(n, ld, co.try_into().unwrap(), false)
541            .unwrap()
542            .solve(&mut bo2);
543        assert_eq!(bo, bo2);
544        assert!(<Quad as Precision<N>>::lu_factor(n, ld, mat(8, 1), false).is_err());
545    }
546
547    /// Double's methods against the double-double types' own, bit for bit: scalars, the dot product (hp's plain sum of
548    /// products), fields and their dual and jet coefficients, and the LU
549    #[test]
550    fn precision_double_matches_its_types() {
551        const B: usize = 2;
552        type D = Double;
553        for s in 0..20 {
554            let q = qd_parts(s);
555            assert_eq!(<D as Precision<N>>::new(&q[..2]), DD::new(q[0], q[1]));
556            assert_eq!(<D as Precision<N>>::new(&q[..1]), DD::new(q[0], 0.0));
557            let r = <D as Precision<N>>::from_raw(&q[..2]);
558            assert_eq!((r.hi, r.lo), (q[0], q[1]));
559            assert_eq!(
560                (<D as Precision<N>>::part(r, 0), <D as Precision<N>>::part(r, 1)),
561                (q[0], q[1])
562            );
563        }
564        let xs: Vec<DD> = (0..13).map(|s| DD::new(qd_parts(s)[0], qd_parts(s)[1])).collect();
565        let mut acc = DD::ZERO;
566        for w in xs.windows(2) {
567            acc = acc + w[0] * w[1];
568        }
569        assert_eq!(<D as Precision<N>>::dot(xs.windows(2).map(|w| (w[0], w[1]))), acc);
570        let cols: Vec<Vec<f64>> = (0..2)
571            .map(|q| (0..8 * B).map(|i| qd_parts(i as u64)[q]).collect())
572            .collect();
573        let f = <D as Precision<N>>::load_field::<B>(&[&cols[0], &cols[1]]);
574        let g = Fl::<DV<N>, B>::load(&cols[0], &cols[1]);
575        let unit = Fl::<DV<N>, B>::unit(3);
576        let du = <D as Precision<N>>::field_dual(&f, &unit);
577        let mut jet = Fl::<Jet<N, 4>, B>::zero();
578        <D as Precision<N>>::field_set_coef(&mut jet, 2, &g);
579        for j in 0..8 * B {
580            assert_eq!(<D as Precision<N>>::field_get(&f, j), g.get(j));
581            assert_eq!(<D as Precision<N>>::field_coef(&du, 0).get(j), g.get(j));
582            assert_eq!(
583                <D as Precision<N>>::field_coef(&du, 1).get(j),
584                if j == 3 { DD::ONE } else { DD::ZERO }
585            );
586            assert_eq!(<D as Precision<N>>::field_coef(&jet, 2).get(j), g.get(j));
587        }
588        let (n, ld) = (5, 8);
589        let hi: Vec<f64> = (0..ld * n)
590            .map(|i| ((i * 7 + 3) % 11) as f64 - 5.0 + if i % ld == i / ld { 9.0 } else { 0.0 })
591            .collect();
592        let lo: Vec<f64> = (0..ld * n).map(|i| qd_parts(i as u64)[1] * 1e-3).collect();
593        let mut b1: Vec<DD> = (0..n)
594            .map(|i| DD::new(qd_parts(i as u64)[0], qd_parts(i as u64)[1]))
595            .collect();
596        let mut b2 = b1.clone();
597        let lu = <D as Precision<N>>::lu_factor(n, ld, vec![hi.clone(), lo.clone()], false).unwrap();
598        <D as Precision<N>>::lu_solve(&lu, &mut b1);
599        DDLu::factor(n, ld, hi, lo, false).unwrap().solve(&mut b2);
600        assert_eq!(b1, b2);
601        assert!(<D as Precision<N>>::lu_factor(n, ld, vec![vec![0.0; ld * n]; 4], false).is_err());
602    }
603
604    /// the control: the same comparison with two parts swapped on input must fail
605    #[test]
606    fn precision_mutated_differs() {
607        let p = qd_parts(3);
608        let swapped = [p[1], p[0], p[2], p[3]];
609        assert_ne!(<Quad as Precision<N>>::from_raw(&swapped), QD { x: p });
610        let o = od_parts(3);
611        let mut so = o;
612        so.swap(0, 1);
613        assert_ne!(<Octo as Precision<N>>::from_raw(&so), OD { x: o });
614        assert_ne!(
615            <Double as Precision<N>>::from_raw(&[p[1], p[0]]),
616            DD { hi: p[0], lo: p[1] }
617        );
618    }
619}
multidouble/precision.rs

multidouble/
precision.rs
