rustnum 0.1.0

//! Polynomial tests: NumPy-computed constants plus mathematical properties.
// @generated fixtures: values computed with NumPy 2.3.2.
const PV_P: [f64; 4] = [2.0, -3.0, 0.5, 7.0];
const PV_X: [f64; 6] = [-2.5, -1.0, 0.0, 0.3, 1.0, 4.75];
const PV_Y: [f64; 6] = [-44.25, 1.5, 7.0, 6.934, 6.5, 156.03125];
const A: [f64; 4] = [1.0, -2.0, 3.0, 4.0];
const B: [f64; 2] = [5.0, 0.5];
const ADD_AB: [f64; 4] = [1.0, -2.0, 8.0, 4.5];
const SUB_AB: [f64; 4] = [1.0, -2.0, -2.0, 3.5];
const SUB_BA: [f64; 4] = [-1.0, 2.0, 2.0, -3.5];
const ADD_CANCEL: [f64; 2] = [0.0, 5.0];
const MUL_AB: [f64; 5] = [5.0, -9.5, 14.0, 21.5, 2.0];
const DIV_AB_Q: [f64; 3] = [0.2, -0.42000000000000004, 0.642];
const DIV_AB_R: [f64; 1] = [3.679];
const U: [f64; 6] = [1.0, 0.0, -3.0, 0.0, 2.0, 5.0];
const V: [f64; 3] = [2.0, 0.0, 1.0];
const DIV_UV_Q: [f64; 4] = [0.5, 0.0, -1.75, 0.0];
const DIV_UV_R: [f64; 2] = [3.75, 5.0];
const DIV_EXACT_Q: [f64; 2] = [1.0, -2.0];
const DIV_EXACT_R: [f64; 1] = [0.0];
const DIV_SMALL_Q: [f64; 1] = [0.0];
const DIV_SMALL_R: [f64; 2] = [3.0, 1.0];
const DER1: [f64; 3] = [3.0, -4.0, 3.0];
const DER2: [f64; 2] = [6.0, -4.0];
const DER4: [f64; 0] = [];
const INT1: [f64; 5] = [0.25, -0.6666666666666666, 1.5, 4.0, 0.0];
const INT2K: [f64; 6] = [0.05, -0.16666666666666666, 0.5, 2.0, 1.0, 2.0];
const INT3K1: [f64; 7] = [
    0.008333333333333333,
    -0.03333333333333333,
    0.125,
    0.6666666666666666,
    2.0,
    4.0,
    4.0,
];
const FIT_X: [f64; 21] = [
    0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
    1.9, 2.0,
];
const FIT_Y: [f64; 21] = [
    0.0,
    0.4955202066613396,
    0.9646424733950355,
    0.8833269096274833,
    1.2320390859672266,
    0.9974949866040544,
    1.1738476308781953,
    1.263209366648874,
    0.7754631805511506,
    0.7273798802338298,
    0.1411200080598672,
    0.042254305856751356,
    -0.04252044329485205,
    -0.5877661591839741,
    -0.5715757724135877,
    -0.977530117665097,
    -0.7961646088358405,
    -0.5258146823277324,
    -0.6727644875559872,
    -0.2506855425976383,
    -0.27941549819892586,
];
const FIT_D1: [f64; 2] = [-0.9009521332404372, 1.091050262879017];
const FIT_D2: [f64; 3] = [
    -0.3685690435813028,
    -0.16381404607783098,
    0.8576232019441914,
];
const FIT_D3: [f64; 4] = [
    2.5133647868452136,
    -7.9086634041169335,
    5.722486284713642,
    -0.0019475551568655875,
];
const FIT_D5: [f64; 6] = [
    -1.31938326272825,
    6.460007713558401,
    -8.542243617816832,
    -0.17994451126405064,
    3.8201444313868405,
    0.07271815198620898,
];
const FIT_D8: [f64; 9] = [
    -3.616552083295404,
    29.76582142673909,
    -102.08832696641937,
    187.21922472862386,
    -194.7686658626967,
    114.16981752329677,
    -39.11182652692485,
    8.789370449009056,
    -0.014867090309481513,
];
const FIT2_X: [f64; 15] = [
    1000.0, 1001.0, 1002.0, 1003.0, 1004.0, 1005.0, 1006.0, 1007.0, 1008.0, 1009.0, 1010.0, 1011.0,
    1012.0, 1013.0, 1014.0,
];
const FIT2_Y: [f64; 15] = [
    3.025, 3.046, 3.029, 3.014, 3.001, 3.03, 3.021, 3.014, 3.009, 3.046, 3.045, 3.046, 3.049,
    3.094, 3.101,
];
const FIT2_D2: [f64; 3] = [
    0.0010258564964420317,
    -2.0617892695485227,
    1038.9718907535496,
];
const R_123_P: [f64; 4] = [1.0, -6.0, 11.0, -6.0];
const R_123: [(f64, f64); 3] = [
    (3.0000000000000018, 0.0),
    (1.999999999999998, 0.0),
    (1.0000000000000002, 0.0),
];
const R_I_P: [f64; 3] = [1.0, 0.0, 1.0];
const R_I: [(f64, f64); 2] = [(-0.0, 1.0), (0.0, -1.0)];
const R_CUBIC_P: [f64; 4] = [1.0, -6.0, 11.0, -6.0];
const R_CUBIC: [(f64, f64); 3] = [
    (3.0000000000000018, 0.0),
    (1.999999999999998, 0.0),
    (1.0000000000000002, 0.0),
];
const R_QUARTIC_P: [f64; 5] = [1.0, 2.0, 3.0, 4.0, 5.0];
const R_QUARTIC: [(f64, f64); 4] = [
    (-1.287815479557648, 0.8578967583284908),
    (-1.287815479557648, -0.8578967583284908),
    (0.28781547955764797, 1.4160930801719085),
    (0.28781547955764797, -1.4160930801719085),
];
const R_DEG6_P: [f64; 7] = [3.0, 0.0, -2.0, 5.0, 1.0, -7.0, 2.0];
const R_DEG6: [(f64, f64); 6] = [
    (-1.1620299333554476, 0.4396771928850095),
    (-1.1620299333554476, -0.4396771928850095),
    (0.5722553754616273, 1.1124339542978319),
    (0.5722553754616273, -1.1124339542978319),
    (0.8578569955636121, 0.0),
    (0.3216921202240273, 0.0),
];
const R_ZEROS_P: [f64; 7] = [0.0, 0.0, 1.0, -3.0, 2.0, 0.0, 0.0];
const R_ZEROS: [(f64, f64); 4] = [(2.0, 0.0), (1.0, 0.0), (0.0, 0.0), (0.0, 0.0)];
const R_UNITY5_P: [f64; 6] = [1.0, 0.0, 0.0, 0.0, 0.0, -1.0];
const R_UNITY5: [(f64, f64); 5] = [
    (-0.8090169943749476, 0.5877852522924731),
    (-0.8090169943749476, -0.5877852522924731),
    (0.30901699437494734, 0.9510565162951533),
    (0.30901699437494734, -0.9510565162951533),
    (0.9999999999999997, 0.0),
];
const R_LINEAR_P: [f64; 2] = [2.0, -4.0];
const R_LINEAR: [(f64, f64); 1] = [(2.0, 0.0)];
const R_WILK10_P: [f64; 11] = [
    1.0,
    -55.0,
    1320.0,
    -18150.0,
    157773.0,
    -902055.0,
    3416930.0,
    -8409500.0,
    12753576.0,
    -10628640.0,
    3628800.0,
];
const R_WILK10: [(f64, f64); 10] = [
    (9.99999999976065, 0.0),
    (9.000000001383428, 0.0),
    (7.9999999967098, 0.0),
    (7.000000004192445, 0.0),
    (5.999999996891773, 0.0),
    (5.000000001359408, 0.0),
    (3.9999999996602487, 0.0),
    (3.0000000000448317, 0.0),
    (1.9999999999973728, 0.0),
    (1.000000000000033, 0.0),
];
const R_MIXED_P: [f64; 8] = [1.0, -1.5, -3.25, 7.125, -17.375, -43.875, -11.875, 18.75];
const R_MIXED: [(f64, f64); 7] = [
    (3.0000000000000018, 0.0),
    (0.9999999999999984, 1.9999999999999978),
    (0.9999999999999984, -1.9999999999999978),
    (-2.0000000000000027, 0.0),
    (-1.0000000000000007, 0.5),
    (-1.0000000000000007, -0.5),
    (0.5000000000000002, 0.0),
];
const R_DEG12_P: [f64; 13] = [
    1.0, -0.5, 2.0, 0.0, -3.0, 1.0, 4.0, -2.0, 0.25, 1.0, -1.0, 2.0, 5.0,
];
const R_DEG12: [(f64, f64); 12] = [
    (0.11359091342389227, 1.7815027535626269),
    (0.11359091342389227, -1.7815027535626269),
    (1.0929335021475517, 0.36563863377332934),
    (1.0929335021475517, -0.36563863377332934),
    (0.7917779380470269, 0.7742975348611193),
    (0.7917779380470269, -0.7742975348611193),
    (0.02470026132855882, 0.9583142109689866),
    (0.02470026132855882, -0.9583142109689866),
    (-0.9563122661959573, 0.3402127430628228),
    (-0.9563122661959573, -0.3402127430628228),
    (-0.816690348751071, 0.5919071182190434),
    (-0.816690348751071, -0.5919071182190434),
];
const R_RAND0_P: [f64; 13] = [
    0.2987455375084699,
    -0.2741378553622176,
    -0.8905918387572742,
    -0.45467078517172255,
    -0.9916465549964624,
    0.060143602597438485,
    1.3402152455545335,
    -0.49220651855132963,
    -0.6204748998199404,
    0.4898420501851982,
    0.35688700816006075,
    0.10541424899789856,
    -0.9304680447082047,
];
const R_RAND0: [(f64, f64); 12] = [
    (2.5129080853978163, 0.0),
    (-1.2439990657833135, 0.0),
    (-1.0864617597787305, 0.38727356908454436),
    (-1.0864617597787305, -0.38727356908454436),
    (-0.0496751036348092, 1.2899204385747876),
    (-0.0496751036348092, -1.2899204385747876),
    (-0.4291421378266405, 0.7119807587619146),
    (-0.4291421378266405, -0.7119807587619146),
    (0.5315749972644248, 0.7289606800728178),
    (0.5315749972644248, -0.7289606800728178),
    (0.858064473713503, 0.25038659653139456),
    (0.858064473713503, -0.25038659653139456),
];
const R_RAND1_P: [f64; 10] = [
    -0.02925182246327349,
    0.0006953031944582878,
    -0.0013442145472850819,
    -0.45761576104021817,
    -1901.222739800844,
    -1.289537739784976,
    -184.17350377917325,
    -235.09113107468127,
    -126.74464814437032,
    2.7126435882170155,
];
const R_RAND1: [(f64, f64); 9] = [
    (11.296966152928595, 11.304939557730288),
    (11.296966152928595, -11.304939557730288),
    (-11.284753727071724, 11.27449696160991),
    (-11.284753727071724, -11.27449696160991),
    (0.3330983505446268, 0.5051056696222311),
    (0.3330983505446268, -0.5051056696222311),
    (-0.3437271748915496, 0.26650490887617434),
    (-0.3437271748915496, -0.26650490887617434),
    (0.020602365366185614, 0.0),
];
const R_RAND2_P: [f64; 8] = [
    -1.5301357655053935,
    -0.47775327603393064,
    -0.9785190780566395,
    -0.8088372394255993,
    1.0608986233860787,
    -0.8075346753318965,
    -0.0325217049455206,
    0.8843898673831739,
];
const R_RAND2: [(f64, f64); 7] = [
    (-0.13857631887736482, 1.1562947394649767),
    (-0.13857631887736482, -1.1562947394649767),
    (-0.8814657933691092, 0.0),
    (-0.8480321826588088, 0.0),
    (0.46488446648749043, 0.7276502821171952),
    (0.46488446648749043, -0.7276502821171952),
    (0.7646523435134465, 0.0),
];
const R_RAND3_P: [f64; 9] = [
    -58.360043274330195,
    -0.11170194958415963,
    1.104641432494806,
    63.781774255061954,
    -12.250558264176934,
    0.00761402303770081,
    13.588234217415376,
    -15.471446781284824,
    0.0008593826880215983,
];
const R_RAND3: [(f64, f64); 8] = [
    (-0.6262493124664997, 0.9495907284966595),
    (-0.6262493124664997, -0.9495907284966595),
    (-0.706075781436964, 0.0),
    (0.16649856594027074, 0.6507244000872888),
    (0.16649856594027074, -0.6507244000872888),
    (0.937850242774563, 0.0),
    (0.6857574684695118, 0.0),
    (5.554907885401621e-05, 0.0),
];
const R_RAND4_P: [f64; 5] = [
    0.07451622877146342,
    0.5766895836701853,
    -0.1887821253507493,
    0.682910267195206,
    -0.06651732014941557,
];
const R_RAND4: [(f64, f64); 4] = [
    (-8.186922260554272, 0.0),
    (0.17425876548165797, 1.033327124670754),
    (0.17425876548165797, -1.033327124670754),
    (0.09929080530935885, 0.0),
];
const R_RAND5_P: [f64; 10] = [
    1438.522591656152,
    -675.6622510056528,
    0.0020313861038960906,
    -0.00046330757653841514,
    1.2726841122583084,
    -0.0011871945278501398,
    -0.05793015965026732,
    -0.1961959728044967,
    898.7638721004078,
    0.001145222007454132,
];
const R_RAND5: [(f64, f64); 9] = [
    (0.9441294267057248, 0.352830584638361),
    (0.9441294267057248, -0.352830584638361),
    (0.42162637261509106, 0.8574666237287917),
    (0.42162637261509106, -0.8574666237287917),
    (-0.3086644835736533, 0.8616893184714663),
    (-0.3086644835736533, -0.8616893184714663),
    (-0.822244783646267, 0.35771319464399426),
    (-0.822244783646267, -0.35771319464399426),
    (-1.274219005331402e-06, 0.0),
];
const R_RAND6_P: [f64; 10] = [
    -0.09728692567008902,
    1.2570149772868198,
    0.6894039005707556,
    -0.32721342022219785,
    -0.3685758940999591,
    -0.25019540051792494,
    1.5235294004561601,
    -0.4280249425728672,
    -0.3036803883647294,
    0.35258906728526535,
];
const R_RAND6: [(f64, f64); 9] = [
    (13.42815562441223, 0.0),
    (-1.1762255124082708, 0.0),
    (-0.49890758150626746, 0.9956287634366944),
    (-0.49890758150626746, -0.9956287634366944),
    (-0.6179921430279388, 0.0),
    (0.7061445924350094, 0.4153344883186751),
    (0.7061445924350094, -0.4153344883186751),
    (0.4361429601535919, 0.5058356006589692),
    (0.4361429601535919, -0.5058356006589692),
];
const R_RAND7_P: [f64; 8] = [
    -197.28422796572255,
    -111.40671431510563,
    -0.011521468038548173,
    -0.4435812229744192,
    0.11661277761902228,
    0.6530885027011638,
    -2.414361300993223,
    0.0006683810232673439,
];
const R_RAND7: [(f64, f64); 7] = [
    (-0.564797008080977, 0.184583018250353),
    (-0.564797008080977, -0.184583018250353),
    (-0.0642793522325968, 0.4480433222074588),
    (-0.0642793522325968, -0.4480433222074588),
    (0.34658714452177897, 0.22147446603692048),
    (0.34658714452177897, -0.22147446603692048),
    (0.0002768562784308821, 0.0),
];
type RootCase = (&'static [f64], &'static [(f64, f64)]);
const ROOT_CASES: [RootCase; 19] = [
    (&R_123_P, &R_123),
    (&R_I_P, &R_I),
    (&R_CUBIC_P, &R_CUBIC),
    (&R_QUARTIC_P, &R_QUARTIC),
    (&R_DEG6_P, &R_DEG6),
    (&R_ZEROS_P, &R_ZEROS),
    (&R_UNITY5_P, &R_UNITY5),
    (&R_LINEAR_P, &R_LINEAR),
    (&R_WILK10_P, &R_WILK10),
    (&R_MIXED_P, &R_MIXED),
    (&R_DEG12_P, &R_DEG12),
    (&R_RAND0_P, &R_RAND0),
    (&R_RAND1_P, &R_RAND1),
    (&R_RAND2_P, &R_RAND2),
    (&R_RAND3_P, &R_RAND3),
    (&R_RAND4_P, &R_RAND4),
    (&R_RAND5_P, &R_RAND5),
    (&R_RAND6_P, &R_RAND6),
    (&R_RAND7_P, &R_RAND7),
];

