multidouble 3.2.0

//! Fields: N = 8 B values of a lane type as B lane vectors (`Fl<T, B>`; a function on a grid of 8 B points, or 8 B
//! independent problems), with the power-series arithmetic that series methods for ODEs need on them: Cauchy products
//! of coefficient sequences (`cs`, `cs_range`, `cs_w`), Horner evaluation (`horner`), and a constant double-double
//! matrix applied to a field (`DMat`, for instance a spectral derivative). Every operation is the lane type's, so code
//! written over `Fl<T, B>` runs with values, dual numbers or quad-double lanes alike (`dual::Lane`).
//!
//! The products accumulate with the lane type's dot accumulator (`Lane::acc`), in the order of the index; for the f64
//! lanes (`VAL_F64` or `TAN_F64`) the B lane vectors' accumulations are interleaved, since a single f64 chain is
//! latency-bound, and each still sums its terms in index order. `cs_into`, `cs_range_into` and `cs_w_into` write the
//! same sums, bitwise, into a caller's field, with the accumulator and the terms borrowed (`Lane::acc_mut`), which
//! for jet lanes saves copying kilobytes per term. `cs_lane_into` and `cs_w_lane_into` are one lane vector of them,
//! for a caller that sums the lane vectors on separate threads. A matrix product is `MatLane::dmat_apply`: the
//! double-double columns through the accumulator by default, the columns' f64 copy for the f64 lanes.
use crate::dd::DD;
use crate::dual::{Counted, Dc, Dtn, Du, Fu, Fun, Lane, Qu, F8};
use crate::dv::{DotAcc, DV};
use crate::jet::{Jet, QJet};
use crate::linalg::DDLu;
use crate::qd::QD;
use crate::qv::QV;
use crate::simd::Simd;
use crate::v8::V8;
use std::ops::{Add, Div, Mul, Neg, Sub};

/// 8 B values of the lane type T as B lane vectors: value j is lane j % 8 of vector j / 8.
#[derive(Clone, Copy, Debug)]
pub struct Fl<T: Lane, const B: usize>(pub [T; B]);

/// The lane type's double-double vector (constants and values of `Fl<T, B>` come as these)
type Dv<T> = DV<<T as Lane>::S>;

macro_rules! fl_binop {
    ($tr:ident, $f:ident, $op:tt) => {
        impl<T: Lane, const B: usize> $tr for Fl<T, B> {
            type Output = Fl<T, B>;
            #[inline(always)]
            fn $f(self, o: Fl<T, B>) -> Fl<T, B> {
                let mut r = self.0;
                for b in 0..B {
                    r[b] = self.0[b] $op o.0[b];
                }
                Fl(r)
            }
        }
    };
}
fl_binop!(Add, add, +);
fl_binop!(Sub, sub, -);
fl_binop!(Mul, mul, *);
fl_binop!(Div, div, /);

impl<T: Lane, const B: usize> Neg for Fl<T, B> {
    type Output = Fl<T, B>;
    #[inline(always)]
    fn neg(self) -> Fl<T, B> {
        let mut r = self.0;
        for b in 0..B {
            r[b] = -self.0[b];
        }
        Fl(r)
    }
}

