multidouble 4.3.0

//! Fields (`field::Fl`): Cauchy products, Horner evaluation and constant-matrix products against scalar references,
//! in double-double (to its rounding, against quad-double sums), in the f64 lanes (bitwise the scalar fused
//! multiply-add chains in index order) and in quad-double; the counting lane's matrix product; and the field solves,
//! with the refinement of the double-double solve in the lane type (to quad-double for `QV`), and the order-by-order
//! solves of jets (bit for bit the loops they replace, and against planted solutions formed in exact arithmetic).
//! Each accuracy check has a control fed a wrong result.
#![allow(clippy::needless_range_loop)]
mod common;
use common::Rng;
use multidouble::dd::DD;
use multidouble::dual::{op, take_counts, Counted, Lane};
use multidouble::exact::Dy;
use multidouble::field::{cs, cs_range, cs_w, horner, DMat, FieldSolve, Fl, MatLane, TangentTail};
use multidouble::linalg::{ld_for, DDLu};
use multidouble::native::{Dc, Dtn, Du, Fu, Fun, Jet, QJet, Qu, F8, QV};
use multidouble::qd::QD;
use multidouble::simd::Native;

type DV = multidouble::native::DV;
const B: usize = 2;
const M: usize = 8 * B;

fn field(r: &mut Rng, scale: f64) -> (Fl<DV, B>, Vec<DD>) {
    let v: Vec<DD> = (0..M)
        .map(|_| DD::new(scale * r.sym(), scale * r.sym() * 1e-17))
        .collect();
    let (h, l): (Vec<f64>, Vec<f64>) = v.iter().map(|d| (d.hi, d.lo)).unzip();
    (Fl::load(&h, &l), v)
}

fn q(d: DD) -> QD {
    QD::from_dd(d)
}

/// Cauchy products in double-double against quad-double sums of the same terms: within 2^-100 of sum |x_j y_{k-j}|
/// (w_j |x_j y_{k-j}| for cs_w). Control: the reference with y shifted by one order must miss that bound.
#[test]
fn cauchy_products_to_double_double() {
    let mut r = Rng(51);
    let n = 12;
    let (xs, xv): (Vec<Fl<DV, B>>, Vec<Vec<DD>>) = (0..n).map(|_| field(&mut r, 1.0)).unzip();
    let (ys, yv): (Vec<Fl<DV, B>>, Vec<Vec<DD>>) = (0..n).map(|_| field(&mut r, 1.0)).unzip();
    let w = |j: usize| (j + 1) as f64;
    let (mut worst, mut ctl) = (0.0f64, f64::INFINITY);
    for k in 0..n {
        let (lo, hi) = (k / 3, k - k / 4);
        let got = [
            cs(&xs, &ys, k),
            cs_range(&xs, &ys, k, lo, hi),
            cs_w(&xs, &ys, k, lo, hi, w),
        ];
        for i in 0..M {
            for (case, g) in got.iter().enumerate() {
                let (a, b, wt): (usize, usize, bool) = match case {
                    0 => (0, k, false),
                    1 => (lo, hi, false),
                    _ => (lo, hi, true),
                };
                let (mut s, mut size, mut sc) = (QD::ZERO, 0.0f64, QD::ZERO);
                for j in a..=b {
                    let f = if wt { w(j) } else { 1.0 };
                    let t = q(xv[j][i]).mul_f(f) * q(yv[k - j][i]);
                    s = s + t;
                    size += t.to_f64().abs();
                    if k > j {
                        sc = sc + q(xv[j][i]).mul_f(f) * q(yv[k - j - 1][i]);
                    }
                }
                let gi = q(g.get(i));
                worst = worst.max((gi - s).to_f64().abs() / size);
                if size > 0.0 {
                    ctl = ctl.min((gi - sc).to_f64().abs() / size);
                }
            }
        }
    }
    let bound = 2f64.powi(-100);
    println!(
        "Cauchy products: worst error {:.2e} of sum |terms| (bound {bound:.1e}); control (shifted) {ctl:.1e}",
        worst
    );
    assert!(worst <= bound, "{worst:e}");
    assert!(ctl > bound, "control: the shifted reference agrees ({ctl:e})");
}

/// In the f64 lanes a Cauchy product is the fused multiply-add chain in index order, bitwise (F8, and Fun<2>'s value
/// and tangents); Horner in double-double is the scalar Horner bitwise.
#[test]
fn f64_lanes_and_horner_are_the_scalar_chains() {
    let mut r = Rng(52);
    let n = 9;
    let xs: Vec<Fl<DV, B>> = (0..n).map(|_| field(&mut r, 1.0).0).collect();
    let ys: Vec<Fl<DV, B>> = (0..n).map(|_| field(&mut r, 1.0).0).collect();
    let f8 = |x: &Fl<DV, B>| Fl::<F8, B>::from_vt(x, &Fl::zero());
    let (x8, y8): (Vec<_>, Vec<_>) = (xs.iter().map(f8).collect(), ys.iter().map(f8).collect());
    for k in 0..n {
        let c = cs(&x8, &y8, k);
        let cw = cs_w(&x8, &y8, k, 0, k, |j| (2 * j + 1) as f64);
        for i in 0..M {
            let (mut a, mut aw) = (0.0f64, 0.0f64);
            for j in 0..=k {
                let (xj, yj) = (xs[j].get(i).hi, ys[k - j].get(i).hi);
                a = xj.mul_add(yj, a);
                aw = (xj * (2 * j + 1) as f64).mul_add(yj, aw);
            }
            assert_eq!(
                c.0[i / 8].0.to_array()[i % 8].to_bits(),
                a.to_bits(),
                "F8 cs k {k} i {i}"
            );
            assert_eq!(
                cw.0[i / 8].0.to_array()[i % 8].to_bits(),
                aw.to_bits(),
                "F8 cs_w k {k} i {i}"
            );
        }
    }
    let h = DD::new(0.3, 1e-18);
    let hv = horner(&xs, h);
    for i in 0..M {
        let mut s = DD::ZERO;
        for k in (0..n).rev() {
            s = s * h + xs[k].get(i);
        }
        let g = hv.get(i);
        assert!(
            g.hi.to_bits() == s.hi.to_bits() && g.lo.to_bits() == s.lo.to_bits(),
            "Horner i {i}"
        );
    }
}

fn matrix(r: &mut Rng) -> (DMat<Native, B>, Vec<DD>) {
    let a: Vec<DD> = (0..M * M).map(|_| DD::new(r.sym(), r.sym() * 1e-17)).collect();
    let (h, l): (Vec<f64>, Vec<f64>) = a.iter().map(|d| (d.hi, d.lo)).unzip();
    (DMat::from_rows(&h, &l), a)
}