use super::kernels::*;

fn close(got: &[f64], want: &[f64], tol: f64, what: &str) {
    assert_eq!(got.len(), want.len(), "{what}: length {got:?} vs {want:?}");
    let scale = 1.0 + want.iter().fold(0.0f64, |m, v| m.max(v.abs()));
    for (i, (g, w)) in got.iter().zip(want).enumerate() {
        assert!(
            (g - w).abs() <= tol * scale,
            "{what}[{i}]: got {g:e}, want {w:e}\n got  {got:?}\n want {want:?}"
        );
    }
}

/// Every wanted root has a distinct computed root within `tol` (order ignored).
fn same_roots(got: &[(f64, f64)], want: &[(f64, f64)], tol: f64, what: &str) {
    assert_eq!(got.len(), want.len(), "{what}: {got:?} vs {want:?}");
    let mut used = vec![false; got.len()];
    for w in want {
        let (best, d) = got
            .iter()
            .enumerate()
            .filter(|(i, _)| !used[*i])
            .map(|(i, g)| (i, (g.0 - w.0).hypot(g.1 - w.1)))
            .fold(
                (usize::MAX, f64::INFINITY),
                |a, b| if b.1 < a.1 { b } else { a },
            );
        let scale = 1.0 + w.0.hypot(w.1);
        assert!(
            d <= tol * scale,
            "{what}: no match for {w:?} (nearest off by {d:e}) in {got:?}"
        );
        used[best] = true;
    }
}

