///|
/// Polynomial coefficients refer to t=(x-center)/scale; arrays stay private.
pub struct Polynomial {
  coefficients : Array[Double]
  center : Double
  scale : Double
}

///|
pub fn Polynomial::coefficients(self : Polynomial) -> Array[Double] {
  self.coefficients.copy()
}

///|
pub fn Polynomial::center(self : Polynomial) -> Double {
  self.center
}

///|
pub fn Polynomial::scale(self : Polynomial) -> Double {
  self.scale
}

///|
pub fn Polynomial::degree(self : Polynomial) -> Int32 {
  self.coefficients.length() - 1
}

///|
/// Evaluate a physical-coordinate derivative (order 0 returns fitted value).
pub fn Polynomial::evaluate(
  self : Polynomial,
  x : Double,
  derivative~ : Int32 = 0,
) -> Double raise SpectralError {
  finite(x, "evaluation coordinate")
  if derivative < 0 || derivative > 8 {
    raise Invalid("derivative order must lie in [0,8]")
  }
  if derivative > self.degree() {
    return 0.0
  }
  let t = (x - self.center) / self.scale
  let mut value = 0.0
  for offset in 0..<(self.coefficients.length() - derivative) {
    let i = self.degree() - offset
    let mut factor = 1.0
    for j in 0.. Polynomial raise SpectralError {
  validate_values(values, 2)
  if axis.length() != values.length() {
    raise Invalid("polynomial lengths differ")
  }
  if degree < 0 || degree > 8 || degree >= values.length() {
    raise Invalid(
      "polynomial degree must be 0..8 and less than selected samples",
    )
  }
  for i in 0.. 0 && axis[i] <= axis[i - 1] {
      raise Invalid("polynomial axis must increase")
    }
  }
  let span = checked_result(axis[axis.length() - 1] - axis[0], "axis span")
  let scale = span / 2.0
  if scale <= 0.0 {
    raise Numerical("polynomial axis scale vanished")
  }
  let center = axis[0] + scale
  let n = values.length()
  let columns = degree + 1
  let basis : Array[Array[Double]] = []
  let upper = Array::makei(columns, fn(_) { Array::make(columns, 0.0) })
  let projected = Array::make(columns, 0.0)
  let normalized = axis.map(fn(x) { (x - center) / scale })
  let powers = Array::make(n, 1.0)
  for j in 0..