use bone_solver::{
    ConstraintSystem, CurveRadius3, DofConfig, Ellipse3Handles, Line3Handle, NewtonConfig,
    Point3Handle, Residual, Triplet, analyze_dof, solve_newton,
};
use bone_types::{Parameter, ParameterIndex, ResidualIndex};

fn p(idx: u32) -> ParameterIndex {
    ParameterIndex::new(idx)
}

fn point3(x: u32, y: u32, z: u32) -> Point3Handle {
    Point3Handle {
        x: p(x),
        y: p(y),
        z: p(z),
    }
}

fn parameters(values: &[f64]) -> Vec<Parameter> {
    values.iter().copied().map(Parameter::new).collect()
}

fn analytic_jacobian(residual: &Residual, params: &[f64]) -> Vec<Triplet> {
    let mut sink = Vec::new();
    residual.jacobian(params, ResidualIndex::new(0), &mut sink);
    sink
}

fn numeric_derivative(residual: &Residual, params: &[f64], col: usize, row: usize) -> f64 {
    let rows = residual.rows();
    let h = 1e-6;
    let perturb = |sign: f64| -> Vec<f64> {
        params
            .iter()
            .enumerate()
            .map(|(i, value)| if i == col { value + sign * h } else { *value })
            .collect()
    };
    let plus = perturb(1.0);
    let minus = perturb(-1.0);
    let eval = |state: &[f64]| -> f64 {
        let mut out = vec![0.0; rows];
        residual.evaluate(state, &mut out);
        out[row]
    };
    (eval(&plus) - eval(&minus)) / (2.0 * h)
}

fn assert_jacobian_matches_finite_difference(residual: &Residual, params: &[f64]) {
    let rows = residual.rows();
    let cols = params.len();
    let triplets = analytic_jacobian(residual, params);
    (0..rows)
        .flat_map(|row| (0..cols).map(move |col| (row, col)))
        .for_each(|(row, col)| {
            let analytic: f64 = triplets
                .iter()
                .filter(|(r, c, _)| r.as_usize() == row && c.as_usize() == col)
                .map(|(_, _, v)| *v)
                .sum();
            let numeric = numeric_derivative(residual, params, col, row);
            let tolerance = 1e-4 * (1.0 + analytic.abs() + numeric.abs());
            assert!(
                (analytic - numeric).abs() < tolerance,
                "row {row} col {col}: analytic {analytic} vs numeric {numeric}",
            );
        });
}

#[test]
fn coincident3_jacobian_matches_finite_difference() {
    let residual = Residual::Coincident3(point3(0, 1, 2), point3(3, 4, 5));
    assert_eq!(residual.rows(), 3);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9]);
}

#[test]
fn distance3_jacobian_matches_finite_difference() {
    let residual = Residual::Distance3 {
        a: point3(0, 1, 2),
        b: point3(3, 4, 5),
        value_mm: 4.0,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9]);
}

#[test]
fn on_plane_jacobian_matches_finite_difference() {
    let residual = Residual::OnPlane {
        point: point3(0, 1, 2),
        normal: [0.3, -0.8, 0.5],
        offset_mm: 1.25,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(&residual, &[1.1, -2.2, 0.7]);
}

#[test]
fn coincident_point_ellipse3_jacobian_matches_finite_difference() {
    let residual = Residual::CoincidentPointEllipse3(Ellipse3Handles {
        point: point3(0, 1, 2),
        center: point3(3, 4, 5),
        major: point3(6, 7, 8),
        minor: point3(9, 10, 11),
    });
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[1.0, 2.0, 1.5, 0.0, 0.0, 0.0, 3.0, 0.0, 4.0, 0.0, 2.0, 0.0],
    );
}

#[test]
fn along_axis_jacobian_matches_finite_difference() {
    let residual = Residual::AlongAxis {
        line: Line3Handle {
            a: point3(0, 1, 2),
            b: point3(3, 4, 5),
        },
        basis: [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
    };
    assert_eq!(residual.rows(), 2);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9]);
}

#[test]
fn perpendicular3_jacobian_matches_finite_difference() {
    let residual = Residual::Perpendicular3(
        Line3Handle {
            a: point3(0, 1, 2),
            b: point3(3, 4, 5),
        },
        Line3Handle {
            a: point3(6, 7, 8),
            b: point3(9, 10, 11),
        },
    );
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[
            0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5, -1.5, 2.5,
        ],
    );
}