// ------------------------------------------------------------------ NumPy parity

#[test]
fn polyval_matches_numpy() {
    let got: Vec<f64> = PV_X.iter().map(|&x| polyval(&PV_P, x)).collect();
    close(&got, &PV_Y, 1e-14, "polyval");
    assert_eq!(polyval(&[], 3.0), 0.0);
    assert_eq!(polyval(&[4.0], 3.0), 4.0);
}

#[test]
fn add_sub_mul_match_numpy() {
    close(&polyadd(&A, &B), &ADD_AB, 1e-15, "polyadd");
    close(&polysub(&A, &B), &SUB_AB, 1e-15, "polysub");
    close(&polysub(&B, &A), &SUB_BA, 1e-15, "polysub reversed");
    close(
        &polyadd(&[1.0, 2.0], &[-1.0, 3.0]),
        &ADD_CANCEL,
        0.0,
        "polyadd keeps leading zero",
    );
    close(&polymul(&A, &B), &MUL_AB, 1e-14, "polymul");
}

#[test]
fn polydiv_matches_numpy() {
    let (q, r) = polydiv(&A, &B);
    close(&q, &DIV_AB_Q, 1e-14, "quotient");
    close(&r, &DIV_AB_R, 1e-14, "remainder");
    let (q, r) = polydiv(&U, &V);
    close(&q, &DIV_UV_Q, 1e-14, "quotient 2");
    close(&r, &DIV_UV_R, 1e-14, "remainder 2");
    let (q, r) = polydiv(&[1.0, -3.0, 2.0], &[1.0, -1.0]);
    close(&q, &DIV_EXACT_Q, 1e-15, "exact quotient");
    close(&r, &DIV_EXACT_R, 1e-15, "exact remainder");
    let (q, r) = polydiv(&[3.0, 1.0], &[1.0, 0.0, 0.0]);
    close(&q, &DIV_SMALL_Q, 0.0, "deg(u) < deg(v) quotient");
    close(&r, &DIV_SMALL_R, 0.0, "deg(u) < deg(v) remainder");
}