/// D f for every lane type: double-double (DV, and Du's value and tangent) within 2^-100 of sum |D_ij f_j| of the
/// quad-double product (control: D transposed misses it); quad-double (QV) within 2^-200; the f64 lanes bitwise the
/// scalar fma chains over the columns' hi parts in column order (F8, Fu, Fun<2>, and Dtn<2>'s tangents, whose value is
/// DV's product bitwise); Counted<DV> bitwise DV's, counted as 8 B x B accc and B finish.
#[test]
fn matrix_products() {
    let mut r = Rng(53);
    let (d, a) = matrix(&mut r);
    let (f, fv) = field(&mut r, 1.0);
    let (t, tv) = field(&mut r, 1.0);
    let exact = |i: usize, v: &[DD], tr: bool| -> (QD, f64) {
        let mut s = QD::ZERO;
        let mut size = 0.0f64;
        for j in 0..M {
            let e = if tr { a[j * M + i] } else { a[i * M + j] };
            let p = q(e) * q(v[j]);
            s = s + p;
            size += p.to_f64().abs();
        }
        (s, size)
    };
    let dv = d.apply(&f);
    let du = d.apply(&Fl::<multidouble::native::Du, B>::from_vt(&f, &t));
    let (mut worst, mut ctl) = (0.0f64, f64::INFINITY);
    for i in 0..M {
        let (s, size) = exact(i, &fv, false);
        let (st, sizet) = exact(i, &tv, false);
        let (sx, _) = exact(i, &fv, true);
        worst = worst.max((q(dv.get(i)) - s).to_f64().abs() / size);
        worst = worst.max((q(du.val().get(i)) - s).to_f64().abs() / size);
        worst = worst.max((q(du.tan().get(i)) - st).to_f64().abs() / sizet);
        ctl = ctl.min((q(dv.get(i)) - sx).to_f64().abs() / size);
    }
    println!("D f in double-double: worst {worst:.2e} of sum |terms|; control (transposed) {ctl:.1e}");
    assert!(worst <= 2f64.powi(-100), "{worst:e}");
    assert!(
        ctl > 2f64.powi(-100),
        "control: the transposed product agrees ({ctl:e})"
    );
    // quad-double
    let fq = Fl::<QV, B>::from_vt(&f, &Fl::zero());
    let dq = d.apply(&fq);
    let mut wq = 0.0f64;
    for i in 0..M {
        let (s, size) = exact(i, &fv, false);
        wq = wq.max((dq.0[i / 8].lane(i % 8) - s).to_f64().abs() / size);
    }
    println!("D f in quad-double: worst {wq:.2e} of sum |terms|");
    assert!(wq <= 2f64.powi(-200), "{wq:e}");
    // the f64 lanes: scalar fma chains over the columns' hi parts
    let chain = |v: &[f64], i: usize| -> f64 {
        let mut s = 0.0f64;
        for j in 0..M {
            s = v[j].mul_add(a[i * M + j].hi, s);
        }
        s
    };
    let (fh, th): (Vec<f64>, Vec<f64>) = (fv.iter().map(|x| x.hi).collect(), tv.iter().map(|x| x.hi).collect());
    let lane = |v: multidouble::native::V8, i: usize| v.to_array()[i % 8].to_bits();
    let r8 = d.apply(&Fl::<F8, B>::from_vt(&f, &Fl::zero()));
    let ru = d.apply(&Fl::<Fu, B>::from_vt(&f, &t));
    let mut fun = Fl::<Fun<2>, B>::from_vt(&f, &Fl::zero());
    let mut dtn = Fl::<Dtn<2>, B>::from_vt(&f, &Fl::zero());
    for b in 0..B {
        for k in 0..2 {
            let tk = if k == 0 {
                t.0[b].hi
            } else {
                f.0[b].hi.mul(multidouble::native::V8::splat(0.5))
            };
            fun.0[b].set_tan64(k, tk);
            dtn.0[b].set_tan64(k, tk);
        }
    }
    let (rn, rt) = (d.apply(&fun), d.apply(&dtn));
    for i in 0..M {
        let (cf, ct) = (chain(&fh, i).to_bits(), chain(&th, i).to_bits());
        let half: Vec<f64> = fh.iter().map(|x| x * 0.5).collect();
        let ch = chain(&half, i).to_bits();
        assert_eq!(lane(r8.0[i / 8].0, i), cf, "F8 i {i}");
        assert_eq!((lane(ru.0[i / 8].v, i), lane(ru.0[i / 8].d, i)), (cf, ct), "Fu i {i}");
        assert_eq!(
            (
                lane(rn.0[i / 8].v, i),
                lane(rn.0[i / 8].d[0], i),
                lane(rn.0[i / 8].d[1], i)
            ),
            (cf, ct, ch),
            "Fun i {i}"
        );
        assert_eq!(
            (lane(rt.0[i / 8].d[0], i), lane(rt.0[i / 8].d[1], i)),
            (ct, ch),
            "Dtn tangents i {i}"
        );
        let (v, w) = (rt.0[i / 8].v.lane(i % 8), dv.get(i));
        assert!(
            v.hi.to_bits() == w.hi.to_bits() && v.lo.to_bits() == w.lo.to_bits(),
            "Dtn value i {i}"
        );
    }
    // the counting lane
    let fc = Fl(f.0.map(Counted));
    take_counts();
    let dc = d.apply(&fc);
    let n = take_counts();
    for i in 0..M {
        let (x, y) = (dc.0[i / 8].0.lane(i % 8), dv.get(i));
        assert!(
            x.hi.to_bits() == y.hi.to_bits() && x.lo.to_bits() == y.lo.to_bits(),
            "Counted i {i}"
        );
    }
    // the default product (apply_dd) run through Counted's own per-operation counts must count the same
    take_counts();
    let dd = d.apply_dd(&fc);
    let nd = take_counts();
    assert_eq!(
        n, nd,
        "Counted::dmat_apply against the default product on counting lanes"
    );
    for i in 0..M {
        let (x, y) = (dd.0[i / 8].0.lane(i % 8), dv.get(i));
        assert!(
            x.hi.to_bits() == y.hi.to_bits() && x.lo.to_bits() == y.lo.to_bits(),
            "apply_dd i {i}"
        );
    }
    assert_eq!(n[op::ACCC], (8 * B * B) as u64);
    assert_eq!(n[op::FIN], B as u64);
    assert_eq!(n.iter().sum::<u64>(), (8 * B * B + B) as u64);
    let _ = <Counted<DV> as MatLane>::dmat_apply::<B>;
}

/// The field solves are DDLu's solves of the stored vector, bitwise
#[test]
fn field_solves_are_the_lu_solves() {
    let mut r = Rng(54);
    let ld = ld_for(M);
    let (mut h, mut l) = (vec![0.0; ld * M], vec![0.0; ld * M]);
    for j in 0..M {
        for i in 0..M {
            let v = DD::new(r.sym() + if i == j { 4.0 } else { 0.0 }, r.sym() * 1e-17);
            h[j * ld + i] = v.hi;
            l[j * ld + i] = v.lo;
        }
    }
    let lu = DDLu::factor(M, ld, h, l, false).unwrap();
    let (f, fv) = field(&mut r, 1.0);
    let x = lu.solve_lanes(&f);
    let mut b = fv.clone();
    lu.solve(&mut b);
    let mut bh: Vec<f64> = fv.iter().map(|v| v.hi).collect();
    lu.solve_hi(&mut bh);
    let xh = lu.solve_lanes_hi(&f);
    for i in 0..M {
        let g = x.get(i);
        assert!(
            g.hi.to_bits() == b[i].hi.to_bits() && g.lo.to_bits() == b[i].lo.to_bits(),
            "solve_lanes i {i}"
        );
        assert!(
            xh.get(i).hi.to_bits() == bh[i].to_bits() && xh.get(i).lo == 0.0,
            "solve_lanes_hi i {i}"
        );
    }
}

/// A DDLu whose order is not 8 B is refused by both field solves (smaller: the rest of the field would pass through
/// unsolved; larger: an index past the end)
fn lu_of_order(n: usize) -> DDLu {
    let ld = ld_for(n);
    let mut h = vec![0.0; ld * n];
    for i in 0..n {
        h[i * ld + i] = 2.0;
    }
    DDLu::factor(n, ld, h, vec![0.0; ld * n], false).unwrap()
}

#[test]
#[should_panic(expected = "FieldSolve: a DDLu of order 8 for a field of 16 values")]
fn field_solve_refuses_a_smaller_lu() {
    let (f, _) = field(&mut Rng(3), 1.0);
    lu_of_order(M - 8).solve_lanes(&f);
}

#[test]
#[should_panic(expected = "FieldSolve: a DDLu of order 24 for a field of 16 values")]
fn field_solve_hi_refuses_a_larger_lu() {
    let (f, _) = field(&mut Rng(3), 1.0);
    lu_of_order(M + 8).solve_lanes_hi(&f);
}

/// a quad-double field from 8 B values
fn qfield(v: &[QD]) -> Fl<QV, B> {
    Fl(std::array::from_fn(|b| {
        QV::from_lanes(v[8 * b..8 * b + 8].try_into().unwrap())
    }))
}

/// the 8 B values of a quad-double field
fn qvalues(f: &Fl<QV, B>) -> Vec<QD> {
    f.0.iter().flat_map(|v| v.lanes()).collect()
}

