pub fn quadratic_roots_batch(
eqns: &[QuadraticEqn],
backend: ComputeBackend,
) -> Vec<Roots<2>>Expand description
Roots of N independent quadratic equations a x^2 + b x + c = 0, on the
chosen backend.
A batched wrapper over QuadraticEqn::roots, which uses a
Kahan-compensated discriminant; the per-equation solver is unchanged and is
not reimplemented here.
§Arguments
eqns— oneQuadraticEqnper lane, dimensionless coefficients.backend— requested execution backend.
§Returns
One Roots<2> per equation, in input order. Two Real roots when the
discriminant is positive; a Complex-tagged (real part, imaginary part) pair
when it is negative; one Real plus one Nan when |a| underflows and the
equation is really linear.
§Determinism
Bit-for-bit identical to calling QuadraticEqn::roots in a loop, on every
backend and at any thread count.
§Verification
Methodology. Every lane’s real roots are substituted back into their own
equation and the residual |a x^2 + b x + c| compared against zero, over
4 096 equations built as (x - r1)(x - r2) with pseudorandom real roots on
[-8, 8) from a fixed-seed xorshift64*; pass criterion: residual
<= 1e-9 and bitwise equality against the per-equation scalar solver.
Result, measured 2026-08-12 by quadratic_batch_residuals_are_analytic in
parallel/tests.rs, release build: worst residual 7.105427e-15 over
8 192 real roots; bitwise equality against the scalar solver held on every
lane. Interpretation: the batching adds nothing to the scalar solver’s own
error, which is what “bit-for-bit identical” means measured rather than
argued.
§Example
use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::math::parallel::quadratic_roots_batch;
use outram_foam_basic_lib::polynomial::{QuadraticEqn, RootType};
// (x - 1)(x - 2) = x^2 - 3x + 2
let eqns = [QuadraticEqn::new(1.0, -3.0, 2.0)];
let roots = quadratic_roots_batch(&eqns, ComputeBackend::CpuMulti);
let mut real: Vec<f64> = (0..2)
.filter(|&k| roots[0].root_type(k) == RootType::Real)
.map(|k| roots[0].get(k))
.collect();
real.sort_by(f64::total_cmp);
assert!((real[0] - 1.0).abs() < 1e-12);
assert!((real[1] - 2.0).abs() < 1e-12);