impl<T: Lane, const B: usize> Fl<T, B> {
    #[inline(always)]
    pub fn zero() -> Self {
        Fl([T::zero(); B])
    }
    #[inline(always)]
    pub fn splat(x: DD) -> Self {
        Fl([T::splat(x); B])
    }
    /// a constant field (zero tangents)
    #[inline(always)]
    pub fn cst(c: &Fl<Dv<T>, B>) -> Self {
        let mut r = [T::zero(); B];
        for b in 0..B {
            r[b] = T::cst(c.0[b]);
        }
        Fl(r)
    }
    #[inline(always)]
    pub fn muls(self, s: DD) -> Self {
        let sv = DV::splat(s);
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b].mulc(sv);
        }
        Fl(r)
    }
    #[inline(always)]
    pub fn mulf(self, s: f64) -> Self {
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b].mulf(s);
        }
        Fl(r)
    }
    /// times a power of two, every part exact (`Lane::mul_pwr2`): bitwise `mulf(s)` for such s
    #[inline(always)]
    pub fn mul_pwr2(self, s: f64) -> Self {
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b].mul_pwr2(s);
        }
        Fl(r)
    }
    /// times a lane-constant scalar of the lane type (e.g. a frequency that carries a tangent)
    #[inline(always)]
    pub fn mul_t(self, s: T) -> Self {
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b] * s;
        }
        Fl(r)
    }
    #[inline(always)]
    pub fn adds(self, s: DD) -> Self {
        let sv = DV::splat(s);
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b].addc(sv);
        }
        Fl(r)
    }
    pub fn val(&self) -> Fl<Dv<T>, B> {
        let mut r = [DV::zero(); B];
        for b in 0..B {
            r[b] = self.0[b].val();
        }
        Fl(r)
    }
    /// times a parameter constant (`Lane::mulp`)
    pub fn mulp(self, p: T::P) -> Self {
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b].mulp(p);
        }
        Fl(r)
    }
    /// plus a parameter constant (`Lane::addp`)
    pub fn addp(self, p: T::P) -> Self {
        let mut r = self.0;
        for b in 0..B {
            r[b] = self.0[b].addp(p);
        }
        Fl(r)
    }
    /// coefficient k of every lane (`Lane::coef`)
    pub fn coef(&self, k: usize) -> Fl<Dv<T>, B> {
        let mut r = [DV::zero(); B];
        for b in 0..B {
            r[b] = self.0[b].coef(k);
        }
        Fl(r)
    }
    /// set coefficient k of every lane (`Lane::set_coef`, jets)
    pub fn set_coef(&mut self, k: usize, c: &Fl<Dv<T>, B>) {
        for b in 0..B {
            self.0[b].set_coef(k, c.0[b]);
        }
    }
    /// coefficient k of every lane in quad-double (`Lane::coef_qd`)
    pub fn coef_qd(&self, k: usize) -> Fl<QV<T::S>, B> {
        let mut r = [QV::zero(); B];
        for b in 0..B {
            r[b] = self.0[b].coef_qd(k);
        }
        Fl(r)
    }
    /// set coefficient k of every lane from quad-double (`Lane::set_coef_qd`, jets)
    pub fn set_coef_qd(&mut self, k: usize, c: &Fl<QV<T::S>, B>) {
        for b in 0..B {
            self.0[b].set_coef_qd(k, c.0[b]);
        }
    }
    pub fn tan(&self) -> Fl<Dv<T>, B> {
        let mut r = [DV::zero(); B];
        for b in 0..B {
            r[b] = self.0[b].tan();
        }
        Fl(r)
    }
    pub fn from_vt(v: &Fl<Dv<T>, B>, t: &Fl<Dv<T>, B>) -> Self {
        let mut r = [T::zero(); B];
        for b in 0..B {
            r[b] = T::from_vt(v.0[b], t.0[b]);
        }
        Fl(r)
    }
    /// f64 tangent k of every value, as an `F8` field (the `TAN_F64` lane types)
    pub fn tan64(&self, k: usize) -> Fl<F8<T::S>, B> {
        let mut r = [F8(V8::zero()); B];
        for b in 0..B {
            r[b] = F8(self.0[b].tan64(k));
        }
        Fl(r)
    }
    /// set f64 tangent k of every value from an `F8` field, the inverse of `tan64`, on the `TAN_F64` lane types
    /// (`Fu`, `Dtn`, `Fun`) with k < `Lane::NT`. Any other lane type, or a k past the last tangent, panics.
    pub fn set_tan64(&mut self, k: usize, t: &Fl<F8<T::S>, B>) {
        assert!(
            T::TAN_F64 && k < T::NT,
            "Fl::set_tan64: f64 tangent {k} of a lane type with {} f64 tangents",
            T::NT
        );
        for b in 0..B {
            self.0[b].set_tan64(k, t.0[b].0);
        }
    }
    /// the unit field e_j, one at value j and zero at the others, as constants of the lane type (zero tangents)
    pub fn unit(j: usize) -> Self {
        assert!(j < 8 * B, "Fl::unit: value {j} of a field of {} values", 8 * B);
        let mut e = [0.0; 8];
        e[j % 8] = 1.0;
        let mut r = Self::zero();
        r.0[j / 8] = T::cst(DV {
            hi: V8::from_array(e),
            lo: V8::zero(),
        });
        r
    }
    /// values v with the tangent e_j (`unit`), the seed of a Jacobian column in forward mode. It is `from_vt(v, e_j)`,
    /// so it sets the tangent of `Du` and `Qu`, the real part of `Dc`'s, order 1 of a jet, and the first f64 tangent
    /// of `Fu`, `Dtn` and `Fun`, whose other tangents `set_tan64` with `unit` seeds. A lane type without a tangent
    /// panics.
    pub fn unit_tangent(v: &Fl<Dv<T>, B>, j: usize) -> Self {
        assert!(T::TANGENT, "Fl::unit_tangent: a lane type without a tangent");
        Self::from_vt(v, &Fl::<Dv<T>, B>::unit(j))
    }
}

impl<S: Simd, const B: usize> Fl<Dc<S>, B> {
    /// the real part of the tangent (`Lane::tan`)
    pub fn tan_re(&self) -> Fl<DV<S>, B> {
        let mut r = [DV::zero(); B];
        for b in 0..B {
            r[b] = self.0[b].r;
        }
        Fl(r)
    }
    /// the imaginary part of the tangent
    pub fn tan_im(&self) -> Fl<DV<S>, B> {
        let mut r = [DV::zero(); B];
        for b in 0..B {
            r[b] = self.0[b].i;
        }
        Fl(r)
    }
    /// a field from its values and the real and imaginary parts of its tangent
    pub fn from_vri(v: &Fl<DV<S>, B>, re: &Fl<DV<S>, B>, im: &Fl<DV<S>, B>) -> Self {
        let mut r = [Dc::zero(); B];
        for b in 0..B {
            r[b] = Dc {
                v: v.0[b],
                r: re.0[b],
                i: im.0[b],
            };
        }
        Fl(r)
    }
}

impl<S: Simd, const B: usize> Fl<Qu<S>, B> {
    /// a field from quad-double values and tangents, exactly (`from_vt` takes them in double-double), for instance the
    /// values with the tangent `Fl::<QV<S>, B>::unit(j)`, a quad-double Jacobian column's seed. `coef_qd(0)` and
    /// `coef_qd(1)` read them back.
    pub fn from_vt_qd(v: &Fl<QV<S>, B>, t: &Fl<QV<S>, B>) -> Self {
        let mut r = [Qu::zero(); B];
        for b in 0..B {
            r[b] = Qu { v: v.0[b], d: t.0[b] };
        }
        Fl(r)
    }
}

impl<S: Simd, const B: usize> Fl<DV<S>, B> {
    /// from 8 B values (hi, lo)
    pub fn load(hi: &[f64], lo: &[f64]) -> Self {
        let mut r = [DV::zero(); B];
        for b in 0..B {
            r[b] = DV::load(&hi[8 * b..8 * b + 8], &lo[8 * b..8 * b + 8]);
        }
        Fl(r)
    }
    pub fn store(&self, hi: &mut [f64], lo: &mut [f64]) {
        for b in 0..B {
            self.0[b].store(&mut hi[8 * b..8 * b + 8], &mut lo[8 * b..8 * b + 8]);
        }
    }
    pub fn get(&self, j: usize) -> DD {
        self.0[j / 8].lane(j % 8)
    }
    /// max |hi| over the values
    pub fn amax(&self) -> f64 {
        let mut m = 0.0f64;
        for b in 0..B {
            m = m.max(self.0[b].hi.abs().hmax());
        }
        m
    }
}