#[test]
fn polyder_and_polyint_match_numpy() {
    close(&polyder(&A, 1), &DER1, 1e-15, "polyder");
    close(&polyder(&A, 2), &DER2, 1e-15, "polyder m=2");
    assert!(polyder(&A, 4).is_empty() && DER4.is_empty());
    assert_eq!(polyder(&A, 0), A);
    close(&polyint(&A, 1, &[]), &INT1, 1e-15, "polyint");
    close(&polyint(&A, 2, &[1.0, 2.0]), &INT2K, 1e-15, "polyint m=2 k");
    close(
        &polyint(&A, 3, &[4.0]),
        &INT3K1,
        1e-15,
        "polyint broadcast k",
    );
    assert_eq!(polyint(&A, 0, &[]), A);
}

#[test]
fn polyfit_matches_numpy() {
    let fits: [(usize, &[f64]); 5] = [
        (1, &FIT_D1),
        (2, &FIT_D2),
        (3, &FIT_D3),
        (5, &FIT_D5),
        (8, &FIT_D8),
    ];
    for (d, want) in fits {
        let got = polyfit(&FIT_X, &FIT_Y, d);
        // Compare in the scaled sense: coefficient c_j multiplies x^(d-j) with x in [0, 2].
        for (j, (g, w)) in got.iter().zip(want).enumerate() {
            let tol = 1e-8 * (1.0 + w.abs());
            assert!(
                (g - w).abs() < tol,
                "deg {d} coef {j}: got {g:e}, want {w:e}"
            );
        }
    }
}

