Skip to main content

cubic_roots_batch

outram_foam_basic_lib::math::parallel

Function cubic_roots_batch 

Source 
pub fn cubic_roots_batch(
    eqns: &[CubicEqn],
    backend: ComputeBackend,
) -> Vec<Roots<3>>
Expand description

Roots of N independent cubic equations a x^3 + b x^2 + c x + d = 0, on the chosen backend.

A batched wrapper over CubicEqn::roots — the depressed-cubic Cardano solver with Kahan-compensated discriminants, unchanged and not reimplemented here.

This is the kernel a cubic-equation-of-state caller wants: a Peng-Robinson or Redlich-Kwong compressibility factor is the root of a cubic in Z, one cubic per cell, and the whole field can be solved in a single call.

§Arguments

§Returns

One Roots<3> per equation, in input order. Three Real roots, or one Real plus a Complex-tagged (real, imaginary) pair, or a triple root, or — when |a| underflows — the quadratic’s roots plus a Nan slot. The tags are the caller’s guide; do not read a slot without checking Roots::root_type.

There is no _into variant, unlike crate::ldu_matrix::parallel’s kernels: Roots has no Default, so a caller-owned output buffer would have to be pre-filled with a dummy value, and the per-lane Cardano solve dominates the allocation anyway.

§Determinism

Bit-for-bit identical to calling CubicEqn::roots in a loop, on every backend and at any thread count. Each lane is one independent closed-form evaluation; no arithmetic crosses lanes.

§Verification

Methodology. Cubics are constructed from known roots as (x - r1)(x - r2)(x - r3) with r1, r2, r3 drawn from a fixed-seed xorshift64* on [-8, 8), expanded to coefficients — so the analytic answer is known exactly, which makes this an oracle rather than a second opinion. 4 096 equations. Two pass criteria: every Real-tagged root satisfies |a x^3 + b x^2 + c x + d| <= 1e-6 (loose because a cubic’s residual is steeply amplified near a triple-ish root), and every returned bit is equal to the per-equation scalar CubicEqn::roots.

Result, measured 2026-08-12 by cubic_batch_matches_analytic_construction in parallel/tests.rs, release build: worst residual 3.979039e-13 over 12 288 real roots, worst distance from the nearest constructed root 2.079981e-10; bitwise equality against the scalar solver held on every lane, and (per poly_batch_bitwise_identical_across_backends) across Serial / CpuMulti at 1, 2, 4 and 8 worker threads.

Interpretation: both figures belong to the scalar Cardano solver, not to the batching — the bitwise-equality half of the test proves the batch adds exactly zero error. That the worst root displacement (2.1e-10) is three orders larger than the worst residual (4.0e-13) is the signature of an ill-conditioned root, not of a bad solve: with roots drawn at random from [-8, 8), some triples come out nearly coincident, and near a near-double root the cubic is flat, so a tiny residual still corresponds to a visibly displaced abscissa. A caller solving a cubic equation of state should expect the same behaviour near the critical point and judge accuracy by the residual.

§Example

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::math::parallel::cubic_roots_batch;
use outram_foam_basic_lib::polynomial::{CubicEqn, RootType};

// (x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x - 6
let eqns = [CubicEqn::new(1.0, -6.0, 11.0, -6.0)];
let roots = cubic_roots_batch(&eqns, ComputeBackend::CpuMulti);

let mut real: Vec<f64> = (0..3)
    .filter(|&k| roots[0].root_type(k) == RootType::Real)
    .map(|k| roots[0].get(k))
    .collect();
real.sort_by(f64::total_cmp);
assert_eq!(real.len(), 3);
assert!((real[0] - 1.0).abs() < 1e-9);
assert!((real[2] - 3.0).abs() < 1e-9);