impl<S: Simd, const B: usize> Fl<QV<S>, B> {
    /// from 8 B quad-double values given as their four parts' arrays (c0 + c1 + c2 + c3 each, `QV::load` per lane
    /// vector: the parts as they are, not renormalized)
    pub fn load_parts(c0: &[f64], c1: &[f64], c2: &[f64], c3: &[f64]) -> Self {
        let mut r = [QV::zero(); B];
        for b in 0..B {
            let o = 8 * b;
            r[b] = QV::load(&c0[o..o + 8], &c1[o..o + 8], &c2[o..o + 8], &c3[o..o + 8]);
        }
        Fl(r)
    }
    /// the 8 B values' four parts into four arrays
    pub fn store_parts(&self, c0: &mut [f64], c1: &mut [f64], c2: &mut [f64], c3: &mut [f64]) {
        for b in 0..B {
            let o = 8 * b;
            self.0[b].store(
                &mut c0[o..o + 8],
                &mut c1[o..o + 8],
                &mut c2[o..o + 8],
                &mut c3[o..o + 8],
            );
        }
    }
    pub fn get(&self, j: usize) -> QD {
        self.0[j / 8].lane(j % 8)
    }
}

/// Cauchy product coefficient `sum_{j=0..k} x[j] y[k-j]`
#[inline(always)]
pub fn cs<T: Lane, const B: usize>(x: &[Fl<T, B>], y: &[Fl<T, B>], k: usize) -> Fl<T, B> {
    if T::VAL_F64 || T::TAN_F64 {
        return cs_il(x, y, k, 0, k, |_| 1.0, false);
    }
    let mut out = [T::zero(); B];
    for b in 0..B {
        let mut a = T::acc0();
        for j in 0..=k {
            a = T::acc(a, x[j].0[b], y[k - j].0[b]);
        }
        out[b] = T::fin(a);
    }
    Fl(out)
}

/// `sum_{j=lo..=hi} x[j] y[k-j]`
#[inline(always)]
pub fn cs_range<T: Lane, const B: usize>(x: &[Fl<T, B>], y: &[Fl<T, B>], k: usize, lo: usize, hi: usize) -> Fl<T, B> {
    if T::VAL_F64 || T::TAN_F64 {
        return cs_il(x, y, k, lo, hi, |_| 1.0, false);
    }
    let mut out = [T::zero(); B];
    for b in 0..B {
        let mut a = T::acc0();
        for j in lo..=hi {
            a = T::acc(a, x[j].0[b], y[k - j].0[b]);
        }
        out[b] = T::fin(a);
    }
    Fl(out)
}

/// `sum_{j=lo..=hi} (w(j) x[j]) y[k-j]` with exact f64 weights w(j) (small integers, as in the coefficients of a
/// derivative)
#[inline(always)]
pub fn cs_w<T: Lane, const B: usize>(
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
    w: impl Fn(usize) -> f64,
) -> Fl<T, B> {
    if T::VAL_F64 || T::TAN_F64 {
        return cs_il(x, y, k, lo, hi, w, true);
    }
    let mut out = [T::zero(); B];
    for b in 0..B {
        let mut a = T::acc0();
        for j in lo..=hi {
            a = T::acc(a, x[j].0[b].mulf(w(j)), y[k - j].0[b]);
        }
        out[b] = T::fin(a);
    }
    Fl(out)
}

/// sum_{j=lo..=hi} (w(j) x[j]) y[k-j] (weights only if weighted) with the B lane vectors' accumulations interleaved,
/// for the f64 lanes, whose single accumulation chains are latency-bound; each lane vector still sums its terms in the
/// order j = lo..=hi.
#[inline(always)]
fn cs_il<T: Lane, const B: usize>(
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
    w: impl Fn(usize) -> f64,
    weighted: bool,
) -> Fl<T, B> {
    let mut a = [T::acc0(); B];
    for j in lo..=hi {
        let (xj, yj) = (&x[j].0, &y[k - j].0);
        if weighted {
            let wj = w(j);
            for b in 0..B {
                a[b] = T::acc(a[b], xj[b].mulf(wj), yj[b]);
            }
        } else {
            for b in 0..B {
                a[b] = T::acc(a[b], xj[b], yj[b]);
            }
        }
    }
    Fl(a.map(T::fin))
}

/// `cs` into out: bitwise `cs(x, y, k)`, with the accumulator and the terms borrowed (`Lane::acc_mut`) and the sum
/// written in place, so a jet lane's sum makes no copy of its kilobytes per term. out is borrowed apart from x and y;
/// to write a series' own order k, split the vector (`split_at_mut`).
///
/// ```
/// use multidouble::dd::DD;
/// use multidouble::field::{cs, cs_into, Fl};
/// use multidouble::native::QJet;
/// let x: Vec<Fl<QJet<3>, 1>> = (1..4).map(|j| Fl::splat(DD::from_f64(j as f64))).collect();
/// let mut out = Fl::zero();
/// cs_into(&mut out, &x, &x, 2);
/// assert_eq!(out.coef_qd(0).get(0), cs(&x, &x, 2).coef_qd(0).get(0));
/// ```
#[inline(always)]
pub fn cs_into<T: Lane, const B: usize>(out: &mut Fl<T, B>, x: &[Fl<T, B>], y: &[Fl<T, B>], k: usize) {
    cs_range_into(out, x, y, k, 0, k);
}

