interpolation.scala

abstract class InterpolationMethod(val s: Seq[(Double, Double)]) extends (Double => Double) {
  protected val x_i = s.map(_._1)
  protected val y_i = s.map(_._2)
}

class Lagrange(s: Seq[(Double, Double)]) extends InterpolationMethod(s) {
  val n = s.size

  private def poly(i: Int)(x: Double): Double = (for (j <- 0 until n) yield if (j != i) ((x - x_i(j)) / (x_i(i) - x_i(j))) else 1).product

  def apply(x: Double): Double = (for (i <- 0 until n) yield y_i(i) * poly(i)(x)).sum

  def this(f: Double => Double, f2: Int => Double, n: Int) = this(for (i <- 1 until n) yield {
    val x = f2(i)
    (x, f(x))
  })

  def this(f: Double => Double, n: Int) = this(f, _.toDouble, n)
}

class Newton(s: Seq[(Double, Double)]) extends InterpolationMethod(s) {
  val c_dividers: scala.collection.mutable.Map[(Int, Int), Double] = scala.collection.mutable.Map.empty

  val n = s.size

  private def poly(i: Int)(x: Double): Double = (for (j <- 0 until i) yield x - x_i(j)).product

  def apply(x: Double): Double = (for (i <- 0 until n) yield divided(0, i) * poly(i)(x)).sum

  private def divided(start: Int, end: Int): Double =
    c_dividers.get((start, end)).getOrElse {
      val tmp = if (start == end) y_i(start) else ((divided(start + 1, end) - divided(start, end - 1)) / (x_i(end) - x_i(start)))
      c_dividers.put((start, end), tmp)
      tmp
    }

  def this(f: Double => Double, f2: Int => Double, n: Int) = this(for (i <- 1 until n) yield {
    val x = f2(i)
    (x, f(x))
  })

  def this(f: Double => Double, n: Int) = this(f, _.toDouble, n)
}