#[test]
fn parallel3_jacobian_matches_finite_difference() {
    let residual = Residual::Parallel3(
        Line3Handle {
            a: point3(0, 1, 2),
            b: point3(3, 4, 5),
        },
        Line3Handle {
            a: point3(6, 7, 8),
            b: point3(9, 10, 11),
        },
    );
    assert_eq!(residual.rows(), 3);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[
            0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5, -1.5, 2.5,
        ],
    );
}

#[test]
fn point_on_line3_jacobian_matches_finite_difference() {
    let residual = Residual::PointOnLine3 {
        point: point3(0, 1, 2),
        line: Line3Handle {
            a: point3(3, 4, 5),
            b: point3(6, 7, 8),
        },
    };
    assert_eq!(residual.rows(), 3);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0],
    );
}

#[test]
fn point_on_line3_pulls_free_point_onto_the_line() {
    let point = point3(0, 1, 2);
    let line = Line3Handle {
        a: point3(3, 4, 5),
        b: point3(6, 7, 8),
    };
    let pins = [
        (line.a.x, 0.0),
        (line.a.y, 0.0),
        (line.a.z, 0.0),
        (line.b.x, 1.0),
        (line.b.y, 0.0),
        (line.b.z, 0.0),
    ]
    .into_iter()
    .map(|(param, target)| Residual::Pin { param, target });
    let residuals: Vec<Residual> = pins
        .chain(std::iter::once(Residual::PointOnLine3 { point, line }))
        .collect();
    let system = ConstraintSystem::new(
        parameters(&[2.0, 3.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]),
        residuals,
    );
    let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("point-on-line over a pinned x-axis line is linear, must converge");
    };
    assert!(out[1].value().abs() < 1e-6, "point y {}", out[1].value());
    assert!(out[2].value().abs() < 1e-6, "point z {}", out[2].value());
    assert!(
        (out[0].value() - 2.0).abs() < 1e-6,
        "free slide stays put {}",
        out[0].value()
    );
}

#[test]
fn point_on_line3_leaves_one_sliding_dof() {
    let point = point3(0, 1, 2);
    let line = Line3Handle {
        a: point3(3, 4, 5),
        b: point3(6, 7, 8),
    };
    let pins = [
        (line.a.x, 0.0),
        (line.a.y, 0.0),
        (line.a.z, 0.0),
        (line.b.x, 1.0),
        (line.b.y, 0.0),
        (line.b.z, 0.0),
    ]
    .into_iter()
    .map(|(param, target)| Residual::Pin { param, target });
    let residuals: Vec<Residual> = pins
        .chain(std::iter::once(Residual::PointOnLine3 { point, line }))
        .collect();
    let system = ConstraintSystem::new(
        parameters(&[2.0, 0.5, -0.5, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]),
        residuals,
    );
    let report = analyze_dof(&system, DofConfig::DEFAULT);
    assert_eq!(
        report.dof().value(),
        1,
        "the point may slide along the line"
    );
    assert!(report.over_constrained().is_empty());
}

#[test]
fn coincident3_solves_second_point_to_pinned_first() {
    let first = point3(0, 1, 2);
    let second = point3(3, 4, 5);
    let pins = [(first.x, 1.0), (first.y, -2.0), (first.z, 3.5)]
        .into_iter()
        .map(|(param, target)| Residual::Pin { param, target });
    let residuals: Vec<Residual> = pins
        .chain(std::iter::once(Residual::Coincident3(second, first)))
        .collect();
    let system = ConstraintSystem::new(parameters(&[0.0, 0.0, 0.0, 9.0, 9.0, 9.0]), residuals);
    let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("coincident over pinned point is linear, must converge");
    };
    assert!(
        (out[3].value() - 1.0).abs() < 1e-6,
        "b.x {}",
        out[3].value()
    );
    assert!(
        (out[4].value() - -2.0).abs() < 1e-6,
        "b.y {}",
        out[4].value()
    );
    assert!(
        (out[5].value() - 3.5).abs() < 1e-6,
        "b.z {}",
        out[5].value()
    );
}