/// A x0 = b: A a diagonally dominant double-double matrix (as the constant matrix for the quad-double residual and
/// as its LU for the solves), x0 a random quad-double field (four parts), b = A x0 formed exactly (dyadic arithmetic)
/// and rounded to quad-double
fn qd_system(r: &mut Rng) -> (DMat<Native, B>, DDLu, Vec<QD>, Vec<QD>) {
    let a: Vec<DD> = (0..M * M)
        .map(|k| DD::new(r.sym() + if k / M == k % M { 4.0 } else { 0.0 }, r.sym() * 1e-17))
        .collect();
    let (h, l): (Vec<f64>, Vec<f64>) = a.iter().map(|d| (d.hi, d.lo)).unzip();
    let ld = ld_for(M);
    let (mut ch, mut cl) = (vec![0.0; ld * M], vec![0.0; ld * M]);
    for i in 0..M {
        for j in 0..M {
            ch[j * ld + i] = h[i * M + j];
            cl[j * ld + i] = l[i * M + j];
        }
    }
    let lu = DDLu::factor(M, ld, ch, cl, false).unwrap();
    let x0: Vec<QD> = (0..M)
        .map(|_| {
            let v = r.sym();
            QD::new(
                v,
                v * 2f64.powi(-53) * r.sym(),
                v * 2f64.powi(-106) * r.sym(),
                v * 2f64.powi(-159) * r.sym(),
            )
        })
        .collect();
    let b: Vec<QD> = (0..M)
        .map(|i| {
            let mut s = Dy::zero();
            for j in 0..M {
                let aij = Dy::from_f64(h[i * M + j]).add(&Dy::from_f64(l[i * M + j]));
                s = s.add(&aij.mul(&x0[j].exact()));
            }
            QD { x: s.greedy4() }
        })
        .collect();
    (DMat::from_rows(&h, &l), lu, x0, b)
}

/// max |x_i - y_i| / max |y_i|, in quad-double
fn qd_err(x: &[QD], y: &[QD]) -> f64 {
    let ym = y.iter().map(|v| v.x[0].abs()).fold(0.0, f64::max);
    x.iter()
        .zip(y)
        .map(|(a, b)| (*a - *b).to_f64().abs())
        .fold(0.0, f64::max)
        / ym
}

/// The field refinement in quad-double lanes is DDLu::solve_refine_qd on the field's values, bitwise (both with the
/// same residual function), for 0 to 3 steps; with no steps it is the double-double field solve of the values, in
/// QV and in DV lanes alike.
#[test]
fn field_refine_is_solve_refine_qd() {
    let mut r = Rng(61);
    let (d, lu, _, b) = qd_system(&mut r);
    let bq = qfield(&b);
    let resid = |x: &Fl<QV, B>| bq - d.apply(x);
    for steps in 0..4 {
        let x = qvalues(&lu.solve_lanes_refine(&bq, steps, resid));
        let y = lu.solve_refine_qd(&b, steps, |x, out| out.copy_from_slice(&qvalues(&resid(&qfield(x)))));
        for i in 0..M {
            assert!(
                x[i].x
                    .iter()
                    .zip(y[i].x.iter())
                    .all(|(p, q)| p.to_bits() == q.to_bits()),
                "steps {steps} i {i}: {:?} vs {:?}",
                x[i],
                y[i]
            );
        }
    }
    let plain = lu.solve_lanes(&bq.val());
    let x0 = lu.solve_lanes_refine(&bq, 0, |_| unreachable!());
    let (f, _) = field(&mut r, 1.0);
    let xd = lu.solve_lanes_refine(&f, 0, |_| unreachable!());
    let fd = lu.solve_lanes(&f);
    for i in 0..M {
        let (p, q) = (x0.0[i / 8].lane(i % 8), plain.get(i));
        assert!(p.x[0].to_bits() == q.hi.to_bits() && p.x[1].to_bits() == q.lo.to_bits() && p.x[2] == 0.0);
        assert!(xd.get(i).hi.to_bits() == fd.get(i).hi.to_bits() && xd.get(i).lo.to_bits() == fd.get(i).lo.to_bits());
    }
}

/// A x = b to quad-double by the field refinement: two steps give x0 back within 1e-60 of max |x0| (about 1e-64 seen;
/// b's rounding to quad-double and A's condition number allow about 1e-62), and three the same. Controls: without
/// refinement the solve misses by more than 1e-40 (double-double accuracy, about 1e-32 seen); with a residual off by
/// 1e-40 relative (A scaled by 1 + 1e-40 in it) the refinement converges to the wrong system and misses 1e-60.
#[test]
fn field_refine_reaches_quad_double() {
    let mut r = Rng(62);
    for _ in 0..4 {
        let (d, lu, x0, b) = qd_system(&mut r);
        let bq = qfield(&b);
        let resid = |x: &Fl<QV, B>| bq - d.apply(x);
        let e2 = qd_err(&qvalues(&lu.solve_lanes_refine(&bq, 2, resid)), &x0);
        let e3 = qd_err(&qvalues(&lu.solve_lanes_refine(&bq, 3, resid)), &x0);
        assert!(
            e2 <= 1e-60 && e3 <= 1e-60,
            "refined: {e2:e} (2 steps), {e3:e} (3 steps)"
        );
        let e0 = qd_err(&qvalues(&lu.solve_lanes_refine(&bq, 0, resid)), &x0);
        assert!(e0 > 1e-40, "control: the unrefined solve is within 1e-40 ({e0:e})");
        let s = QV::splat(QD::from_f64(1.0).add_f(1e-40));
        let wrong = |x: &Fl<QV, B>| bq - Fl(d.apply(x).0.map(|v| v * s));
        let ew = qd_err(&qvalues(&lu.solve_lanes_refine(&bq, 3, wrong)), &x0);
        assert!(ew > 1e-60, "control: a residual off by 1e-40 passes ({ew:e})");
    }
}

#[test]
#[should_panic(expected = "FieldSolve: a DDLu of order 8 for a field of 16 values")]
fn field_refine_refuses_a_smaller_lu() {
    let f = qfield(&[QD::from_f64(1.0); M]);
    lu_of_order(M - 8).solve_lanes_refine(&f, 2, |x| *x);
}

/// A quad-double field's four parts' arrays go in and out unchanged (`Fl::<QV, B>::load_parts` and `store_parts`: not
/// renormalized, so parts that overlap stay as they are), value j at lane j % 8 of vector j / 8 as `QV::from_lanes`
/// places them; `get` reads them back, and `Fl::load_parts` with its type inferred from its use compiles (one
/// `load_parts` among the field types). Control: the test parts overlap, so a load that renormalized would move them.
#[test]
fn qd_field_load_store_get() {
    let mut r = Rng(63);
    // overlapping parts (each about 2^-20 of the one before, far above half its ulp)
    let c: [Vec<f64>; 4] = std::array::from_fn(|k| (0..M).map(|_| r.sym() * 2f64.powi(-20 * k as i32)).collect());
    let f = Fl::<QV, B>::load_parts(&c[0], &c[1], &c[2], &c[3]);
    let mut out: [Vec<f64>; 4] = std::array::from_fn(|_| vec![0.0; M]);
    let [o0, o1, o2, o3] = &mut out;
    f.store_parts(o0, o1, o2, o3);
    for k in 0..4 {
        assert!(
            c[k].iter().zip(&out[k]).all(|(a, b)| a.to_bits() == b.to_bits()),
            "part {k}"
        );
    }
    let lanes: Vec<QD> = (0..M)
        .map(|j| QD {
            x: std::array::from_fn(|k| c[k][j]),
        })
        .collect();
    let g = qfield(&lanes);
    for j in 0..M {
        let (q, v) = (f.get(j), g.0[j / 8].lane(j % 8));
        assert!(
            (0..4).all(|k| q.x[k].to_bits() == c[k][j].to_bits() && v.x[k].to_bits() == c[k][j].to_bits()),
            "value {j}"
        );
    }
    let renorm = QD::new(c[0][0], c[1][0], c[2][0], c[3][0]);
    assert!(
        (0..4).any(|k| renorm.x[k].to_bits() != c[k][0].to_bits()),
        "control: the parts do not overlap"
    );
    // the type inferred from the use, as hp's qsolve calls it: one load_parts on a field type only
    let inferred: Fl<QV, B> = Fl::load_parts(&c[0], &c[1], &c[2], &c[3]);
    assert!((0..M).all(|j| inferred.get(j).x == f.get(j).x), "the inferred load");
}

/// The refinement is for values only: a lane type with tangents (Du here) is refused, not given a zero derivative
#[test]
#[should_panic(expected = "FieldSolve::solve_lanes_refine refines values only")]
fn field_refine_refuses_tangent_lanes() {
    let lu = lu_of_order(M);
    let rhs = Fl::<multidouble::native::Du, B>([multidouble::native::Du::splat(DD::from_f64(1.0)); B]);
    lu.solve_lanes_refine(&rhs, 1, |x| rhs - *x);
}

