///|
/// Solve (W + lambda D2' D2) z = W y using pentadiagonal Cholesky.
/// The penalty is in sample-index coordinates; use a common uniform axis for
/// physically comparable fits. Positive weights keep the system definite.
pub fn whittaker(
  values : Array[Double],
  weights : Array[Double],
  lambda : Double,
) -> Array[Double] raise SpectralError {
  validate_values(values, 3)
  validate_iteration_options(iteration_options(lambda~))
  let n = values.length()
  if weights.length() != n {
    raise Invalid("weights length differs")
  }
  let diagonal = Array::make(n, 0.0)
  let first = Array::make(n - 1, 0.0)
  let second = Array::make(n - 2, 0.0)
  let rhs = Array::make(n, 0.0)
  for i in 0.. 1.0e8 {
      raise Invalid("weights must lie in [1e-8,1e8]")
    }
    diagonal[i] = weights[i]
    rhs[i] = checked_result(values[i] * weights[i], "weighted sample")
  }
  // Assemble each [1,-2,1] outer product, including short-array boundaries.
  for i in 0..<(n - 2) {
    diagonal[i] = diagonal[i] + lambda
    diagonal[i + 1] = diagonal[i + 1] + 4.0 * lambda
    diagonal[i + 2] = diagonal[i + 2] + lambda
    first[i] = first[i] - 2.0 * lambda
    first[i + 1] = first[i + 1] - 2.0 * lambda
    second[i] = second[i] + lambda
  }
  let pivot = Array::make(n, 0.0)
  let lower_one = Array::make(n, 0.0)
  let lower_two = Array::make(n, 0.0)
  for i in 0..= 2 {
      lower_two[i] = second[i - 2] / pivot[i - 2]
    }
    if i >= 1 {
      lower_one[i] = (first[i - 1] - lower_two[i] * lower_one[i - 1]) /
        pivot[i - 1]
    }
    let d = diagonal[i] -
      lower_one[i] * lower_one[i] -
      lower_two[i] * lower_two[i]
    if d <= 0.0 || d.is_nan() || d.is_inf() {
      raise Numerical(
        "Whittaker factorization lost positive definiteness; reduce lambda or weight contrast",
      )
    }
    pivot[i] = d.sqrt()
  }
  let forward = Array::make(n, 0.0)
  for i in 0..= 1 {
      value = value - lower_one[i] * forward[i - 1]
    }
    if i >= 2 {
      value = value - lower_two[i] * forward[i - 2]
    }
    forward[i] = checked_result(value / pivot[i], "forward substitution")
  }
  let result = Array::make(n, 0.0)
  for k in 0.. Double raise SpectralError {
  validate_values(values, 3)
  validate_values(baseline, 3)
  validate_iteration_options(iteration_options(lambda~))
  if baseline.length() != values.length() || weights.length() != values.length() {
    raise Invalid("objective lengths differ")
  }
  let mut objective = 0.0
  for i in 0..= 2 {
      let curvature = baseline[i] - 2.0 * baseline[i - 1] + baseline[i - 2]
      objective = objective + lambda * curvature * curvature
    }
  }
  checked_result(objective, "objective")
}