/// `cs_range` into out: bitwise `cs_range(x, y, k, lo, hi)`, in place as `cs_into`
#[inline(always)]
pub fn cs_range_into<T: Lane, const B: usize>(
    out: &mut Fl<T, B>,
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
) {
    if T::VAL_F64 || T::TAN_F64 {
        cs_il_into(out, x, y, k, lo, hi, |_| 1.0, false);
    } else {
        for b in 0..B {
            cs_lane_into(&mut out.0[b], b, x, y, k, lo, hi);
        }
    }
}

/// `cs_w` into out: bitwise `cs_w(x, y, k, lo, hi, w)`, in place as `cs_into`
#[inline(always)]
pub fn cs_w_into<T: Lane, const B: usize>(
    out: &mut Fl<T, B>,
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
    w: impl Fn(usize) -> f64,
) {
    if T::VAL_F64 || T::TAN_F64 {
        cs_il_into(out, x, y, k, lo, hi, w, true);
    } else {
        for b in 0..B {
            cs_w_lane_into(&mut out.0[b], b, x, y, k, lo, hi, &w);
        }
    }
}

/// Lane vector b of `cs_range_into`, `sum_{j=lo..=hi} x[j].0[b] y[k-j].0[b]` into out: bitwise `cs_range`'s lane
/// vector b on every lane type, the f64 lanes included (their interleaved sums take each vector's terms in the order
/// j = lo..=hi too), so a caller can sum the B lane vectors on separate threads.
#[inline(always)]
pub fn cs_lane_into<T: Lane, const B: usize>(
    out: &mut T,
    b: usize,
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
) {
    let mut a = T::acc0();
    for j in lo..=hi {
        T::acc_mut(&mut a, &x[j].0[b], &y[k - j].0[b]);
    }
    *out = T::fin(a);
}

/// Lane vector b of `cs_w_into`, `sum_{j=lo..=hi} (w(j) x[j].0[b]) y[k-j].0[b]` into out: bitwise `cs_w`'s lane
/// vector b on every lane type, as `cs_lane_into`
#[allow(clippy::too_many_arguments)]
#[inline(always)]
pub fn cs_w_lane_into<T: Lane, const B: usize>(
    out: &mut T,
    b: usize,
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
    w: impl Fn(usize) -> f64,
) {
    let mut a = T::acc0();
    for j in lo..=hi {
        let xw = x[j].0[b].mulf(w(j));
        T::acc_mut(&mut a, &xw, &y[k - j].0[b]);
    }
    *out = T::fin(a);
}

/// `cs_il` into out, with the accumulators and the terms borrowed (`Lane::acc_mut`): bitwise `cs_il`
#[allow(clippy::too_many_arguments)]
#[inline(always)]
fn cs_il_into<T: Lane, const B: usize>(
    out: &mut Fl<T, B>,
    x: &[Fl<T, B>],
    y: &[Fl<T, B>],
    k: usize,
    lo: usize,
    hi: usize,
    w: impl Fn(usize) -> f64,
    weighted: bool,
) {
    let mut a = [T::acc0(); B];
    for j in lo..=hi {
        let (xj, yj) = (&x[j].0, &y[k - j].0);
        if weighted {
            let wj = w(j);
            for b in 0..B {
                let xw = xj[b].mulf(wj);
                T::acc_mut(&mut a[b], &xw, &yj[b]);
            }
        } else {
            for b in 0..B {
                T::acc_mut(&mut a[b], &xj[b], &yj[b]);
            }
        }
    }
    for b in 0..B {
        out.0[b] = T::fin(a[b]);
    }
}

/// Horner evaluation `sum_k s[k] h^k`
#[inline(always)]
pub fn horner<T: Lane, const B: usize>(s: &[Fl<T, B>], h: DD) -> Fl<T, B> {
    let mut out = Fl::<T, B>::zero();
    for k in (0..s.len()).rev() {
        out = out.muls(h) + s[k];
    }
    out
}

/// A constant (8 B) x (8 B) double-double matrix stored by columns (column j as a field over the rows), with a copy
/// of the columns' hi parts for the f64 lanes.
#[derive(Clone)]
pub struct DMat<S: Simd, const B: usize> {
    pub cols: Vec<Fl<DV<S>, B>>,
    /// `cols64[j][b] = cols[j].0[b].hi`
    pub cols64: Vec<[V8<S>; B]>,
}

impl<S: Simd, const B: usize> DMat<S, B> {
    /// from row-major (8 B, 8 B) hi / lo arrays
    pub fn from_rows(hi: &[f64], lo: &[f64]) -> Self {
        let m = 8 * B;
        let mut cols = Vec::with_capacity(m);
        for j in 0..m {
            let ch: Vec<f64> = (0..m).map(|i| hi[i * m + j]).collect();
            let cl: Vec<f64> = (0..m).map(|i| lo[i * m + j]).collect();
            cols.push(Fl::<DV<S>, B>::load(&ch, &cl));
        }
        let cols64 = cols.iter().map(|c| c.0.map(|v| v.hi)).collect();
        DMat { cols, cols64 }
    }

    /// D f in the lane type's precision (`MatLane::dmat_apply`)
    #[inline(always)]
    pub fn apply<T: MatLane<S = S>>(&self, f: &Fl<T, B>) -> Fl<T, B> {
        T::dmat_apply(self, f)
    }