/// A jet is a lane type with tangents, so the refinement refuses it as well
#[test]
#[should_panic(expected = "FieldSolve::solve_lanes_refine refines values only")]
fn field_refine_refuses_a_jet() {
    let lu = lu_of_order(M);
    let rhs = Fl::<QJet<3>, B>([QJet::<3>::splat(DD::ONE); B]);
    lu.solve_lanes_refine(&rhs, 1, |x| rhs - *x);
}

/// The refinement with a tangent refuses a double-double jet, whose orders past 1 it would return as zero
/// (`Lane::from_vt` sets orders 0 and 1 only)
#[test]
#[should_panic(expected = "FieldSolve::solve_lanes_refine_t refines one tangent: a jet would lose its orders past 1")]
fn field_refine_t_refuses_a_jet() {
    let lu = lu_of_order(M);
    let rhs = Fl::<Jet<3>, B>([Jet::<3>::splat(DD::ONE); B]);
    lu.solve_lanes_refine_t(&lu, &rhs, 1, |x| rhs - *x);
}

/// and a quad-double jet likewise
#[test]
#[should_panic(expected = "FieldSolve::solve_lanes_refine_t refines one tangent: a jet would lose its orders past 1")]
fn field_refine_t_refuses_a_quad_double_jet() {
    let lu = lu_of_order(M);
    let rhs = Fl::<QJet<3>, B>([QJet::<3>::splat(DD::ONE); B]);
    lu.solve_lanes_refine_t(&lu, &rhs, 2, |x| rhs - *x);
}

// ---------------------------------------------------------------------------------------------------------------
// tan64, the complex-tangent fields, the refinement with a tangent, and the early end of a series
// ---------------------------------------------------------------------------------------------------------------

/// `Fl::tan64(k)` is lane by lane the lanes' `tan64(k)`, bit for bit, as an `F8` field
#[test]
fn field_tan64_is_the_lanes_tangent() {
    let mut r = Rng(71);
    let (v, _) = field(&mut r, 1.0);
    let mut f = Fl::<Fun<2>, B>::from_vt(&v, &Fl::zero());
    for b in 0..B {
        for k in 0..2 {
            let a: [f64; 8] = std::array::from_fn(|_| r.sym());
            f.0[b].set_tan64(k, multidouble::v8::V8::load(&a));
        }
    }
    for k in 0..2 {
        let t: Fl<F8, B> = f.tan64(k);
        for b in 0..B {
            assert_eq!(
                t.0[b].0.to_array().map(f64::to_bits),
                f.0[b].tan64(k).to_array().map(f64::to_bits)
            );
        }
    }
    assert_ne!(
        f.tan64(0).0[0].0.to_array(),
        f.tan64(1).0[0].0.to_array(),
        "control: the two tangents agree"
    );
}

/// the bits of a field of any lane type: coefficients 0 and 1 (value and tangent) and a jet's others in quad-double
/// (`Lane::coef_qd`, exact), and the f64 tangents
fn lane_bits<T: Lane<S = Native>>(f: &Fl<T, B>) -> Vec<u64> {
    let mut w: Vec<u64> = (0..T::NJ.max(2))
        .flat_map(|k| qvalues(&f.coef_qd(k)).into_iter().flat_map(|q| q.x.map(f64::to_bits)))
        .collect();
    for k in 0..T::NT {
        w.extend(f.tan64(k).0.iter().flat_map(|t| t.0.to_array().map(f64::to_bits)));
    }
    w
}

/// an f64 field of 8 B random values
fn f8_field(r: &mut Rng) -> Fl<F8, B> {
    Fl(std::array::from_fn(|_| {
        multidouble::dual::F8(multidouble::v8::V8::from_array(std::array::from_fn(|_| r.sym())))
    }))
}

/// `Fl::set_tan64` is the lane-by-lane `Lane::set_tan64` loop it replaces (hp's psolve writes a solved f64 tangent
/// so, `out.0[b].set_tan64(k, fd.0[b].hi)`) bit for bit, on `Fu`, `Dtn<3>` and `Fun<2>`, and `tan64` reads back what
/// it set, with the value and the other tangents unchanged. Control: the loop on the next tangent is caught.
#[test]
fn field_set_tan64_is_the_lane_loop() {
    fn check<T: Lane<S = Native>>(r: &mut Rng) {
        let (v, _) = field(r, 1.0);
        let (d, _) = field(r, 1.0);
        let mut f = Fl::<T, B>::from_vt(&v, &d);
        for k in 0..T::NT {
            let t = f8_field(r);
            let (before, mut g) = (f, f);
            for b in 0..B {
                g.0[b].set_tan64(k, t.0[b].0);
            }
            f.set_tan64(k, &t);
            assert_eq!(lane_bits(&f), lane_bits(&g), "the lane loop, tangent {k}");
            let back: Vec<u64> = f
                .tan64(k)
                .0
                .iter()
                .flat_map(|x| x.0.to_array().map(f64::to_bits))
                .collect();
            let want: Vec<u64> = t.0.iter().flat_map(|x| x.0.to_array().map(f64::to_bits)).collect();
            assert_eq!(back, want, "tan64 after set_tan64, tangent {k}");
            for o in (0..T::NT).filter(|&o| o != k) {
                assert_eq!(lane_bits(&f.tan64(o)), lane_bits(&before.tan64(o)), "tangent {o} moved");
            }
            assert_eq!(lane_bits(&f.val()), lane_bits(&before.val()), "the value moved");
            if k + 1 < T::NT {
                let mut wrong = before;
                for b in 0..B {
                    wrong.0[b].set_tan64(k + 1, t.0[b].0);
                }
                assert_ne!(lane_bits(&f), lane_bits(&wrong), "control: the next tangent");
            }
        }
    }
    let mut r = Rng(75);
    check::<Fu>(&mut r);
    check::<Dtn<3>>(&mut r);
    check::<Fun<2>>(&mut r);
}

#[test]
#[should_panic(expected = "Fl::set_tan64: f64 tangent 0 of a lane type with 0 f64 tangents")]
fn field_set_tan64_refuses_a_double_double_tangent() {
    let mut f = Fl::<Du, B>::zero();
    f.set_tan64(0, &Fl::zero());
}

/// `Fu`'s own `set_tan64` ignores k, so the field's check is what refuses its tangent 1
#[test]
#[should_panic(expected = "Fl::set_tan64: f64 tangent 1 of a lane type with 1 f64 tangents")]
fn field_set_tan64_refuses_a_tangent_past_the_last() {
    let mut f = Fl::<Fu, B>::zero();
    f.set_tan64(1, &Fl::zero());
}

