///|
pub fn trapezoid_integral_uniform(
  values : ArrayView[Double],
  spacing : Double,
) -> Double {
  if values.length() < 2 {
    0.0
  } else {
    let mut total = 0.5 * (values[0] + values[values.length() - 1])
    for i in 1..<(values.length() - 1) {
      total = total + values[i]
    }
    total * spacing
  }
}

///|
pub fn simpson_integral(values : ArrayView[Double], spacing : Double) -> Double {
  if values.length() < 2 {
    0.0
  } else if values.length() == 2 {
    trapezoid_integral_uniform(values, spacing)
  } else {
    let mut total = values[0] + values[values.length() - 1]
    for i in 1..<(values.length() - 1) {
      total = total + values[i] * (if i % 2 == 0 { 2.0 } else { 4.0 })
    }
    total * spacing / 3.0
  }
}

///|
pub fn finite_difference_first(
  values : ArrayView[Double],
  spacing : Double,
) -> Array[Double] {
  finite_difference_series(values, spacing)
}

///|
pub fn finite_difference_second(
  values : ArrayView[Double],
  spacing : Double,
) -> Array[Double] {
  let first = finite_difference_series(values, spacing)
  finite_difference_series(first, spacing)
}

///|
pub fn central_interpolate(
  first : Double,
  second : Double,
  fraction : Double,
) -> Double {
  lerp(first, second, fraction)
}

///|
pub fn polynomial_basis(x : Double, order : Int) -> Array[Double] {
  Array::makei(order.max(0) + 1, fn(i) { @math.pow(x, i.to_double()) })
}

///|
pub fn polynomial_evaluate(
  coefficients : ArrayView[Double],
  x : Double,
) -> Double {
  let mut result = 0.0
  for coefficient in coefficients.rev() {
    result = result * x + coefficient
  }
  result
}

///|
pub fn polynomial_derivative(coefficients : ArrayView[Double]) -> Array[Double] {
  if coefficients.length() <= 1 {
    [0.0]
  } else {
    Array::makei(coefficients.length() - 1, fn(i) {
      coefficients[i + 1] * (i + 1).to_double()
    })
  }
}

///|
pub fn polynomial_fit_constant(values : ArrayView[Double]) -> Array[Double] {
  [mean(values)]
}

///|
pub fn residual_sum_of_squares(
  actual : ArrayView[Double],
  predicted : ArrayView[Double],
) -> Double {
  let mut total = 0.0
  let count = actual.length().min(predicted.length())
  for i in 0.. Double {
  let count = actual.length().min(predicted.length())
  if count == 0 {
    0.0
  } else {
    let mut total = 0.0
    for i in 0.. Double {
  safe_ratio((actual - predicted).abs(), actual.abs().max(1.0e-30), 0.0)
}

///|
pub fn convergence_order(
  errors : ArrayView[Double],
  step_ratio : Double,
) -> Double {
  if errors.length() < 2 || step_ratio <= 0.0 || step_ratio == 1.0 {
    0.0
  } else {
    approximate_log(errors[0] / errors[errors.length() - 1]) /
    approximate_log(step_ratio)
  }
}

///|
pub fn geometric_schedule(
  start : Double,
  factor : Double,
  count : Int,
) -> Array[Double] {
  Array::makei(count.max(0), fn(i) { start * @math.pow(factor, i.to_double()) })
}

///|
pub fn linear_schedule(
  start : Double,
  stop : Double,
  count : Int,
) -> Array[Double] {
  linspace(start, stop, count)
}

///|
pub fn stable_time_step(
  cell_width : Double,
  speed : Double,
  cfl : Double,
) -> Double {
  safe_ratio(cell_width * cfl, speed.abs(), cell_width * cfl)
}

///|
pub fn clamp_time_step(
  value : Double,
  minimum : Double,
  maximum : Double,
) -> Double {
  clamp(value, minimum, maximum)
}

///|
pub fn numerical_epsilon(scale : Double) -> Double {
  1.0e-12 * scale.abs().max(1.0)
}