    /// D f with the double-double columns, accumulated by `Lane::accc` (column j's contribution in the order j)
    #[inline(always)]
    pub fn apply_dd<T: Lane<S = S>>(&self, f: &Fl<T, B>) -> Fl<T, B> {
        let mut acc = [T::acc0(); B];
        for jb in 0..B {
            for jl in 0..8 {
                let s = f.0[jb].bcast(jl);
                let col = &self.cols[8 * jb + jl];
                for b in 0..B {
                    acc[b] = T::accc(acc[b], s, col.0[b]);
                }
            }
        }
        let mut out = [T::zero(); B];
        for b in 0..B {
            out[b] = T::fin(acc[b]);
        }
        Fl(out)
    }
}

/// A lane type with the product of a constant matrix: `DMat::apply_dd` by default (the double-double columns through
/// the dot accumulator); the f64 lanes multiply by the columns' f64 copy with fused multiply-adds, `Dtn<N>` its value
/// by the double-double columns and its tangents by the f64 copy.
pub trait MatLane: Lane {
    #[inline(always)]
    fn dmat_apply<const B: usize>(d: &DMat<Self::S, B>, f: &Fl<Self, B>) -> Fl<Self, B> {
        d.apply_dd(f)
    }
}

impl<S: Simd> MatLane for DV<S> {}
impl<S: Simd> MatLane for Du<S> {}
impl<S: Simd> MatLane for Dc<S> {}
impl<S: Simd> MatLane for QV<S> {}
impl<S: Simd> MatLane for crate::ov::OV<S> {}
impl<S: Simd> MatLane for Qu<S> {}
impl<S: Simd, const N: usize> MatLane for Jet<S, N> {}
impl<S: Simd, const N: usize> MatLane for QJet<S, N> {}

impl<S: Simd> MatLane for F8<S> {
    #[inline(always)]
    fn dmat_apply<const B: usize>(d: &DMat<S, B>, f: &Fl<Self, B>) -> Fl<Self, B> {
        let mut acc = [V8::zero(); B];
        for jb in 0..B {
            for jl in 0..8 {
                let s = f.0[jb].0.bcast(jl);
                let col = &d.cols64[8 * jb + jl];
                for b in 0..B {
                    acc[b] = s.fma(col[b], acc[b]);
                }
            }
        }
        Fl(acc.map(F8))
    }
}

impl<S: Simd> MatLane for Fu<S> {
    #[inline(always)]
    fn dmat_apply<const B: usize>(d: &DMat<S, B>, f: &Fl<Self, B>) -> Fl<Self, B> {
        let mut av = [V8::zero(); B];
        let mut ad = [V8::zero(); B];
        for jb in 0..B {
            for jl in 0..8 {
                let sv = f.0[jb].v.bcast(jl);
                let sd = f.0[jb].d.bcast(jl);
                let col = &d.cols64[8 * jb + jl];
                for b in 0..B {
                    av[b] = sv.fma(col[b], av[b]);
                    ad[b] = sd.fma(col[b], ad[b]);
                }
            }
        }
        let mut out = [Fu::zero(); B];
        for b in 0..B {
            out[b] = Fu { v: av[b], d: ad[b] };
        }
        Fl(out)
    }
}

impl<S: Simd, const N: usize> MatLane for Dtn<S, N> {
    #[inline(always)]
    fn dmat_apply<const B: usize>(dm: &DMat<S, B>, f: &Fl<Self, B>) -> Fl<Self, B> {
        let mut av = [DotAcc::zero(); B];
        let mut ad = [[V8::zero(); N]; B];
        for jb in 0..B {
            for jl in 0..8 {
                let sv = f.0[jb].v.bcast(jl);
                let sd = f.0[jb].d.map(|x| x.bcast(jl));
                let col = &dm.cols[8 * jb + jl];
                let col64 = &dm.cols64[8 * jb + jl];
                for b in 0..B {
                    av[b] = av[b].add_prod(sv, col.0[b]);
                    for i in 0..N {
                        ad[b][i] = sd[i].fma(col64[b], ad[b][i]);
                    }
                }
            }
        }
        let mut out = [Self::zero(); B];
        for b in 0..B {
            out[b] = Dtn {
                v: av[b].finish(),
                d: ad[b],
            };
        }
        Fl(out)
    }
}

impl<S: Simd, const N: usize> MatLane for Fun<S, N> {
    #[inline(always)]
    fn dmat_apply<const B: usize>(dm: &DMat<S, B>, f: &Fl<Self, B>) -> Fl<Self, B> {
        let mut av = [V8::zero(); B];
        let mut ad = [[V8::zero(); N]; B];
        for jb in 0..B {
            for jl in 0..8 {
                let sv = f.0[jb].v.bcast(jl);
                let sd = f.0[jb].d.map(|x| x.bcast(jl));
                let col64 = &dm.cols64[8 * jb + jl];
                for b in 0..B {
                    av[b] = sv.fma(col64[b], av[b]);
                    for i in 0..N {
                        ad[b][i] = sd[i].fma(col64[b], ad[b][i]);
                    }
                }
            }
        }
        let mut out = [Self::zero(); B];
        for b in 0..B {
            out[b] = Fun { v: av[b], d: ad[b] };
        }
        Fl(out)
    }
}

/// The wrapped type's product, bitwise; counted as the default product would be (8 B x B `accc` and B `fin`).
impl<L: MatLane> MatLane for Counted<L> {
    #[inline(always)]
    fn dmat_apply<const B: usize>(d: &DMat<Self::S, B>, f: &Fl<Self, B>) -> Fl<Self, B> {
        crate::dual::count_n(crate::dual::op::ACCC, (8 * B * B) as u64);
        crate::dual::count_n(crate::dual::op::FIN, B as u64);
        Fl(L::dmat_apply(d, &Fl(f.0.map(|x| x.0))).0.map(Counted))
    }
}

