Function normal_equations 

fn normal_equations(
    jacobian: &[Vec<f64>],
    residuals: &[f64],
    p: usize,
) -> (Vec<Vec<f64>>, Vec<f64>)

JᵀJ and Jᵀr, both p-sized — the matrices this method inverts are the parameter count square, never the bond count. JᵀJ is symmetric, so one triangle is accumulated and mirrored.