#[test]
fn polyfit_handles_badly_scaled_abscissae() {
    let got = polyfit(&FIT2_X, &FIT2_Y, 2);
    // The fitted curve (not the raw coefficients, which cancel heavily) must agree with NumPy's.
    for &x in &FIT2_X {
        let g = polyval(&got, x);
        let w = polyval(&FIT2_D2, x);
        assert!((g - w).abs() < 1e-9, "x={x}: {g} vs {w}");
    }
}

#[test]
fn roots_match_numpy_as_sets() {
    for (i, (p, want)) in ROOT_CASES.iter().enumerate() {
        let got = roots(p);
        // Wilkinson-10 is ill-conditioned: NumPy itself is only good to ~1e-6.
        let tol = if i == 8 { 1e-5 } else { 1e-8 };
        same_roots(&got, want, tol, &format!("case {i} {p:?}"));
    }
}

#[test]
fn roots_order_matches_numpy_on_simple_cases() {
    // (x-1)(x-2)(x-3) -> [3, 2, 1] exactly like np.roots.
    let r = roots(&[1.0, -6.0, 11.0, -6.0]);
    close(
        &r.iter().map(|c| c.0).collect::<Vec<_>>(),
        &[3.0, 2.0, 1.0],
        1e-12,
        "re",
    );
    assert!(r.iter().all(|c| c.1 == 0.0));
    // Trailing zero roots are appended last.
    let r = roots(&[0.0, 0.0, 1.0, -3.0, 2.0, 0.0, 0.0]);
    assert_eq!(r.len(), 4);
    assert!((r[0].0 - 2.0).abs() < 1e-12 && (r[1].0 - 1.0).abs() < 1e-12);
    assert_eq!((r[2], r[3]), ((0.0, 0.0), (0.0, 0.0)));
    // Conjugate pairs: positive imaginary part first.
    let r = roots(&[1.0, 0.0, 1.0]);
    assert!(r[0].1 > 0.0 && r[1].1 < 0.0 && (r[0].1 + r[1].1).abs() < 1e-15);
    // Every fixture case (incl. random and badly scaled ones) comes out in NumPy's order.
    for (i, (p, want)) in ROOT_CASES.iter().enumerate() {
        let got = roots(p);
        assert_eq!(got.len(), want.len(), "case {i}");
        for (k, (g, w)) in got.iter().zip(want.iter()).enumerate() {
            let tol = if i == 8 { 1e-5 } else { 1e-7 };
            assert!(
                (g.0 - w.0).hypot(g.1 - w.1) < tol * (1.0 + w.0.hypot(w.1)),
                "case {i} root {k}: got {g:?}, NumPy {w:?}"
            );
        }
    }
}