#[test]
fn along_axis_snaps_free_line_parallel_to_x_axis() {
    let line = Line3Handle {
        a: point3(0, 1, 2),
        b: point3(3, 4, 5),
    };
    let pins = [
        (line.a.x, 0.0),
        (line.a.y, 0.0),
        (line.a.z, 0.0),
        (line.b.x, 5.0),
    ]
    .into_iter()
    .map(|(param, target)| Residual::Pin { param, target });
    let residuals: Vec<Residual> = pins
        .chain(std::iter::once(Residual::AlongAxis {
            line,
            basis: [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
        }))
        .collect();
    let system = ConstraintSystem::new(parameters(&[0.0, 0.0, 0.0, 5.0, 2.7, -1.3]), residuals);
    let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("along-axis snap to x must converge");
    };
    let dir = (
        out[3].value() - out[0].value(),
        out[4].value() - out[1].value(),
        out[5].value() - out[2].value(),
    );
    assert!(dir.1.abs() < 1e-6, "direction y component {}", dir.1);
    assert!(dir.2.abs() < 1e-6, "direction z component {}", dir.2);
    assert!(
        dir.0.abs() > 1e-3,
        "direction x component must stay nonzero {}",
        dir.0
    );
}

#[test]
fn free_3d_point_has_three_dof_and_pins_remove_all() {
    let point = point3(0, 1, 2);
    let free = ConstraintSystem::new(
        parameters(&[0.5, -1.0, 2.0]),
        vec![Residual::OnPlane {
            point,
            normal: [0.0, 0.0, 1.0],
            offset_mm: 2.0,
        }],
    );
    let free_report = analyze_dof(&free, DofConfig::DEFAULT);
    assert_eq!(free_report.dof().value(), 2);

    let unconstrained = ConstraintSystem::new(parameters(&[0.5, -1.0, 2.0]), vec![]);
    let unconstrained_report = analyze_dof(&unconstrained, DofConfig::DEFAULT);
    assert_eq!(unconstrained_report.dof().value(), 3);
    assert_eq!(unconstrained_report.under_constrained().len(), 3);

    let pinned = ConstraintSystem::new(
        parameters(&[0.5, -1.0, 2.0]),
        vec![
            Residual::Pin {
                param: point.x,
                target: 0.5,
            },
            Residual::Pin {
                param: point.y,
                target: -1.0,
            },
            Residual::Pin {
                param: point.z,
                target: 2.0,
            },
        ],
    );
    let pinned_report = analyze_dof(&pinned, DofConfig::DEFAULT);
    assert_eq!(pinned_report.dof().value(), 0);
    assert!(pinned_report.under_constrained().is_empty());
    assert!(pinned_report.over_constrained().is_empty());
}

fn line3(a: (u32, u32, u32), b: (u32, u32, u32)) -> Line3Handle {
    Line3Handle {
        a: point3(a.0, a.1, a.2),
        b: point3(b.0, b.1, b.2),
    }
}

#[test]
fn angular_between_lines3_jacobian_matches_finite_difference() {
    let residual = Residual::AngularBetweenLines3 {
        a: line3((0, 1, 2), (3, 4, 5)),
        b: line3((6, 7, 8), (9, 10, 11)),
        angle_rad: core::f64::consts::FRAC_PI_3,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[
            0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5, -1.5, 2.5,
        ],
    );
}

#[test]
fn equal_distance3_jacobian_matches_finite_difference() {
    let residual = Residual::EqualDistance3 {
        center: point3(0, 1, 2),
        a: point3(3, 4, 5),
        b: point3(6, 7, 8),
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0],
    );
}

#[test]
fn orthogonal_spokes3_jacobian_matches_finite_difference() {
    let residual = Residual::OrthogonalSpokes3 {
        center: point3(0, 1, 2),
        a: point3(3, 4, 5),
        b: point3(6, 7, 8),
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0],
    );
}