/// Solves of a `DDLu` of order 8 B for a field: `solve_lanes` in double-double, `solve_lanes_hi` in plain f64 with
/// the hi parts of the factors (`DDLu::solve_hi`; the lo parts out zero), `solve_lanes_refine`, the double-double
/// solve refined in the field's lane type, and `solve_lanes_jet` and `solve_lanes_jet_refine`, a jet's system solved
/// order by order in double-double and in quad-double. All panic if the order is not 8 B.
pub trait FieldSolve {
    fn solve_lanes<S: Simd, const B: usize>(&self, f: &Fl<DV<S>, B>) -> Fl<DV<S>, B>;
    fn solve_lanes_hi<S: Simd, const B: usize>(&self, f: &Fl<DV<S>, B>) -> Fl<DV<S>, B>;
    /// x with A x = rhs to the precision of the lane type T, by iterative refinement with the residual in that
    /// precision (Wilkinson; Moler, J. ACM 14 (1967) 316, doi:10.1145/321386.321394). This is the LU of A or of a
    /// double-double approximation of it; `resid(x)` returns rhs - A x formed by the caller in T's arithmetic, with
    /// its own A. x starts as the double-double solve of rhs's values (`Lane::val`), and each of the `steps` steps
    /// adds the solve of the residual's values (`Lane::from_vt` with a zero tangent, then `+`). With the LU of A to
    /// double-double, each step multiplies the error by about cond(A) 2^-104, so for quad-double lanes (`QV`) two
    /// steps reach quad-double where cond(A) is below about 10^10. For `QV` this is `DDLu::solve_refine_qd` on the
    /// field's values, bitwise.
    ///
    /// Values only: the solves and corrections carry no tangent, so a lane type with tangents (`Lane::TANGENT`: `Du`,
    /// `Fu`, `Dtn`, `Fun`) panics rather than being given a zero derivative. A jet has tangents too and panics, and
    /// `solve_lanes_jet_refine` refines its orders.
    ///
    /// ```
    /// use multidouble::dd::DD;
    /// use multidouble::field::{FieldSolve, Fl};
    /// use multidouble::linalg::DDLu;
    /// use multidouble::native::QV;
    /// let lu = DDLu::factor(8, 8, (0..64).map(|k| if k % 9 == 0 { 2.0 } else { 0.0 }).collect(), vec![0.0; 64], false)
    ///     .unwrap();
    /// let rhs = Fl::<QV, 1>([QV::splat_f(1.0)]);
    /// let x = lu.solve_lanes_refine(&rhs, 2, |x| rhs - x.muls(DD::from_f64(2.0)));
    /// assert_eq!(x.get(0).to_f64(), 0.5);
    /// ```
    ///
    /// The same with `Du` lanes (a value and a tangent) panics:
    ///
    /// ```should_panic
    /// use multidouble::dd::DD;
    /// use multidouble::dual::Lane;
    /// use multidouble::field::{FieldSolve, Fl};
    /// use multidouble::linalg::DDLu;
    /// use multidouble::native::Du;
    /// let lu = DDLu::factor(8, 8, (0..64).map(|k| if k % 9 == 0 { 2.0 } else { 0.0 }).collect(), vec![0.0; 64], false)
    ///     .unwrap();
    /// let rhs = Fl::<Du, 1>([Du::splat(DD::from_f64(1.0))]);
    /// let x = lu.solve_lanes_refine(&rhs, 2, |x| rhs - x.muls(DD::from_f64(2.0)));
    /// ```
    fn solve_lanes_refine<T: Lane, const B: usize>(
        &self,
        rhs: &Fl<T, B>,
        steps: usize,
        resid: impl FnMut(&Fl<T, B>) -> Fl<T, B>,
    ) -> Fl<T, B>;
    /// `solve_lanes_refine` for a lane type with a tangent (`Du`, `Qu`): the value and the tangent are refined
    /// together. The value solves use `self`, the tangent solves `lu_t`, the LU of the differentiated system's matrix
    /// (`self` again where that is the value system's matrix), and `resid(x)` returns rhs - A x formed by the caller
    /// in T's arithmetic for the value and for the tangent, so that the tangent's residual is that of the
    /// differentiated system (rhs_d - A_d x_v - A x_d, with the caller's A_d). The first solve of a tangent ignores
    /// A_d (its right-hand side is rhs's tangent), so a tangent takes one step more than a value to reach the same
    /// precision. Without a tangent (`Lane::TANGENT` false) it is `solve_lanes_refine`, bitwise. A complex tangent
    /// (`Dc`) is carried by its real part only (`Lane::tan`).
    ///
    /// A jet (`Lane::NJ` > 0) panics. `Lane::from_vt` sets only its orders 0 and 1, so the orders past 1 would come
    /// back zero, and `solve_lanes_jet` and `solve_lanes_jet_refine` solve a jet order by order instead.
    fn solve_lanes_refine_t<T: Lane, const B: usize>(
        &self,
        lu_t: &Self,
        rhs: &Fl<T, B>,
        steps: usize,
        resid: impl FnMut(&Fl<T, B>) -> Fl<T, B>,
    ) -> Fl<T, B>
    where
        Self: Sized;
    /// x(h) with A(h) x(h) = rhs(h) for a jet lane (`Lane::NJ` > 0), with A(h) = A_0 + sum_j A_j h^j and this the LU
    /// of A_0. The orders are solved in turn in double-double (`Lane::coef`, `Lane::set_coef`), x_k = A_0^-1 r_k with
    /// r_k = rhs_k - sum_{j=1..k} A_j x_{k-j}. The caller forms the coupling. `couple(j, x, r)` returns r - A_j x,
    /// called for j = 1, ..., k in that order with x = x_{k-j}, so the caller also fixes the order of the operations
    /// and with it the bits. A quad-double jet's coefficients are rounded to double-double here, and
    /// `solve_lanes_jet_refine` keeps them. A lane type without jet coefficients panics.
    ///
    /// ```
    /// use multidouble::dd::DD;
    /// use multidouble::dual::Lane;
    /// use multidouble::field::{FieldSolve, Fl};
    /// use multidouble::linalg::DDLu;
    /// use multidouble::native::Jet;
    /// // (2 + h) x(h) = 1, so x_k = (-1)^k / 2^(k+1)
    /// let lu = DDLu::factor(8, 8, (0..64).map(|k| if k % 9 == 0 { 2.0 } else { 0.0 }).collect(), vec![0.0; 64], false)
    ///     .unwrap();
    /// let kappa = Fl::<Jet<3>, 1>([Jet::var(DD::from_f64(2.0))]);
    /// let rhs = Fl::<Jet<3>, 1>([Jet::splat(DD::ONE)]);
    /// let x = lu.solve_lanes_jet(&rhs, |j, x, r| r - kappa.coef(j) * *x);
    /// assert_eq!([0, 1, 2].map(|k| x.coef(k).get(0).to_f64()), [0.5, -0.25, 0.125]);
    /// ```
    fn solve_lanes_jet<T: Lane, const B: usize>(
        &self,
        rhs: &Fl<T, B>,
        couple: impl FnMut(usize, &Fl<DV<T::S>, B>, Fl<DV<T::S>, B>) -> Fl<DV<T::S>, B>,
    ) -> Fl<T, B>;
    /// `solve_lanes_jet` in quad-double, for a quad-double jet (`crate::jet::QJet`). The coefficients are read and
    /// set in quad-double (`Lane::coef_qd`, `Lane::set_coef_qd`), and `couple(j, x, r)` returns r - A_j x in
    /// quad-double. Each order's solve is `solve_lanes_refine` on `QV` lanes with `steps` steps, its residual
    /// `resid0(r, x)` = r - A_0 x formed by the caller in quad-double with its own A_0 (r the order's right-hand side).
    /// When this is the LU of A_0 to double-double, two steps reach quad-double where cond(A_0) is below about 10^10,
    /// as for `solve_lanes_refine`. A double-double jet's coefficients are rounded back to double-double. A lane type
    /// without jet coefficients panics.
    ///
    /// ```
    /// use multidouble::dd::DD;
    /// use multidouble::dual::Lane;
    /// use multidouble::field::{FieldSolve, Fl};
    /// use multidouble::linalg::DDLu;
    /// use multidouble::native::QJet;
    /// use multidouble::qd::QD;
    /// // (3 + h) x(h) = 1, so x_1 = -1/9, here to quad-double
    /// let lu = DDLu::factor(8, 8, (0..64).map(|k| if k % 9 == 0 { 3.0 } else { 0.0 }).collect(), vec![0.0; 64], false)
    ///     .unwrap();
    /// let kappa = Fl::<QJet<3>, 1>([QJet::var(DD::from_f64(3.0))]);
    /// let rhs = Fl::<QJet<3>, 1>([QJet::splat(DD::ONE)]);
    /// let k0 = kappa.coef_qd(0);
    /// let x = lu.solve_lanes_jet_refine(&rhs, 2, |j, x, r| r - kappa.coef_qd(j) * *x, |r, x| *r - k0 * *x);
    /// let x1 = x.coef_qd(1).get(0);
    /// assert!((x1 + QD::ONE / QD::from_f64(9.0)).to_f64().abs() < 1e-60);
    /// ```
    fn solve_lanes_jet_refine<T: Lane, const B: usize>(
        &self,
        rhs: &Fl<T, B>,
        steps: usize,
        couple: impl FnMut(usize, &Fl<QV<T::S>, B>, Fl<QV<T::S>, B>) -> Fl<QV<T::S>, B>,
        resid0: impl FnMut(&Fl<QV<T::S>, B>, &Fl<QV<T::S>, B>) -> Fl<QV<T::S>, B>,
    ) -> Fl<T, B>;
}

