Skip to main content

quadratic_roots_batch

Function quadratic_roots_batch 

Source
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 — one QuadraticEqn per 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);