#[test]
fn equal_distance3_drives_two_spokes_to_one_radius() {
    let center = point3(0, 1, 2);
    let a = point3(3, 4, 5);
    let b = point3(6, 7, 8);
    let system = ConstraintSystem::new(
        parameters(&[0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 9.0, 0.0]),
        vec![
            Residual::Pin {
                param: center.x,
                target: 0.0,
            },
            Residual::Pin {
                param: center.y,
                target: 0.0,
            },
            Residual::Pin {
                param: center.z,
                target: 0.0,
            },
            Residual::Pin {
                param: a.x,
                target: 4.0,
            },
            Residual::Pin {
                param: a.y,
                target: 0.0,
            },
            Residual::Pin {
                param: a.z,
                target: 0.0,
            },
            Residual::Pin {
                param: b.x,
                target: 0.0,
            },
            Residual::Pin {
                param: b.z,
                target: 0.0,
            },
            Residual::EqualDistance3 { center, a, b },
        ],
    );
    let Ok(solved) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("equal-distance over a pinned spoke is well-posed");
    };
    let bx = solved[b.x.as_usize()].value();
    let by = solved[b.y.as_usize()].value();
    let bz = solved[b.z.as_usize()].value();
    let radius_b = (bx * bx + by * by + bz * bz).sqrt();
    assert!(
        (radius_b - 4.0).abs() < 1e-6,
        "second spoke matches the first radius: {radius_b}"
    );
    assert!(by.abs() > 1e-3, "the free axis stays non-degenerate: {by}");
}

fn explicit3(center: (u32, u32, u32), radius: u32) -> CurveRadius3 {
    CurveRadius3::Explicit {
        center: point3(center.0, center.1, center.2),
        radius: p(radius),
    }
}

fn spoke3(center: (u32, u32, u32), spoke: (u32, u32, u32)) -> CurveRadius3 {
    CurveRadius3::FromSpoke {
        center: point3(center.0, center.1, center.2),
        spoke: point3(spoke.0, spoke.1, spoke.2),
    }
}

#[test]
fn radius_curve3_explicit_jacobian_matches_finite_difference() {
    let residual = Residual::RadiusCurve3 {
        curve: explicit3((0, 1, 2), 3),
        value_mm: 2.5,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4]);
}

#[test]
fn radius_curve3_spoke_jacobian_matches_finite_difference() {
    let residual = Residual::RadiusCurve3 {
        curve: spoke3((0, 1, 2), (3, 4, 5)),
        value_mm: 2.5,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9]);
}

#[test]
fn coincident_point_curve3_explicit_jacobian_matches_finite_difference() {
    let residual = Residual::CoincidentPointCurve3 {
        point: point3(0, 1, 2),
        curve: explicit3((3, 4, 5), 6),
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 2.5]);
}

#[test]
fn coincident_point_curve3_spoke_jacobian_matches_finite_difference() {
    let residual = Residual::CoincidentPointCurve3 {
        point: point3(0, 1, 2),
        curve: spoke3((3, 4, 5), (6, 7, 8)),
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0],
    );
}

#[test]
fn equal_radius3_mixed_jacobian_matches_finite_difference() {
    let residual = Residual::EqualRadius3(explicit3((0, 1, 2), 3), spoke3((4, 5, 6), (7, 8, 9)));
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5],
    );
}

#[test]
fn equal_length3_jacobian_matches_finite_difference() {
    let residual =
        Residual::EqualLength3(line3((0, 1, 2), (3, 4, 5)), line3((6, 7, 8), (9, 10, 11)));
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[
            0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5, -1.5, 2.5,
        ],
    );
}

#[test]
fn tangent_line_curve3_explicit_jacobian_matches_finite_difference() {
    let residual = Residual::TangentLineCurve3 {
        line: line3((0, 1, 2), (3, 4, 5)),
        curve: explicit3((6, 7, 8), 9),
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 1.5],
    );
}

#[test]
fn tangent_line_curve3_spoke_jacobian_matches_finite_difference() {
    let residual = Residual::TangentLineCurve3 {
        line: line3((0, 1, 2), (3, 4, 5)),
        curve: spoke3((6, 7, 8), (9, 10, 11)),
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[
            0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5, -1.5, 2.5,
        ],
    );
}

#[test]
fn tangent_curve_curve3_jacobian_matches_finite_difference() {
    let residual = Residual::TangentCurveCurve3 {
        a: explicit3((0, 1, 2), 3),
        b: spoke3((4, 5, 6), (7, 8, 9)),
        internal: false,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 1.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5],
    );
}

#[test]
fn tangent_curve_curve3_internal_jacobian_matches_finite_difference() {
    let residual = Residual::TangentCurveCurve3 {
        a: explicit3((0, 1, 2), 3),
        b: spoke3((4, 5, 6), (7, 8, 9)),
        internal: true,
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[0.3, -1.7, 2.1, 1.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5],
    );
}

