{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# MATH50003 (2022–23)\n",
    "# Lab 2: Interval arithmetic"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This lab explores the usage of rounding modes for floating point arithmetic and how they\n",
    "can be used to compute _rigorous_ bounds on mathematical constants such as ℯ.\n",
    "The key idea is _interval arithmetic_.\n",
    "That is recall the set operations\n",
    "$$\n",
    "A + B = \\{x + y : x ∈ A, y ∈ B\\}, AB = \\{xy : x ∈ A, y ∈ B\\}, A/B = \\{x/y : x ∈ A, y ∈ B\\}\n",
    "$$\n",
    "\n",
    "We will use floating point arithmetic to construct approximate set operations ⊕, ⊗ so that\n",
    "$$\n",
    "  A + B ⊆ A ⊕ B, AB ⊆ A ⊗ B, A/B ⊆ A ⊘ B\n",
    "$$\n",
    "thereby a complicated algorithm can be run on sets and the true result is guaranteed to be\n",
    "a subset of the output. E.g. we can do $ℯ = {\\rm exp}(1) ∈ {\\rm exp}([1,1]) ⊆ {\\rm exp}^{\\rm FP}([1,1])$\n",
    "where ${\\rm exp}^{\\rm FP}$ is implemented using $⊕$ and $⊗$.\n",
    "\n",
    "This will be consist of the following:\n",
    "1. The finite Taylor series $\\exp x ≈ ∑_{k=0}^n x^k/k!$ where each operation is now\n",
    "   an interval operation\n",
    "2. A bound on $∑_{k=n+1}^∞ x^k/k!$ that we capture in the returned result\n",
    "\n",
    "\n",
    "In what follows, the starred (⋆) problems are meant to be done with pen-and-paper.\n",
    "We need the following packages:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SetRounding, Test"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "-----\n",
    "**Problem 5.1⋆** For intervals $A = [a,b]$ and $B = [c,d]$ such that $0 ∉ A,B$\n",
    " and integer $n ≠ 0$,\n",
    "deduce formulas for the minimum and maximum of $A/n$, $A+B$ and $AB$.\n",
    "-----\n",
    "We want to implement floating point variants such that, for $S = [a,b] + [c,d]$\n",
    " $P = [a,b] * [c,d]$, and $D = [a,b]/n$ for an integer $n$,\n",
    "$$\n",
    "\\begin{align*}\n",
    "[a,b] ⊕ [c,d] &:= [{\\rm fl}^{\\rm down}(\\min S), {\\rm fl}^{\\rm up}(\\max S)] \\\\\n",
    "[a,b] ⊗ [c,d] &:= [{\\rm fl}^{\\rm down}(\\min P), {\\rm fl}^{\\rm up}(\\max P)] \\\\\n",
    "[a,b] ⊘ n &:= [{\\rm fl}^{\\rm down}(\\min D), {\\rm fl}^{\\rm up}(\\max D)]\n",
    "\\end{align*}\n",
    "$$\n",
    "This guarantees $S ⊆ [a,b] ⊕ [c,d]$, $P ⊆ [a,b] ⊗ [c,d]$, and\n",
    "$D ⊆ [a,b] ⊘ n$.\n",
    "In other words, if $x ∈ [a,b]$ and\n",
    "$y ∈ [c,d]$ then $x +y ∈ [a,b] ⊕ [c,d]$, and we thereby have  bounds on $x + y$.\n",
    "\n",
    "-------\n",
    "**Problem 5.2**  Use the formulae from Problem 5.1 to complete (by replacing the `# TODO: …` comments with code)\n",
    " the following implementation of an\n",
    "`Interval`\n",
    "so that `+`, `-`, and `/` implement $⊕$, $⊖$, and $⊘$ as defined above."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "\u001b[33m\u001b[1mTest Broken\u001b[22m\u001b[39m\n  Expression: Interval(-1.2, -1.1) * Interval(-3.1, -2.1) ≡ Interval(2.31, 3.72)"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "# Interval(a,b) represents the closed interval [a,b]\n",
    "# We use templating so that T can be e.g. a `BigFloat`\n",
    "struct Interval{T}\n",
    "    a::T\n",
    "    b::T\n",
    "end\n",
    "\n",
    "# We will overload *, +, -, / to use interval arithmetic.\n",
    "# We also need to support `one` and `in`\n",
    "import Base: *, +, -, /, one, in\n",
    "\n",
    "# create an interval corresponding to [1,1]\n",
    "one(x::Interval) = Interval(one(x.a), one(x.b))\n",
    "\n",
    "# Support x in Interval(a,b)\n",
    "in(x, y::Interval) = y.a ≤ x ≤ y.b\n",
    "\n",
    "# Following should implement ⊕\n",
    "function +(x::Interval, y::Interval)\n",
    "    T = promote_type(typeof(x.a), typeof(x.b))\n",
    "    a = setrounding(T, RoundDown) do\n",
    "        # TODO: lower bound\n",
    "\n",
    "    end\n",
    "    b = setrounding(T, RoundUp) do\n",
    "        # TODO: upper bound\n",
    "\n",
    "    end\n",
    "    Interval(a, b)\n",
    "end\n",
    "\n",
    "# following example was the non-associative example but now we have bounds\n",
    "@test_broken Interval(1.1,1.1) + Interval(1.2,1.2) + Interval(1.3,1.3) ≡ Interval(3.5999999999999996, 3.6000000000000005)\n",
    "\n",
    "# Following should implement ⊘\n",
    "function /(x::Interval, n::Integer)\n",
    "    T = typeof(x.a)\n",
    "    if iszero(n)\n",
    "        error(\"Dividing by zero not support\")\n",
    "    end\n",
    "    a = setrounding(T, RoundDown) do\n",
    "        # TODO: lower bound\n",
    "\n",
    "    end\n",
    "    b = setrounding(T, RoundUp) do\n",
    "        # TODO: upper bound\n",
    "\n",
    "    end\n",
    "    Interval(a, b)\n",
    "end\n",
    "\n",
    "@test_broken Interval(1.0,2.0)/3 ≡ Interval(0.3333333333333333, 0.6666666666666667)\n",
    "@test_broken Interval(1.0,2.0)/(-3) ≡ Interval(-0.6666666666666667, -0.3333333333333333)\n",
    "\n",
    "# Following should implement ⊗\n",
    "function *(x::Interval, y::Interval)\n",
    "    T = promote_type(typeof(x.a), typeof(x.b))\n",
    "    if 0 in x || 0 in y\n",
    "        error(\"Multiplying with intervals containing 0 not supported.\")\n",
    "    end\n",
    "    if x.a > x.b || y.a > y.b\n",
    "        error(\"Empty intervals not supported.\")\n",
    "    end\n",
    "    a = setrounding(T, RoundDown) do\n",
    "        # TODO: lower bound\n",
    "\n",
    "    end\n",
    "    b = setrounding(T, RoundUp) do\n",
    "        # TODO: upper bound\n",
    "\n",
    "    end\n",
    "    Interval(a, b)\n",
    "end\n",
    "\n",
    "@test_broken Interval(1.1, 1.2) * Interval(2.1, 3.1) ≡ Interval(2.31, 3.72)\n",
    "@test_broken Interval(-1.2, -1.1) * Interval(2.1, 3.1) ≡ Interval(-3.72, -2.31)\n",
    "@test_broken Interval(1.1, 1.2) * Interval(-3.1, -2.1) ≡ Interval(-3.72, -2.31)\n",
    "@test_broken Interval(-1.2, -1.1) * Interval(-3.1, -2.1) ≡ Interval(2.31, 3.72)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "-----"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The following function  computes the first `n+1` terms of the Taylor series of $\\exp(x)$:\n",
    "$$\n",
    "\\sum_{k=0}^n {x^k \\over k!}\n",
    "$$\n",
    "(similar to the one seen in lectures)."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "exp_t (generic function with 1 method)"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "function exp_t(x, n)\n",
    "    ret = one(x) # 1 of same type as x\n",
    "    s = one(x)\n",
    "    for k = 1:n\n",
    "        s = s/k * x\n",
    "        ret = ret + s\n",
    "    end\n",
    "    ret\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "-----"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Problem 5.3⋆** Bound the tail of the Taylor series for ${\\rm e}^x$ assuming $|x| ≤ 1$.\n",
    "(Hint: ${\\rm e}^x ≤ 3$ for $x ≤ 1$.)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "------\n",
    "\n",
    "**Problem 5.4** Use the bound\n",
    "to write a function `exp_bound` which computes ${\\rm e}^x$ with rigorous error bounds, that is\n",
    "so that when applied to an interval $[a,b]$ it returns an interval that is\n",
    "guaranteed to contain the interval $[{\\rm e}^a, {\\rm e}^b]$."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "\u001b[33m\u001b[1mTest Broken\u001b[22m\u001b[39m\n  Expression: e_int.b - e_int.a ≤ 1.0e-13"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "#\n",
    "\n",
    "function exp_bound(x::Interval, n)\n",
    "    # TODO: Return an Interval such that exp(x) is guaranteed to be a subset\n",
    "\n",
    "end\n",
    "\n",
    "e_int = exp_bound(Interval(1.0,1.0), 20)\n",
    "@test_broken exp(big(1)) in e_int\n",
    "@test_broken exp(big(-1)) in exp_bound(Interval(-1.0,-1.0), 20)\n",
    "@test_broken e_int.b - e_int.a ≤ 1E-13 # we want our bounds to be sharp"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "------\n",
    "**Problem 5.5** Use `big` and `setprecision` to compute ℯ to a 1000 decimal digits with\n",
    "rigorous error bounds."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "#"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.3-pre.0"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.3-pre.0",
   "language": "julia"
  }
 },
 "nbformat": 4
}