/// The unit seeds are hp's hand-built ones bit for bit: `unit_tangent(v, j)` is `from_vt(v, t)` with t loaded with a
/// one at value j (hp's model::jac_column), on every lane type with a tangent; `set_tan64(k, &unit(j))` on a field
/// from `from_vt(v, 0)` is the lane vector's `set_tan64` with e_{j % 8} (hp's jac_columns_n), for three tangents at
/// once; `from_vt_qd(v, &unit(j))` is the `Qu` field built from the quad-double values with the one at value j
/// (hp's qsolve::jac_column_q), and gives its values and tangents back (`coef_qd`). `unit(j)` reads one at j and zero
/// elsewhere. Control: a seed one value off is caught.
#[test]
fn field_unit_seeds_are_the_hand_built_ones() {
    fn seed<T: Lane<S = Native>>(v: &Fl<DV, B>, j: usize) -> Fl<T, B> {
        let mut th = vec![0.0; M];
        let tl = vec![0.0; M];
        th[j] = 1.0;
        Fl::<T, B>::from_vt(v, &Fl::load(&th, &tl))
    }
    fn check<T: Lane<S = Native>>(v: &Fl<DV, B>) {
        for j in [0, 5, 8, M - 1] {
            assert_eq!(
                lane_bits(&Fl::<T, B>::unit_tangent(v, j)),
                lane_bits(&seed::<T>(v, j)),
                "value {j}"
            );
        }
        assert_ne!(
            lane_bits(&Fl::<T, B>::unit_tangent(v, 3)),
            lane_bits(&seed::<T>(v, 4)),
            "control"
        );
    }
    let mut r = Rng(76);
    let (v, _) = field(&mut r, 1.0);
    check::<Du>(&v);
    check::<Dc>(&v);
    check::<Qu>(&v);
    check::<Fu>(&v);
    check::<Dtn<2>>(&v);
    check::<Fun<3>>(&v);
    check::<Jet<4>>(&v);
    check::<QJet<4>>(&v);
    assert!((0..M).all(|i| Fl::<Dc, B>::unit_tangent(&v, 5).tan_im().get(i) == DD::ZERO));
    for j in 0..M {
        let (e, q) = (Fl::<DV, B>::unit(j), Fl::<QV, B>::unit(j));
        for i in 0..M {
            let one: f64 = if i == j { 1.0 } else { 0.0 };
            assert!(e.get(i).hi.to_bits() == one.to_bits() && e.get(i).lo.to_bits() == 0);
            assert!(q.get(i).x[0].to_bits() == one.to_bits() && q.get(i).x[1..].iter().all(|p| p.to_bits() == 0));
        }
    }
    // three f64 tangents seeded on three values
    let tasks = [(0, 2), (1, 9), (2, M - 1)];
    let mut f = Fl::<Dtn<3>, B>::from_vt(&v, &Fl::zero());
    let mut g = f;
    for &(k, j) in &tasks {
        let mut e = [0.0; 8];
        e[j % 8] = 1.0;
        g.0[j / 8].set_tan64(k, multidouble::v8::V8::from_array(e));
        f.set_tan64(k, &Fl::<F8, B>::unit(j));
    }
    assert_eq!(lane_bits(&f), lane_bits(&g));
    let mut wrong = Fl::<Dtn<3>, B>::from_vt(&v, &Fl::zero());
    wrong.set_tan64(1, &Fl::<F8, B>::unit(10));
    assert_ne!(lane_bits(&f), lane_bits(&wrong), "control");
    // a quad-double Jacobian seed
    let vq: Vec<QD> = (0..M).map(|_| jet_coef(&mut r, 0.0, 1.0, true)).collect();
    let tq: Vec<QD> = (0..M).map(|_| jet_coef(&mut r, 0.0, 1.0, true)).collect();
    let (fv, ft) = (qfield(&vq), qfield(&tq));
    let both = Fl::<Qu, B>::from_vt_qd(&fv, &ft);
    assert!(same_qd(&qvalues(&both.coef_qd(0)), &vq) && same_qd(&qvalues(&both.coef_qd(1)), &tq));
    assert!(!same_qd(&qvalues(&both.coef_qd(1)), &vq), "control");
    for j in [0, 9] {
        let mut hand = [Qu::zero(); B];
        for bb in 0..B {
            hand[bb].v = fv.0[bb];
        }
        let mut e = [QD::ZERO; 8];
        e[j % 8] = QD::ONE;
        hand[j / 8].d = QV::from_lanes(&e);
        let s = Fl::<Qu, B>::from_vt_qd(&fv, &Fl::<QV, B>::unit(j));
        assert_eq!(lane_bits(&s), lane_bits(&Fl(hand)), "Qu seed on value {j}");
    }
}

#[test]
#[should_panic(expected = "Fl::unit_tangent: a lane type without a tangent")]
fn field_unit_tangent_refuses_a_value_lane() {
    Fl::<QV, B>::unit_tangent(&Fl::zero(), 0);
}

#[test]
#[should_panic(expected = "Fl::unit: value 16 of a field of 16 values")]
fn field_unit_refuses_a_value_past_the_field() {
    Fl::<DV, B>::unit(M);
}

/// a `Dc` field from its three parts gives them back bit for bit, and `from_vt` sets the imaginary part to zero
#[test]
fn dc_field_parts_round_trip() {
    let mut r = Rng(72);
    let (v, _) = field(&mut r, 1.0);
    let (re, _) = field(&mut r, 1.0);
    let (im, _) = field(&mut r, 1.0);
    let f = Fl::<Dc, B>::from_vri(&v, &re, &im);
    let same = |a: &Fl<DV, B>, b: &Fl<DV, B>| {
        (0..M).all(|i| a.get(i).hi.to_bits() == b.get(i).hi.to_bits() && a.get(i).lo.to_bits() == b.get(i).lo.to_bits())
    };
    assert!(same(&f.val(), &v) && same(&f.tan_re(), &re) && same(&f.tan_im(), &im) && same(&f.tan(), &re));
    assert!(!same(&f.tan_im(), &re), "control");
    let g = Fl::<Dc, B>::from_vt(&v, &re);
    assert!(same(&g.tan_re(), &re) && same(&g.tan_im(), &Fl::zero()));
}

fn qu_field(v: &Fl<QV, B>, d: &Fl<QV, B>) -> Fl<Qu, B> {
    Fl(std::array::from_fn(|b| Qu { v: v.0[b], d: d.0[b] }))
}

fn qu_val(f: &Fl<Qu, B>) -> Vec<QD> {
    qvalues(&Fl(f.0.map(|q| q.v)))
}

fn qu_tan(f: &Fl<Qu, B>) -> Vec<QD> {
    qvalues(&Fl(f.0.map(|q| q.d)))
}

fn same_qd(x: &[QD], y: &[QD]) -> bool {
    x.iter()
        .zip(y)
        .all(|(a, b)| a.x.iter().zip(b.x.iter()).all(|(p, q)| p.to_bits() == q.to_bits()))
}

/// With the tangent system the value system's (A_d = 0, lu_t the same LU), `solve_lanes_refine_t` on `Qu` is
/// `solve_lanes_refine` on the value and on the tangent separately, bit for bit, for 0 to 3 steps; without a tangent
/// (`QV`) it is `solve_lanes_refine` bit for bit. Control: the tangent against the value's refinement must differ.
#[test]
fn field_refine_t_is_refine_on_value_and_tangent() {
    let mut r = Rng(73);
    let (d, lu, x0, b) = qd_system(&mut r);
    let (bq, tq) = (qfield(&b), qfield(&x0));
    let rhs = qu_field(&bq, &tq);
    for steps in 0..4 {
        let x = lu.solve_lanes_refine_t(&lu, &rhs, steps, |x: &Fl<Qu, B>| rhs - d.apply(x));
        let xv = lu.solve_lanes_refine(&bq, steps, |x: &Fl<QV, B>| bq - d.apply(x));
        let xt = lu.solve_lanes_refine(&tq, steps, |x: &Fl<QV, B>| tq - d.apply(x));
        assert!(same_qd(&qu_val(&x), &qvalues(&xv)), "value, {steps} steps");
        assert!(same_qd(&qu_tan(&x), &qvalues(&xt)), "tangent, {steps} steps");
        assert!(!same_qd(&qu_tan(&x), &qvalues(&xv)), "control, {steps} steps");
        let y = lu.solve_lanes_refine_t(&lu, &bq, steps, |x: &Fl<QV, B>| bq - d.apply(x));
        assert!(same_qd(&qvalues(&y), &qvalues(&xv)), "QV, {steps} steps");
    }
}

/// A(t) x(t) = b with A(t) = A + t E: the tangent of x is -A^-1 E x. With `Qu` lanes whose residual carries E x_v
/// on the tangent, three refinement steps give the tangent to 1e-58 of the reference (two quad-double refinements:
/// x, then -A^-1 E x), the value to 1e-60 as `solve_lanes_refine` does. Controls: one step leaves the tangent above
/// 1e-40 (the first tangent solve ignores E); a residual with E scaled by 1 + 1e-30 converges to the wrong tangent.
#[test]
fn field_refine_t_reaches_quad_double_with_a_matrix_tangent() {
    let mut r = Rng(74);
    let (d, lu, _, b) = qd_system(&mut r);
    let (e, _) = matrix(&mut r);
    let bq = qfield(&b);
    let xv = lu.solve_lanes_refine(&bq, 3, |x: &Fl<QV, B>| bq - d.apply(x));
    let ex = -e.apply(&xv);
    let t_ref = lu.solve_lanes_refine(&ex, 3, |x: &Fl<QV, B>| ex - d.apply(x));
    let rhs = qu_field(&bq, &Fl::<QV, B>::zero());
    let resid = |scale: QD| {
        let (e, d, rhs) = (&e, &d, &rhs);
        move |x: &Fl<Qu, B>| {
            let ev = e.apply(&Fl::<QV, B>(x.0.map(|q| q.v)));
            let s = QV::splat(scale);
            *rhs - d.apply(x)
                - Fl(std::array::from_fn(|k| Qu {
                    v: QV::zero(),
                    d: ev.0[k] * s,
                }))
        }
    };
    let x3 = lu.solve_lanes_refine_t(&lu, &rhs, 3, resid(QD::ONE));
    let (ev, et) = (
        qd_err(&qu_val(&x3), &qvalues(&xv)),
        qd_err(&qu_tan(&x3), &qvalues(&t_ref)),
    );
    assert!(ev <= 1e-60 && et <= 1e-58, "value {ev:e}, tangent {et:e}");
    let x1 = lu.solve_lanes_refine_t(&lu, &rhs, 1, resid(QD::ONE));
    let e1 = qd_err(&qu_tan(&x1), &qvalues(&t_ref));
    assert!(e1 > 1e-40, "control: one step gives the tangent to {e1:e}");
    let xw = lu.solve_lanes_refine_t(&lu, &rhs, 3, resid(QD::from_f64(1.0).add_f(1e-30)));
    let ew = qd_err(&qu_tan(&xw), &qvalues(&t_ref));
    assert!(ew > 1e-58, "control: E scaled by 1 + 1e-30 passes ({ew:e})");
}

