///|
/// Compensated sum reduces cancellation without changing sample order.
pub fn stable_sum(values : Array[Double]) -> Double raise SpectralError {
  validate_values(values, 1)
  let mut sum = 0.0
  let mut correction = 0.0
  for value in values {
    let adjusted = value - correction
    let next = sum + adjusted
    correction = next - sum - adjusted
    sum = next
  }
  checked_result(sum, "sum")
}

///|
fn checked_result(x : Double, label : String) -> Double raise SpectralError {
  if x.is_nan() || x.is_inf() {
    raise Numerical(label + " overflow or nonfinite result")
  }
  x
}

///|
pub fn mean(values : Array[Double]) -> Double raise SpectralError {
  stable_sum(values) / values.length().to_double()
}

///|
/// Sample standard deviation, using a second pass and scaled differences.
pub fn sample_deviation(values : Array[Double]) -> Double raise SpectralError {
  validate_values(values, 2)
  let center = mean(values)
  let mut scale = 0.0
  let mut squares = 1.0
  for value in values {
    let d = (value - center).abs()
    if d > 0.0 {
      if scale < d {
        let ratio = scale / d
        squares = 1.0 + squares * ratio * ratio
        scale = d
      } else {
        let ratio = d / scale
        squares = squares + ratio * ratio
      }
    }
  }
  checked_result(
    scale * (squares / (values.length() - 1).to_double()).sqrt(),
    "deviation",
  )
}

///|
pub fn root_mean_square(values : Array[Double]) -> Double raise SpectralError {
  validate_values(values, 1)
  let mut scale = 0.0
  for v in values {
    if v.abs() > scale {
      scale = v.abs()
    }
  }
  if scale == 0.0 {
    return 0.0
  }
  let mut squares = 0.0
  for v in values {
    let r = v / scale
    squares = squares + r * r
  }
  checked_result(scale * (squares / values.length().to_double()).sqrt(), "rms")
}

///|
/// Linear interpolation quantile at p*(n-1), on an owned sorted copy.
pub fn quantile(
  values : Array[Double],
  probability : Double,
) -> Double raise SpectralError {
  validate_values(values, 1)
  finite(probability, "probability")
  if probability < 0.0 || probability > 1.0 {
    raise Invalid("probability outside [0,1]")
  }
  let sorted = values.copy()
  sorted.sort()
  let position = probability * (sorted.length() - 1).to_double()
  let lower = position.to_int()
  if lower == sorted.length() - 1 {
    return sorted[lower]
  }
  let fraction = position - lower.to_double()
  sorted[lower] * (1.0 - fraction) + sorted[lower + 1] * fraction
}

///|
/// Median absolute deviation; no Gaussian scale factor is applied.
pub fn median_absolute_deviation(
  values : Array[Double],
) -> Double raise SpectralError {
  let median = quantile(values, 0.5)
  quantile(values.map(fn(v) { (v - median).abs() }), 0.5)
}

///|
fn relative_difference(
  next : Array[Double],
  previous : Array[Double],
) -> Double raise SpectralError {
  let delta = Array::makei(next.length(), fn(i) { next[i] - previous[i] })
  let denominator = root_mean_square(previous)
  if denominator <= 1.0e-100 {
    return root_mean_square(delta)
  }
  checked_result(root_mean_square(delta) / denominator, "relative change")
}

///|
/// Relative RMS difference, with absolute RMS fallback for a zero reference.
pub fn relative_rms(
  next : Array[Double],
  previous : Array[Double],
) -> Double raise SpectralError {
  validate_values(next, 1)
  validate_values(previous, 1)
  if next.length() != previous.length() {
    raise Invalid("relative RMS lengths differ")
  }
  relative_difference(next, previous)
}