Calculating pi rigorously with floating-point arithmetic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rigorous bounds on $\\pi$ with floating-point arithmetic \n",
    "\n",
    "[David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders/)\n",
    "\n",
    "Department of Physics, Faculty of Sciences, National University of Mexico (UNAM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the book [*Validated Numerics*](http://press.princeton.edu/titles/9488.html) (Princeton, 2011) by [Warwick Tucker](http://www2.math.uu.se/~warwick/CAPA/warwick/warwick.html), we find *rigorous* (i.e., guaranteed, or *validated*) bounds on $\\pi$ using floating-point arithmetic and directed rounding.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the sum\n",
    "\n",
    "$$ S := \\sum_{n=1}^\\infty \\frac{1}{n^2}$$.\n",
    "\n",
    "It is [known](http://en.wikipedia.org/wiki/Basel_problem) that the exact value is $S = \\frac{\\pi^2}{6}$.\n",
    "Thus, if we can calculate rigorous bounds on $S$, then we can find rigorous bounds on $\\pi$. \n",
    "\n",
    "The idea is to split $S$ up into two parts, $S = S_N + T_N$, with\n",
    "$ S_N := \\sum_{n=1}^N \\frac{1}{n^2}$\n",
    "and $T_N := S - S_N = \\sum_{n=N+1}^\\infty n^{-2}$.\n",
    "We will evalute $S_N$ numerically and find an analytical bound for $T_N$.\n",
    "\n",
    "By bounding $T_N$ using integrals from below and above, we can see that $\\frac{1}{N+1} \\le T_N \\le \\frac{1}{N}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$S_N$ may be calculated easily by summing either forwards or backwards:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "reverse_sum (generic function with 2 methods)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function forward_sum(N, T=Float64)\n",
    "    total = zero(T)\n",
    "\n",
    "    for i in 1:N\n",
    "        total += one(T) / (i^2)\n",
    "    end\n",
    "\n",
    "    total\n",
    "end\n",
    "\n",
    "function reverse_sum(N, T=Float64)\n",
    "    total = zero(T)\n",
    "\n",
    "    for i in N:-1:1\n",
    "        total += one(T) / (i^2)\n",
    "    end\n",
    "    \n",
    "    total\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find *rigorous* bounds for $S_N$, we use \"directed rounding\", that is, we round downwards for the lower bound and  upwards for the upper bound:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.644933066959796"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6\n",
    "\n",
    "lowerbound_S_N = \n",
    "    with_rounding(Float64, RoundDown) do\n",
    "        forward_sum(N)\n",
    "    end\n",
    "\n",
    "upperbound_S_N = \n",
    "    with_rounding(Float64, RoundUp) do\n",
    "        forward_sum(N)\n",
    "        \n",
    "    end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We incorporate the respective bound on $T_N$ to obtain the bounds on $S$, and hence on $\\pi$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.1415926534833463,3.1415926536963346)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6\n",
    "\n",
    "lower_pi =\n",
    "    with_rounding(Float64, RoundDown) do\n",
    "        lower_bound = forward_sum(N) + 1/(N+1)\n",
    "        sqrt(6 * lower_bound)\n",
    "    end\n",
    "\n",
    "upper_pi = \n",
    "    with_rounding(Float64, RoundUp) do\n",
    "        upper_bound = forward_sum(N) + 1/N\n",
    "        sqrt(6 * upper_bound)\n",
    "    end\n",
    "\n",
    "lower_pi, upper_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may check that the true value of $\\pi$ is indeed contained in the interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lower_pi < big(pi) < upper_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may repeat the calculation summing in the opposite direction for a more accurate answer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.1415926535893144,3.141592653590272,true)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6\n",
    "\n",
    "lower_pi =\n",
    "    with_rounding(Float64, RoundDown) do\n",
    "        lower_bound = reverse_sum(N) + 1/(N+1)\n",
    "        sqrt(6 * lower_bound)\n",
    "    end\n",
    "\n",
    "upper_pi = \n",
    "    with_rounding(Float64, RoundUp) do\n",
    "        upper_bound = reverse_sum(N) + 1/N\n",
    "        sqrt(6 * upper_bound)\n",
    "    end\n",
    "\n",
    "lower_pi, upper_pi, lower_pi < big(pi) < upper_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `BigFloat`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is not clear (to me) if the standard `sqrt` operation on `Float64`s respects the rounding mode. We may use `BigFloat`s  to guarantee a correctly-rounded `sqrt` function, although the calculation is much slower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.1415926535893144e+00 with 53 bits of precision,3.1415926535902718e+00 with 53 bits of precision,true)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set_bigfloat_precision(53)  # same precision as `Float64`\n",
    "\n",
    "N = 10^6\n",
    "\n",
    "lower_pi =\n",
    "    with_rounding(BigFloat, RoundDown) do\n",
    "    lower_bound = reverse_sum(N, BigFloat) + one(BigFloat)/(N+1)\n",
    "        sqrt(6 * lower_bound)\n",
    "    end\n",
    "\n",
    "upper_pi = \n",
    "    with_rounding(BigFloat, RoundUp) do\n",
    "        upper_bound = reverse_sum(N, BigFloat) + one(BigFloat)/N\n",
    "        sqrt(6 * upper_bound)\n",
    "    end\n",
    "\n",
    "lower_pi, upper_pi, lower_pi < big(pi) < upper_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note that the result seems to be the same as the result with `Float64`, so perhaps the `sqrt` function does respect the rounding mode.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In principal, we could attain arbitrarily good precision with higher-precision `BigFloat`s, but the result is hampered by the slow convergence of the series."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.3.7-pre",
   "language": "julia",
   "name": "julia 0.3"
  },
  "language_info": {
   "name": "julia",
   "version": "0.3.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}