/// the order after which `TangentTail` ends a series whose tangent terms are 0.75 (r s)^k (s = 2) on two fields: the
/// first k with run orders in a row of (r s)^k <= tol, and at least kmin; None if never within 60 orders
fn tail_order<T: Lane>(mk: impl Fn(f64) -> Fl<T, B>, r: f64, tol: f64, run: usize, kmin: usize) -> Option<usize> {
    let mut tail = TangentTail::new(r, tol, run, kmin);
    (0..60).find(|&k| {
        let (f, g) = (mk(0.75 * 2f64.powi(k as i32)), mk(0.5 * 2f64.powi(k as i32)));
        tail.done(k, &[&f, &g])
    })
}

/// value 1, tangent t (every lane, the last f64 tangent where the lane has them)
fn with_tangent<T: Lane>(t: f64) -> Fl<T, B> {
    let mut f = Fl::<T, B>::from_vt(&Fl::splat(DD::ONE), &Fl::splat(DD::from_f64(t)));
    if T::TAN_F64 {
        for b in 0..B {
            f.0[b].set_tan64(T::NT - 1, multidouble::v8::V8::splat(t));
        }
    }
    f
}

/// Terms 0.75 2^-k (r = 1/4), tol 2^-10, run 3: the run starts at k = 10 and the series ends after k = 12, on the
/// double-double tangents (`Du`), the f64 tangents (`Fu`, `Dtn<2>`, `Fun<2>`) and the quad-double ones (`Qu`); kmin 20
/// moves the end to k = 20. Controls: a value lane type (`DV`) never ends; r = 1/2 (constant terms) never ends;
/// tol 0 never ends; a tangent seeded 2^20 larger on the last order breaks the run (the series ends 3 orders later).
#[test]
fn tangent_tail_ends_after_the_predicted_order() {
    let (r, tol) = (0.25, 2f64.powi(-10));
    assert_eq!(tail_order(with_tangent::<Du>, r, tol, 3, 2), Some(12));
    assert_eq!(tail_order(with_tangent::<Fu>, r, tol, 3, 2), Some(12));
    assert_eq!(tail_order(with_tangent::<Dtn<2>>, r, tol, 3, 2), Some(12));
    assert_eq!(tail_order(with_tangent::<Fun<2>>, r, tol, 3, 2), Some(12));
    assert_eq!(tail_order(with_tangent::<Qu>, r, tol, 3, 2), Some(12));
    assert_eq!(tail_order(with_tangent::<Du>, r, tol, 3, 20), Some(20));
    assert_eq!(tail_order(with_tangent::<Du>, r, tol, 1, 2), Some(10));
    assert_eq!(tail_order(with_tangent::<DV>, r, tol, 3, 2), None, "control: values");
    assert_eq!(
        tail_order(with_tangent::<Du>, 0.5, tol, 3, 2),
        None,
        "control: constant terms"
    );
    assert_eq!(tail_order(with_tangent::<Du>, r, 0.0, 3, 2), None, "control: tol 0");
    let mut tail = TangentTail::new(r, tol, 3, 2);
    let mut end = None;
    for k in 0..60 {
        let f = with_tangent::<Du>(0.75 * 2f64.powi(k as i32) * if k == 12 { 2f64.powi(20) } else { 1.0 });
        if tail.done(k, &[&f]) {
            end = Some(k);
            break;
        }
    }
    assert_eq!(
        end,
        Some(15),
        "control: a large term at k = 12 does not restart the run"
    );
}

// ---------------------------------------------------------------------------------------------------------------
// the order-by-order jet solves
// ---------------------------------------------------------------------------------------------------------------

/// A jet system (omega(h) D + diag kappa(h)) x(h) = rhs(h), the form of a regular-horizon series' systems: D a
/// constant matrix with entries below 0.05, omega(h) a series constant over the values with omega_0 in [1, 2), and
/// kappa(h) a field with kappa_0 in [2, 3), so that A_0 = omega_0 D + diag kappa_0 is diagonally dominant. The
/// coefficients are quad-doubles with all four parts in use, or double-doubles. `lu` is the LU of A_0 rounded to
/// double-double.
struct JetSystem {
    d: DMat<Native, B>,
    /// D by rows
    dd: Vec<DD>,
    omega: Vec<QD>,
    kappa: Vec<Vec<QD>>,
    lu: DDLu,
}

/// c + w u (u uniform in [-1, 1)) with all four parts in use, or rounded to double-double
fn jet_coef(r: &mut Rng, c: f64, w: f64, qd: bool) -> QD {
    let h = c + w * r.sym();
    let a = h.abs();
    let q = QD::from_dd(DD::new(h, r.sym() * a * 2f64.powi(-53)))
        .add_dd(DD::new(r.sym() * a * 2f64.powi(-106), r.sym() * a * 2f64.powi(-159)));
    if qd {
        q
    } else {
        QD::from_dd(q.to_dd())
    }
}

/// n coefficients of a field, each value in [-1, 1)
fn jet_coefs(r: &mut Rng, n: usize, qd: bool) -> Vec<Vec<QD>> {
    (0..n)
        .map(|_| (0..M).map(|_| jet_coef(r, 0.0, 1.0, qd)).collect())
        .collect()
}

fn jet_system(r: &mut Rng, n: usize, qd: bool) -> JetSystem {
    let dd: Vec<DD> = (0..M * M)
        .map(|_| DD::new(0.05 * r.sym(), 0.05 * r.sym() * 1e-17))
        .collect();
    let (h, l): (Vec<f64>, Vec<f64>) = dd.iter().map(|d| (d.hi, d.lo)).unzip();
    let omega: Vec<QD> = (0..n)
        .map(|k| {
            if k == 0 {
                jet_coef(r, 1.5, 0.5, qd)
            } else {
                jet_coef(r, 0.0, 1.0, qd)
            }
        })
        .collect();
    let kappa: Vec<Vec<QD>> = (0..n)
        .map(|k| {
            (0..M)
                .map(|_| {
                    if k == 0 {
                        jet_coef(r, 2.5, 0.5, qd)
                    } else {
                        jet_coef(r, 0.0, 1.0, qd)
                    }
                })
                .collect()
        })
        .collect();
    let ld = ld_for(M);
    let (mut ch, mut cl) = (vec![0.0; ld * M], vec![0.0; ld * M]);
    for i in 0..M {
        for j in 0..M {
            let mut a = omega[0].to_dd() * dd[i * M + j];
            if i == j {
                a = a + kappa[0][i].to_dd();
            }
            ch[j * ld + i] = a.hi;
            cl[j * ld + i] = a.lo;
        }
    }
    let lu = DDLu::factor(M, ld, ch, cl, false).unwrap();
    JetSystem {
        d: DMat::from_rows(&h, &l),
        dd,
        omega,
        kappa,
        lu,
    }
}

/// a jet constant over the values (`Lane::set_coef_qd`, so rounded to double-double on `Jet`)
fn jet_scalar<T: Lane<S = Native>>(c: &[QD]) -> T {
    let mut t = T::zero();
    for (k, &q) in c.iter().enumerate() {
        t.set_coef_qd(k, QV::splat(q));
    }
    t
}