#[test]
fn in_plane3_jacobian_matches_finite_difference() {
    let residual = Residual::InPlane3 {
        point: point3(0, 1, 2),
        origin: point3(3, 4, 5),
        normal: [0.0, 0.0, 1.0],
    };
    assert_eq!(residual.rows(), 1);
    assert_jacobian_matches_finite_difference(&residual, &[0.3, -1.7, 2.1, 4.4, -0.6, 5.9]);
}

#[test]
fn symmetric3_jacobian_matches_finite_difference() {
    let residual = Residual::Symmetric3 {
        a: point3(0, 1, 2),
        b: point3(3, 4, 5),
        axis: Line3Handle {
            a: point3(6, 7, 8),
            b: point3(9, 10, 11),
        },
    };
    assert_eq!(residual.rows(), 4);
    assert_jacobian_matches_finite_difference(
        &residual,
        &[
            0.3, -1.7, 2.1, 4.4, -0.6, 5.9, 1.0, 2.0, -3.0, 0.5, -1.5, 2.5,
        ],
    );
}

#[test]
fn tangent_curve_curve3_drives_centers_to_radius_sum() {
    let a = explicit3((0, 1, 2), 3);
    let b = explicit3((4, 5, 6), 7);
    let CurveRadius3::Explicit {
        center: ca,
        radius: ra,
    } = a
    else {
        unreachable!()
    };
    let CurveRadius3::Explicit {
        center: cb,
        radius: rb,
    } = b
    else {
        unreachable!()
    };
    let system = ConstraintSystem::new(
        parameters(&[0.0, 0.0, 0.0, 2.0, 9.0, 0.0, 0.0, 3.0]),
        vec![
            Residual::Pin {
                param: ca.x,
                target: 0.0,
            },
            Residual::Pin {
                param: ca.y,
                target: 0.0,
            },
            Residual::Pin {
                param: ca.z,
                target: 0.0,
            },
            Residual::Pin {
                param: ra,
                target: 2.0,
            },
            Residual::Pin {
                param: cb.y,
                target: 0.0,
            },
            Residual::Pin {
                param: cb.z,
                target: 0.0,
            },
            Residual::Pin {
                param: rb,
                target: 3.0,
            },
            Residual::TangentCurveCurve3 {
                a,
                b,
                internal: false,
            },
        ],
    );
    let Ok(solved) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("external tangency over pinned radii is well-posed");
    };
    let dx = solved[cb.x.as_usize()].value();
    assert!(
        (dx.abs() - 5.0).abs() < 1e-6,
        "centers separate by the radius sum: {dx}"
    );
}

#[test]
fn internal_tangent_curve_curve3_drives_centers_to_radius_difference() {
    let a = explicit3((0, 1, 2), 3);
    let b = explicit3((4, 5, 6), 7);
    let CurveRadius3::Explicit {
        center: ca,
        radius: ra,
    } = a
    else {
        unreachable!()
    };
    let CurveRadius3::Explicit {
        center: cb,
        radius: rb,
    } = b
    else {
        unreachable!()
    };
    let system = ConstraintSystem::new(
        parameters(&[0.0, 0.0, 0.0, 5.0, 4.0, 0.0, 0.0, 2.0]),
        vec![
            Residual::Pin {
                param: ca.x,
                target: 0.0,
            },
            Residual::Pin {
                param: ca.y,
                target: 0.0,
            },
            Residual::Pin {
                param: ca.z,
                target: 0.0,
            },
            Residual::Pin {
                param: ra,
                target: 5.0,
            },
            Residual::Pin {
                param: cb.y,
                target: 0.0,
            },
            Residual::Pin {
                param: cb.z,
                target: 0.0,
            },
            Residual::Pin {
                param: rb,
                target: 2.0,
            },
            Residual::TangentCurveCurve3 {
                a,
                b,
                internal: true,
            },
        ],
    );
    let Ok(solved) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("internal tangency over pinned radii is well-posed");
    };
    let dx = solved[cb.x.as_usize()].value();
    assert!(
        (dx.abs() - 3.0).abs() < 1e-6,
        "centers separate by the radius difference: {dx}"
    );
}
