Integration.scala

object Integration {
  def integrate(fn: Double => Double,
                interval: (Double, Double),
                precision: Double = 0.0001): Double = {
    def areaBetween(start: Double, end: Double): Double = (fn(end) + fn(start)) / 2  * (end - start)

    def doIntegrate(start: Double, end: Double, area: Double): Double = {
      val mid = (start + end)/2
      val left = areaBetween(start, mid)
      val right = areaBetween(mid, end)

      if(left + right - area < precision) left + right
      else doIntegrate(start, mid, left) + doIntegrate(mid, end, right)
    }

    doIntegrate(interval._1, interval._2, areaBetween(interval._1, interval._2))
  }


  def main(args: Array[String]) = {
    val unitCircle = {x: Double => math.sqrt(1 - x * x)}
    val pi = integrate(unitCircle, (0, 1)) * 4

    println(pi)
  }
}


object IntegrationNDimensional {

  // R^k
  trait Euclidean[T] {
    def split(start: T, end: T): List[(T, T)]
    def measure(start: T, end: T): Double
  }

  def euclidean[T](s: (T, T) => List[(T, T)])(m: (T, T) => Double) = new Euclidean[T] {
    override def split(start: T, end: T): List[(T, T)] = s(start, end)
    override def measure(start: T, end: T): Double = m(start, end)
  }

  implicit val realLine: Euclidean[Double] = euclidean[Double]{
    (start, end) => val mid = (start + end) / 2; List((start, mid), (mid, end))
  }{
    (start, end) => math.abs(end - start)
  }

  implicit def moreDimension[T, K](implicit tE: Euclidean[T], kE: Euclidean[K]): Euclidean[(T, K)] = euclidean[(T, K)] { (start, end) =>
    for {
      (lS, lE) <- tE.split(start._1, end._1)
      (rS, rE) <- kE.split(start._2, end._2)
    } yield ((lS, rS), (lE, rE))
  }{(start, end) => tE.measure(start._1, end._1) * kE.measure(start._2, end._2)}

  def integrate[D, R](fn: D => Double,
                      interval: (D, D),
                      delta: Double = 0.0001)(implicit dE: Euclidean[D]): Double = {
    def sumIn(start: D, end: D): Double =
      (fn(start) + fn(end))/2  * dE.measure(start, end)

    def doIntegrate(start: D, end: D, prev: Double): Double = {
      val splits = dE.split(start, end)
        .map{ case (s, e) => (s, e) -> sumIn(s, e) }
      val total = splits.map(_._2).sum

      if(math.abs(total - prev) < delta) total
      else splits.map {
        case ((s, e), sum) => doIntegrate(s, e, sum)
      }.sum
    }

    doIntegrate(interval._1, interval._2, sumIn(interval._1, interval._2))
  }


  def main(args: Array[String]) = {
    val unitCircle: Double => Double = {x: Double => math.sqrt(1 - x * x)}
    val unitSphere: ((Double, Double)) => Double = {case (x: Double, y: Double) => if(1 - x * x - y * y > 0) math.sqrt(1 - x * x - y * y) else 0}
    val pi = integrate(unitCircle, (0.0, 1.0), 0.00001) * 4

    val piFromSphere = integrate(unitSphere, ((0.0, 0.0), (1.0, 1.0)), 0.0000001) * 6

    println(pi)
    println(piFromSphere)
  }
}

How can we do this for improper Reimann integrals ? How about Lebesgue integrals?