/// a jet field with the values c[k] as its coefficient k
fn jet_field<T: Lane<S = Native>>(c: &[Vec<QD>]) -> Fl<T, B> {
    let mut f = Fl::<T, B>::zero();
    for (k, v) in c.iter().enumerate() {
        f.set_coef_qd(k, &qfield(v));
    }
    f
}

/// the bits of every coefficient in quad-double (exact for `Jet`), order by order
fn jet_bits<T: Lane<S = Native>>(f: &Fl<T, B>) -> Vec<u64> {
    (0..T::NJ)
        .flat_map(|k| qvalues(&f.coef_qd(k)).into_iter().flat_map(|q| q.x.map(f64::to_bits)))
        .collect()
}

fn dv_bits(f: &Fl<DV, B>) -> Vec<u64> {
    (0..M)
        .flat_map(|i| [f.get(i).hi.to_bits(), f.get(i).lo.to_bits()])
        .collect()
}

/// FNV-1a
fn fold_bits(h: u64, w: &[u64]) -> u64 {
    w.iter().fold(h, |h, &w| (h ^ w).wrapping_mul(0x0000_0100_0000_01b3))
}

/// psolve's double-double jet branch (choptuik-large-d, hp/src/model.rs, multidouble 2.0.0), the code
/// `solve_lanes_jet` is to replace bit for bit
fn psolve_jet<T: Lane<S = Native>>(
    d: &DMat<Native, B>,
    lu: &DDLu,
    kappa: &Fl<T, B>,
    rhs: &Fl<T, B>,
    omega: T,
) -> Fl<T, B> {
    let mut f = Fl::<T, B>::zero();
    let mut fs: Vec<Fl<DV, B>> = Vec::with_capacity(T::NJ);
    for k in 0..T::NJ {
        let mut r = rhs.coef(k);
        for j in 1..=k {
            let dfv = d.apply::<DV>(&fs[k - j]);
            let od = Fl::<DV, B>([omega.coef(j); B]);
            r = (r - dfv * od) - kappa.coef(j) * fs[k - j];
        }
        let fk = lu.solve_lanes(&r);
        f.set_coef(k, &fk);
        fs.push(fk);
    }
    f
}

/// psolve's quad-double jet branch, likewise for `solve_lanes_jet_refine`
fn psolve_qjet<T: Lane<S = Native>>(
    d: &DMat<Native, B>,
    lu: &DDLu,
    kappa: &Fl<T, B>,
    rhs: &Fl<T, B>,
    omega: T,
    steps: usize,
) -> Fl<T, B> {
    let om: Vec<QV> = (0..T::NJ).map(|j| omega.coef_qd(j)).collect();
    let k0 = kappa.coef_qd(0);
    let mut f = Fl::<T, B>::zero();
    let mut fs: Vec<Fl<QV, B>> = Vec::with_capacity(T::NJ);
    for k in 0..T::NJ {
        let mut r = rhs.coef_qd(k);
        for j in 1..=k {
            r = (r - d.apply(&fs[k - j]).mul_t(om[j])) - kappa.coef_qd(j) * fs[k - j];
        }
        let fk = lu.solve_lanes_refine(&r, steps, |x| r - (d.apply(x).mul_t(om[0]) + k0 * *x));
        f.set_coef_qd(k, &fk);
        fs.push(fk);
    }
    f
}

/// `solve_lanes_jet` with psolve's coupling
fn solve_dd<T: Lane<S = Native>>(s: &JetSystem, kappa: &Fl<T, B>, rhs: &Fl<T, B>, omega: T) -> Fl<T, B> {
    s.lu.solve_lanes_jet(rhs, |j, x, r| {
        (r - s.d.apply::<DV>(x) * Fl([omega.coef(j); B])) - kappa.coef(j) * *x
    })
}

/// `solve_lanes_jet_refine` with psolve's coupling and residual
fn solve_qd<T: Lane<S = Native>>(s: &JetSystem, kappa: &Fl<T, B>, rhs: &Fl<T, B>, omega: T, steps: usize) -> Fl<T, B> {
    let om: Vec<QV> = (0..T::NJ).map(|j| omega.coef_qd(j)).collect();
    let k0 = kappa.coef_qd(0);
    s.lu.solve_lanes_jet_refine(
        rhs,
        steps,
        |j, x, r| (r - s.d.apply(x).mul_t(om[j])) - kappa.coef_qd(j) * *x,
        |r, x| *r - (s.d.apply(x).mul_t(om[0]) + k0 * *x),
    )
}

/// One random system of N orders through both solves and both of psolve's loops, bit for bit, the outputs folded
/// into h. Returns whether the coupling reassociated, r - (omega_j D x + kappa_j x), changed a bit (double-double,
/// quad-double).
fn psolve_case<const N: usize>(r: &mut Rng, h: &mut [u64; 2]) -> [bool; 2] {
    let s = jet_system(r, N, false);
    let (kappa, omega) = (jet_field::<Jet<N>>(&s.kappa), jet_scalar::<Jet<N>>(&s.omega));
    let c = jet_coefs(r, N, false);
    let rhs = jet_field::<Jet<N>>(&c);
    let x = solve_dd(&s, &kappa, &rhs, omega);
    assert_eq!(
        jet_bits(&x),
        jet_bits(&psolve_jet(&s.d, &s.lu, &kappa, &rhs, omega)),
        "Jet<{N}>"
    );
    h[0] = fold_bits(h[0], &jet_bits(&x));
    // a quad-double jet of the same double-double data is solved in double-double alike
    let (kq, oq) = (jet_field::<QJet<N>>(&s.kappa), jet_scalar::<QJet<N>>(&s.omega));
    assert_eq!(
        jet_bits(&solve_dd(&s, &kq, &jet_field::<QJet<N>>(&c), oq)),
        jet_bits(&x),
        "QJet<{N}> in dd"
    );
    let re = s.lu.solve_lanes_jet(&rhs, |j, x, r| {
        r - (s.d.apply::<DV>(x) * Fl([omega.coef(j); B]) + kappa.coef(j) * *x)
    });
    let moved_dd = jet_bits(&re) != jet_bits(&x);
    let s = jet_system(r, N, true);
    let (kappa, omega) = (jet_field::<QJet<N>>(&s.kappa), jet_scalar::<QJet<N>>(&s.omega));
    let rhs = jet_field::<QJet<N>>(&jet_coefs(r, N, true));
    for steps in [0, 2, 3] {
        let x = solve_qd(&s, &kappa, &rhs, omega, steps);
        let want = psolve_qjet(&s.d, &s.lu, &kappa, &rhs, omega, steps);
        assert_eq!(jet_bits(&x), jet_bits(&want), "QJet<{N}>, {steps} steps");
        h[1] = fold_bits(h[1], &jet_bits(&x));
    }
    let om: Vec<QV> = (0..N).map(|j| omega.coef_qd(j)).collect();
    let k0 = kappa.coef_qd(0);
    let re = s.lu.solve_lanes_jet_refine(
        &rhs,
        2,
        |j, x, r| r - (s.d.apply(x).mul_t(om[j]) + kappa.coef_qd(j) * *x),
        |r, x| *r - (s.d.apply(x).mul_t(om[0]) + k0 * *x),
    );
    [
        moved_dd,
        jet_bits(&re) != jet_bits(&solve_qd(&s, &kappa, &rhs, omega, 2)),
    ]
}

/// The jet solves with psolve's coupling are psolve's own jet loops (choptuik-large-d's hp/src/model.rs, written out
/// above) bit for bit, on `Jet<N>` and on `QJet<N>` with 0, 2 and 3 refinement steps, N = 2, 5 and 8, so that hp can
/// call them without a change in its results; a quad-double jet in `solve_lanes_jet` is the double-double one on its
/// rounded coefficients. The outputs' bits are pinned (hashes recorded with the first implementation; every backend
/// gives these bits). Control: the coupling reassociated, r - (omega_j D x + kappa_j x), changes the bits.
#[test]
fn jet_solves_are_the_psolve_loops_bitwise() {
    let mut r = Rng(81);
    let mut h = [0xcbf2_9ce4_8422_2325u64; 2];
    let mut moved = [false; 2];
    for _ in 0..3 {
        for (i, m) in [
            psolve_case::<2>(&mut r, &mut h),
            psolve_case::<5>(&mut r, &mut h),
            psolve_case::<8>(&mut r, &mut h),
        ]
        .concat()
        .into_iter()
        .enumerate()
        {
            moved[i % 2] |= m;
        }
    }
    let want: [u64; 2] = [0x2f44_6d53_8800_2abe, 0x2e99_109e_8c6e_ed24];
    assert!(
        h == want,
        "pinned: {:#018x} {:#018x}, want {:#018x} {:#018x}",
        h[0],
        h[1],
        want[0],
        want[1]
    );
    assert!(
        moved == [true, true],
        "control: a reassociated coupling gives the same bits ({moved:?})"
    );
}