#[test]
fn degenerate_root_inputs() {
    assert!(roots(&[]).is_empty());
    assert!(roots(&[5.0]).is_empty());
    assert!(roots(&[0.0, 0.0]).is_empty());
    assert_eq!(roots(&[2.0, -4.0]), vec![(2.0, 0.0)]);
    assert_eq!(roots(&[1.0, 0.0, 0.0]), vec![(0.0, 0.0), (0.0, 0.0)]);
}

// ------------------------------------------------------------------ properties

#[test]
fn polyfit_recovers_exact_polynomials() {
    let truth: [&[f64]; 5] = [
        &[2.5],
        &[-1.0, 4.0],
        &[3.0, -2.0, 0.5],
        &[0.5, 1.0, -2.0, 3.0],
        &[1.0, -3.0, 0.25, 2.0, -1.0, 0.5],
    ];
    for coef in truth {
        let deg = coef.len() - 1;
        let x: Vec<f64> = (0..30).map(|k| -1.5 + 0.17 * k as f64).collect();
        let y: Vec<f64> = x.iter().map(|&t| polyval(coef, t)).collect();
        let fit = polyfit(&x, &y, deg);
        close(&fit, coef, 1e-9, &format!("coefficients deg {deg}"));
        for (&xi, &yi) in x.iter().zip(&y) {
            assert!((polyval(&fit, xi) - yi).abs() < 1e-9 * (1.0 + yi.abs()));
        }
    }
}