impl FieldSolve for DDLu {
    fn solve_lanes<S: Simd, const B: usize>(&self, f: &Fl<DV<S>, B>) -> Fl<DV<S>, B> {
        let m = 8 * B;
        assert_eq!(
            self.n, m,
            "FieldSolve: a DDLu of order {} for a field of {m} values",
            self.n
        );
        let mut h = vec![0.0; m];
        let mut l = vec![0.0; m];
        f.store(&mut h, &mut l);
        let mut b: Vec<DD> = (0..m).map(|i| DD { hi: h[i], lo: l[i] }).collect();
        self.solve(&mut b);
        for i in 0..m {
            h[i] = b[i].hi;
            l[i] = b[i].lo;
        }
        Fl::<DV<S>, B>::load(&h, &l)
    }

    fn solve_lanes_hi<S: Simd, const B: usize>(&self, f: &Fl<DV<S>, B>) -> Fl<DV<S>, B> {
        let m = 8 * B;
        assert_eq!(
            self.n, m,
            "FieldSolve: a DDLu of order {} for a field of {m} values",
            self.n
        );
        let mut h = vec![0.0; m];
        let mut l = vec![0.0; m];
        f.store(&mut h, &mut l);
        self.solve_hi(&mut h);
        l.iter_mut().for_each(|v| *v = 0.0);
        Fl::<DV<S>, B>::load(&h, &l)
    }