/// a planted solution x(h) and rhs(h) = A(h) x(h) formed exactly (dyadic arithmetic) and rounded to quad-double, or
/// to double-double, sum_{j=0..k} (omega_j D + diag kappa_j) x_{k-j} at order k
fn jet_planted(r: &mut Rng, s: &JetSystem, n: usize, qd: bool) -> (Vec<Vec<QD>>, Vec<Vec<QD>>) {
    let x = jet_coefs(r, n, qd);
    let rhs = (0..n)
        .map(|k| {
            (0..M)
                .map(|i| {
                    let mut acc = Dy::zero();
                    for j in 0..=k {
                        let xs = &x[k - j];
                        let mut dx = Dy::zero();
                        for l in 0..M {
                            let dil = Dy::from_f64(s.dd[i * M + l].hi).add(&Dy::from_f64(s.dd[i * M + l].lo));
                            dx = dx.add(&dil.mul(&xs[l].exact()));
                        }
                        acc = acc.add(&s.omega[j].exact().mul(&dx));
                        acc = acc.add(&s.kappa[j][i].exact().mul(&xs[i].exact()));
                    }
                    let g = acc.greedy4();
                    if qd {
                        QD { x: g }
                    } else {
                        QD::from_dd(DD::new(g[0], g[1]))
                    }
                })
                .collect()
        })
        .collect();
    (x, rhs)
}

/// worst over the orders of max_i |x_k,i - planted_k,i| / max_i |planted_k,i|
fn jet_err<T: Lane<S = Native>>(f: &Fl<T, B>, planted: &[Vec<QD>]) -> f64 {
    (0..T::NJ)
        .map(|k| qd_err(&qvalues(&f.coef_qd(k)), &planted[k]))
        .fold(0.0, f64::max)
}

/// Planted solutions of jet systems against the solves with psolve's coupling: `Jet<8>` within 1e-29 of the
/// planted coefficients (double-double rounding of the right-hand side and the solves, carried through the orders;
/// 5e-32 seen), also through the refined solve, and `QJet<8>` with two refinement steps within 1e-60 (1e-64 seen).
/// Controls: the coupling with its sign flipped, r + A_j x, and the coupling without its last term, j = k skipped,
/// miss the double-double bound; no refinement (steps 0) misses the quad-double bound, and so does the quad-double
/// jet through the double-double solve.
#[test]
fn jet_solves_reach_the_planted_solution() {
    let mut r = Rng(82);
    let (mut worst_dd, mut worst_qd) = (0.0f64, 0.0f64);
    for _ in 0..3 {
        let s = jet_system(&mut r, 8, false);
        let (x0, c) = jet_planted(&mut r, &s, 8, false);
        let (kappa, omega, rhs) = (
            jet_field::<Jet<8>>(&s.kappa),
            jet_scalar::<Jet<8>>(&s.omega),
            jet_field::<Jet<8>>(&c),
        );
        let e = jet_err(&solve_dd(&s, &kappa, &rhs, omega), &x0);
        assert!(e <= 1e-29, "Jet<8>: {e:e}");
        worst_dd = worst_dd.max(e);
        let flip = s.lu.solve_lanes_jet(&rhs, |j, x, r| {
            (r + s.d.apply::<DV>(x) * Fl([omega.coef(j); B])) + kappa.coef(j) * *x
        });
        let e = jet_err(&flip, &x0);
        assert!(e > 1e-29, "control: the coupling's sign flipped passes ({e:e})");
        // per order k the calls are j = 1, ..., k, so the last skips A_k x_0
        let mut calls = 0;
        let skip = s.lu.solve_lanes_jet(&rhs, |j, x, r| {
            calls += 1;
            if j > 1 && calls % j == 0 {
                return r;
            }
            (r - s.d.apply::<DV>(x) * Fl([omega.coef(j); B])) - kappa.coef(j) * *x
        });
        let e = jet_err(&skip, &x0);
        assert!(e > 1e-29, "control: the coupling without its last term passes ({e:e})");
        // a double-double jet through the refined solve, its coefficients rounded back to double-double
        let e = jet_err(&solve_qd(&s, &kappa, &rhs, omega, 2), &x0);
        assert!(e <= 1e-29, "Jet<8> refined: {e:e}");

        let s = jet_system(&mut r, 8, true);
        let (x0, c) = jet_planted(&mut r, &s, 8, true);
        let (kappa, omega, rhs) = (
            jet_field::<QJet<8>>(&s.kappa),
            jet_scalar::<QJet<8>>(&s.omega),
            jet_field::<QJet<8>>(&c),
        );
        let e = jet_err(&solve_qd(&s, &kappa, &rhs, omega, 2), &x0);
        assert!(e <= 1e-60, "QJet<8>: {e:e}");
        worst_qd = worst_qd.max(e);
        let e0 = jet_err(&solve_qd(&s, &kappa, &rhs, omega, 0), &x0);
        assert!(e0 > 1e-60, "control: the unrefined solve passes ({e0:e})");
        let ed = jet_err(&solve_dd(&s, &kappa, &rhs, omega), &x0);
        assert!(ed > 1e-60, "control: the double-double solve passes ({ed:e})");
    }
    println!("worst error: Jet<8> {worst_dd:.1e}, QJet<8> {worst_qd:.1e}");
}

/// `solve_lanes_jet` on `Jet<2>` is the solve of a system with a tangent on `Du` lanes, value x_v = A_0^-1 r_v and
/// tangent x_d = A_0^-1 (r_d - omega_d D x_v - kappa_d x_v), bit for bit. Control: the jet's coupling with the sign
/// of omega_1 D x flipped differs.
#[test]
fn jet_solve_of_order_one_is_the_dual_solve() {
    let mut r = Rng(83);
    for _ in 0..3 {
        let s = jet_system(&mut r, 2, false);
        let c = jet_coefs(&mut r, 2, false);
        let (kappa, omega, rhs) = (
            jet_field::<Jet<2>>(&s.kappa),
            jet_scalar::<Jet<2>>(&s.omega),
            jet_field::<Jet<2>>(&c),
        );
        let x = solve_dd(&s, &kappa, &rhs, omega);
        let kd = Fl::<Du, B>::from_vt(&kappa.coef(0), &kappa.coef(1));
        let od = Du::from_vt(omega.coef(0), omega.coef(1));
        let rd = Fl::<Du, B>::from_vt(&rhs.coef(0), &rhs.coef(1));
        let fv = s.lu.solve_lanes(&rd.val());
        let rt = (rd.tan() - s.d.apply::<DV>(&fv) * Fl([od.tan(); B])) - kd.tan() * fv;
        let fd = s.lu.solve_lanes(&rt);
        assert_eq!(dv_bits(&x.coef(0)), dv_bits(&fv));
        assert_eq!(dv_bits(&x.coef(1)), dv_bits(&fd));
        let flip = s.lu.solve_lanes_jet(&rhs, |j, x, r| {
            (r + s.d.apply::<DV>(x) * Fl([omega.coef(j); B])) - kappa.coef(j) * *x
        });
        assert_ne!(dv_bits(&flip.coef(1)), dv_bits(&fd), "control");
    }
}

/// The jet solves refuse a lane type without jet coefficients
#[test]
#[should_panic(expected = "FieldSolve::solve_lanes_jet solves jets")]
fn jet_solve_refuses_a_dual_lane() {
    let lu = lu_of_order(M);
    let rhs = Fl::<Du, B>([Du::splat(DD::ONE); B]);
    lu.solve_lanes_jet(&rhs, |_, _, r| r);
}

#[test]
#[should_panic(expected = "FieldSolve::solve_lanes_jet_refine solves jets")]
fn jet_solve_refine_refuses_a_value_lane() {
    let lu = lu_of_order(M);
    let rhs = qfield(&[QD::ONE; M]);
    lu.solve_lanes_jet_refine(&rhs, 2, |_, _, r| r, |r, _| *r);
}