#[test]
fn polyfit_through_exactly_determined_points_interpolates() {
    let x = [0.0, 1.0, 2.0, 3.0];
    let y = [1.0, -1.0, 4.0, 0.5];
    let fit = polyfit(&x, &y, 3);
    for (&xi, &yi) in x.iter().zip(&y) {
        assert!((polyval(&fit, xi) - yi).abs() < 1e-10);
    }
}

#[test]
fn polyfit_rank_deficient_still_passes_through_the_points() {
    // Two points, quadratic: underdetermined; any exact solution interpolates.
    let x = [1.0, 2.0];
    let y = [3.0, 5.0];
    let fit = polyfit(&x, &y, 2);
    assert_eq!(fit.len(), 3);
    for (&xi, &yi) in x.iter().zip(&y) {
        assert!((polyval(&fit, xi) - yi).abs() < 1e-10);
    }
    // Repeated abscissae.
    let fit = polyfit(&[2.0, 2.0, 2.0], &[1.0, 1.0, 1.0], 1);
    assert!((polyval(&fit, 2.0) - 1.0).abs() < 1e-10);
}

#[test]
fn roots_of_products_of_linear_factors() {
    let sets: [&[f64]; 4] = [
        &[1.0, 2.0, 3.0],
        &[-4.0, 0.5, 2.0, 7.0],
        &[0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
        &[-3.0, -2.0, -1.0, 1.0, 2.0, 3.0, 5.0],
    ];
    for rs in sets {
        let mut p = vec![1.0];
        for &r in rs {
            p = polymul(&p, &[1.0, -r]);
        }
        let got = roots(&p);
        let want: Vec<(f64, f64)> = rs.iter().map(|&r| (r, 0.0)).collect();
        same_roots(&got, &want, 1e-9, &format!("{rs:?}"));
    }
}

#[test]
fn roots_are_zeros_of_the_polynomial() {
    let p: [f64; 8] = [1.0, -0.5, 2.0, 0.0, -3.0, 1.0, 4.0, -2.0];
    let scale: f64 = p.iter().map(|c| c.abs()).sum();
    for (re, im) in roots(&p) {
        // Evaluate p at the complex root by Horner.
        let (mut ar, mut ai) = (0.0f64, 0.0f64);
        for &c in &p {
            let (nr, ni) = (ar * re - ai * im + c, ar * im + ai * re);
            ar = nr;
            ai = ni;
        }
        assert!(ar.hypot(ai) < 1e-9 * scale, "p({re}+{im}i) = {ar}+{ai}i");
    }
}

#[test]
fn aberth_fallback_agrees_with_qr() {
    for p in [
        &[1.0, -6.0, 11.0, -6.0][..],
        &[1.0, 0.0, 1.0],
        &[1.0, 2.0, 3.0, 4.0, 5.0],
        &[3.0, 0.0, -2.0, 5.0, 1.0, -7.0, 2.0],
        &[1.0, 0.0, 0.0, 0.0, 0.0, -1.0],
    ] {
        same_roots(&roots_aberth(p), &roots(p), 1e-9, &format!("{p:?}"));
    }
}

#[test]
fn division_identity_u_equals_q_times_v_plus_r() {
    let u = [2.0, -1.0, 0.5, 3.0, -4.0, 1.0];
    let v = [1.5, 0.0, -2.0];
    let (q, r) = polydiv(&u, &v);
    let back = polyadd(&polymul(&q, &v), &r);
    // Compare after aligning (polyadd keeps leading zeros, so lengths agree here).
    close(&back, &u, 1e-12, "q*v + r");
}

#[test]
fn derivative_of_integral_is_identity() {
    let p = [4.0, -3.0, 2.0, 0.0, 7.0];
    for m in 1..4 {
        close(
            &polyder(&polyint(&p, m, &[]), m),
            &p,
            1e-14,
            &format!("m={m}"),
        );
    }
}

#[test]
fn horner_agrees_with_product_form() {
    let rs = [-2.0, 0.5, 3.0];
    let mut p = vec![2.0];
    for &r in &rs {
        p = polymul(&p, &[1.0, -r]);
    }
    for x in [-3.0, -0.1, 0.0, 1.7, 5.0] {
        let want = 2.0 * (x - rs[0]) * (x - rs[1]) * (x - rs[2]);
        assert!((polyval(&p, x) - want).abs() < 1e-12 * (1.0 + want.abs()));
    }
}

#[test]
#[should_panic(expected = "same length")]
fn polyfit_length_mismatch_panics() {
    let _ = polyfit(&[1.0, 2.0], &[1.0], 1);
}

#[test]
#[should_panic(expected = "length 1 or >m")]
fn polyint_too_few_constants_panics() {
    let _ = polyint(&[1.0, 2.0], 3, &[1.0, 2.0]);
}