    fn solve_lanes_refine<T: Lane, const B: usize>(
        &self,
        rhs: &Fl<T, B>,
        steps: usize,
        mut resid: impl FnMut(&Fl<T, B>) -> Fl<T, B>,
    ) -> Fl<T, B> {
        // a run-time check, not a const block: generic callers (a lane-generic solve that refines only for some
        // lane types) instantiate this for every lane type, so a const assertion would fail their builds
        assert!(
            !T::TANGENT,
            "FieldSolve::solve_lanes_refine refines values only: a lane type with tangents would lose them"
        );
        let z = Fl::<Dv<T>, B>::zero();
        let mut x = Fl::<T, B>::from_vt(&self.solve_lanes(&rhs.val()), &z);
        for _ in 0..steps {
            let r = resid(&x);
            x = x + Fl::<T, B>::from_vt(&self.solve_lanes(&r.val()), &z);
        }
        x
    }

    fn solve_lanes_refine_t<T: Lane, const B: usize>(
        &self,
        lu_t: &Self,
        rhs: &Fl<T, B>,
        steps: usize,
        mut resid: impl FnMut(&Fl<T, B>) -> Fl<T, B>,
    ) -> Fl<T, B> {
        // run-time, as in solve_lanes_refine
        assert!(
            T::NJ == 0,
            "FieldSolve::solve_lanes_refine_t refines one tangent: a jet would lose its orders past 1 (solve_lanes_jet)"
        );
        let z = Fl::<Dv<T>, B>::zero();
        let tsolve = |r: &Fl<T, B>| if T::TANGENT { lu_t.solve_lanes(&r.tan()) } else { z };
        let mut x = Fl::<T, B>::from_vt(&self.solve_lanes(&rhs.val()), &tsolve(rhs));
        for _ in 0..steps {
            let r = resid(&x);
            x = x + Fl::<T, B>::from_vt(&self.solve_lanes(&r.val()), &tsolve(&r));
        }
        x
    }

    fn solve_lanes_jet<T: Lane, const B: usize>(
        &self,
        rhs: &Fl<T, B>,
        mut couple: impl FnMut(usize, &Fl<DV<T::S>, B>, Fl<DV<T::S>, B>) -> Fl<DV<T::S>, B>,
    ) -> Fl<T, B> {
        // run-time, as in solve_lanes_refine
        assert!(
            T::NJ > 0,
            "FieldSolve::solve_lanes_jet solves jets: a lane type without jet coefficients has no orders"
        );
        let mut x = Fl::<T, B>::zero();
        let mut xs: Vec<Fl<DV<T::S>, B>> = Vec::with_capacity(T::NJ);
        for k in 0..T::NJ {
            let mut r = rhs.coef(k);
            for j in 1..=k {
                r = couple(j, &xs[k - j], r);
            }
            let xk = self.solve_lanes(&r);
            x.set_coef(k, &xk);
            xs.push(xk);
        }
        x
    }

    fn solve_lanes_jet_refine<T: Lane, const B: usize>(
        &self,
        rhs: &Fl<T, B>,
        steps: usize,
        mut couple: impl FnMut(usize, &Fl<QV<T::S>, B>, Fl<QV<T::S>, B>) -> Fl<QV<T::S>, B>,
        mut resid0: impl FnMut(&Fl<QV<T::S>, B>, &Fl<QV<T::S>, B>) -> Fl<QV<T::S>, B>,
    ) -> Fl<T, B> {
        assert!(
            T::NJ > 0,
            "FieldSolve::solve_lanes_jet_refine solves jets: a lane type without jet coefficients has no orders"
        );
        let mut x = Fl::<T, B>::zero();
        let mut xs: Vec<Fl<QV<T::S>, 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 = couple(j, &xs[k - j], r);
            }
            let xk = self.solve_lanes_refine(&r, steps, |y| resid0(&r, y));
            x.set_coef_qd(k, &xk);
            xs.push(xk);
        }
        x
    }
}

// ---------------------------------------------------------------------------------------------------------------
// early end of a series for the tangents
// ---------------------------------------------------------------------------------------------------------------

/// Early end of a power series for the tangents of a dual lane type, when they are needed to less precision than the
/// values whose accuracy sets the series' order (the f64 tangents of a Jacobian in f64, or double-double tangents
/// that would otherwise run to the values' order). Call `done(k, &[c_k of each field])` for k = 0, 1, 2, ... in
/// turn, with the coefficients of order k; it returns true once the tangent terms max |c_k| r^k (r the evaluation
/// radius; the maximum over the fields, lanes and tangents given, f64 tangents through `Lane::tan64`, else the hi
/// part of `Lane::tan`, so a complex tangent (`Dc`) by its real part) of `run` consecutive orders are all below
/// `tol` times the largest tangent term so far, and k >= kmin. The values are not looked at, so a value lane type
/// never ends early.
pub struct TangentTail {
    rk: f64,
    r: f64,
    big: f64,
    run: usize,
    tol: f64,
    run_len: usize,
    kmin: usize,
}

impl TangentTail {
    pub fn new(r: f64, tol: f64, run: usize, kmin: usize) -> Self {
        TangentTail {
            rk: 1.0,
            r,
            big: 0.0,
            run: 0,
            tol,
            run_len: run,
            kmin,
        }
    }

    /// record the coefficients of order k; true when the series may end after it
    pub fn done<T: Lane, const B: usize>(&mut self, k: usize, f: &[&Fl<T, B>]) -> bool {
        let mut a = 0.0f64;
        for x in f {
            for b in 0..B {
                if T::TAN_F64 {
                    for t in 0..T::NT {
                        a = a.max(x.0[b].tan64(t).abs().hmax());
                    }
                } else {
                    a = a.max(x.0[b].tan().hi.abs().hmax());
                }
            }
        }
        let t = a * self.rk;
        self.rk *= self.r;
        self.big = self.big.max(t);
        self.run = if t <= self.tol * self.big { self.run + 1 } else { 0 };
        self.big > 0.0 && k >= self.kmin && self.run >= self.run_len
    }
}