///|
fn same_parts(x : Array[Double], y : Array[Double]) -> Unit raise CodaError {
  if x.length() != y.length() {
    raise Shape("composition dimension mismatch")
  }
}

///|
/// Aitchison addition, implemented in log space.
pub fn perturb(
  x : Array[Double],
  y : Array[Double],
) -> Array[Double] raise CodaError {
  same_parts(x, y)
  let a = clr(x)
  let b = clr(y)
  clr_inverse(Array::makei(a.length(), i => checked(a[i] + b[i])))
}

///|
pub fn perturb_inverse(
  x : Array[Double],
  y : Array[Double],
) -> Array[Double] raise CodaError {
  same_parts(x, y)
  let a = clr(x)
  let b = clr(y)
  clr_inverse(Array::makei(a.length(), i => checked(a[i] - b[i])))
}

///|
pub fn power(
  x : Array[Double],
  exponent : Double,
) -> Array[Double] raise CodaError {
  if !finite(exponent) {
    raise Invalid("non-finite exponent")
  }
  let a = clr(x)
  clr_inverse(a.map(v => checked(v * exponent)))
}

///|
pub fn inner(x : Array[Double], y : Array[Double]) -> Double raise CodaError {
  same_parts(x, y)
  let a = clr(x)
  let b = clr(y)
  let mut s = 0.0
  for i in 0.. Double raise CodaError {
  same_parts(x, y)
  let a = clr(x)
  let b = clr(y)
  let mut s = 0.0
  for i in 0.. Array[Double] raise CodaError {
  same_parts(x, y)
  if !finite(t) || t.abs() > 1000.0 {
    raise Invalid("interpolation parameter outside -1000..1000")
  }
  let a = clr(x)
  let b = clr(y)
  clr_inverse(Array::makei(a.length(), i => a[i] + t * (b[i] - a[i])))
}