Extended Precision.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python for Everyone!<br/>[Oregon Curriculum Network](http://4dsolutions.net/ocn/)\n",
    "\n",
    "## Extended Precision with the Native Decimal Type\n",
    "### With LaTeX and Generator Functions\n",
    "\n",
    "<img src=\"https://c8.staticflickr.com/6/5691/30269841575_8bea763a54.jpg\" alt=\"TAOCP\" \n",
    "style=\"width: 50%; height: 50%\"/>\n",
    "\n",
    "The Python Standard Library provides a *decimal* module containing the class Decimal.  A Decimal object behaves according to the base 10 algorithms we learn in school.  \n",
    "\n",
    "The precision i.e. number of decimal places to which computations are carried out, is set globally, or within the scope of a context manager.\n",
    "\n",
    "Note also that Jupyter Notebooks are able to render LaTeX when commanded to do so with the %%latex magic command.  As a first example, here is an expression for the mathematical constant *e*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "e = lim_{n \\to \\infty} (1 + 1/n)^n\n",
       "\\end{align}"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "\\begin{align}\n",
    "e = lim_{n \\to \\infty} (1 + 1/n)^n\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.718281828459045\n",
      "3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "from math import e, pi\n",
    "print(e)  # as a floating point number\n",
    "print(pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets show setting precision to a thousand places within a scope defined by decimal.localcontext.  We set precision internally to the scope.  By default, precision is to 28 places.  We set n to 1 followed by 102 zeros, so a very large number.  The resulting computation matches a published value for *e* to 100 decimal places."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import decimal\n",
    "\n",
    "with decimal.localcontext() as ctx:  # context manager\n",
    "    ctx.prec = 1000\n",
    "    n = decimal.Decimal(1e102)\n",
    "    e = (1 + 1/n) ** n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "e_1000_places = \"\"\"2.7182818284590452353602874713526624977572470936999595749669\n",
    "6762772407663035354759457138217852516642742746639193200305992181741359662904357\n",
    "2900334295260595630738132328627943490763233829880753195251019011573834187930702\n",
    "1540891499348841675092447614606680822648001684774118537423454424371075390777449\n",
    "9206955170276183860626133138458300075204493382656029760673711320070932870912744\n",
    "3747047230696977209310141692836819025515108657463772111252389784425056953696770\n",
    "7854499699679468644549059879316368892300987931277361782154249992295763514822082\n",
    "6989519366803318252886939849646510582093923982948879332036250944311730123819706\n",
    "8416140397019837679320683282376464804295311802328782509819455815301756717361332\n",
    "0698112509961818815930416903515988885193458072738667385894228792284998920868058\n",
    "2574927961048419844436346324496848756023362482704197862320900216099023530436994\n",
    "1849146314093431738143640546253152096183690888707016768396424378140592714563549\n",
    "0613031072085103837505101157477041718986106873969655212671546889570350354\"\"\"\n",
    "e_1000_places = e_1000_places.replace(\"\\n\",str())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str(e)[2:103] == e_1000_places[2:103]  # skipping \"2.\" and going to 100 decimals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the context below, the value of n starts at 10 and then gets two more zeros every time around the loop.  \n",
    "\n",
    "The *yield* keyword is similar to *return* in handing back an object, however a generator function then pauses to pick up where it left off when nudged by next(), (which triggers \\_\\_next\\_\\_ internally).\n",
    "\n",
    "Generator functions do not forget their internal state as they advance through next values.\n",
    "\n",
    "Note that when a Decimal type object operates with an integer, that integer is coerced (cast) as a Decimal object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with decimal.localcontext() as ctx:  # context manager\n",
    "    ctx.prec = 1000\n",
    "    def converge():  # generator function\n",
    "        n = decimal.Decimal('10')\n",
    "        while True:\n",
    "            yield (1 + 1/n) ** n\n",
    "            n = n * 100 # two more zeros\n",
    "    f = converge()\n",
    "    for _ in range(9):\n",
    "        next(f) # f.__next__() <--- not quite like Python 2.x (f.next())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('2.718281828459045235224373380')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = next(f)\n",
    "r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str(r)[:20] == e_1000_places[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://sciencenode.org/img/img_2012/stamp.JPG\" alt=\"Ramanujan Postage Stamp\" \n",
    "style=\"width: 50%; height: 50%\"/>\n",
    "\n",
    "The fancier LaTeX below renders a famous equation by Ramanujan, which has been shown to converge to 1/π and therefore π very quickly, relative to many other algorithms.  I don't think anyone understands how some random guy could think up such a miraculaous equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "\\frac{1}{\\pi} = \\frac{2\\sqrt{2}}{9801} \\sum^\\infty_{k=0} \\frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}\n",
       "\\end{align}                                                     "
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "\\begin{align}\n",
    "\\frac{1}{\\pi} = \\frac{2\\sqrt{2}}{9801} \\sum^\\infty_{k=0} \\frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}\n",
    "\\end{align}                                                     "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Lets break the above equation into three components, a constant A and a product of Bn and Cn that accumulate as a sum.  The latter two may be written as generators, with n increasing by one each time they're fed to next(). An outermost generator ties them all together and takes the reciprocal to get an approximation for π itself. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align*}\n",
       "\\frac{1}{\\pi} = \\frac{2\\sqrt{2}}{9801}\\sum_{n = 0}^{\\infty}\\frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}99^{4n}} = A\\sum_{n = 0}^{\\infty}B_{n}C_{n},\\\\\n",
       "A = \\frac{2\\sqrt{2}}{9801},\\,\\\\ \n",
       "B_{n} = \\frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}},\\,\\\\\n",
       "C_{n} = \\frac{1}{99^{4n}}\n",
       "\\end{align*}"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "\\begin{align*}\n",
    "\\frac{1}{\\pi} = \\frac{2\\sqrt{2}}{9801}\\sum_{n = 0}^{\\infty}\\frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}99^{4n}} = A\\sum_{n = 0}^{\\infty}B_{n}C_{n},\\\\\n",
    "A = \\frac{2\\sqrt{2}}{9801},\\,\\\\ \n",
    "B_{n} = \\frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}},\\,\\\\\n",
    "C_{n} = \\frac{1}{99^{4n}}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('0.0002885855652225477091728780173879600201142070962915923014338699597981292681281720311907739076273118198')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from math import factorial\n",
    "from decimal import Decimal as D\n",
    "decimal.getcontext().prec=100\n",
    "A = (2 * D('2').sqrt()) / 9801\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def B():\n",
    "    n = 0\n",
    "    while True:\n",
    "        numerator = factorial(4 * n) * (D(1103) + 26390 * n) \n",
    "        denominator = (4 ** (4*n))*(factorial(n))**4\n",
    "        yield numerator / denominator\n",
    "        n += 1\n",
    "\n",
    "def C():\n",
    "    n = 0\n",
    "    while True:\n",
    "        yield 1 / (D('99')**(4*n))\n",
    "        n += 1\n",
    "        \n",
    "def Pi():\n",
    "    Bn = B()\n",
    "    Cn = C()\n",
    "    the_sum = 0\n",
    "    while True:\n",
    "        the_sum += next(Bn) * next(Cn) \n",
    "        yield 1/(A * the_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pi = Pi()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.141592653589793238462649065702758898156677480462334781168399595644739794558841580205059234965983146')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next(pi)\n",
    "next(pi)\n",
    "next(pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pi_1000_places = \"\"\"3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944\n",
    "5923078164 0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172\n",
    "5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461\n",
    "2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273\n",
    "7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652\n",
    "1384146951 9415116094 3305727036 5759591953 0921861173 8193261179 3105118548 0744623799\n",
    "6274956735 1885752724 8912279381 8301194912 9833673362 4406566430 8602139494 6395224737\n",
    "1907021798 6094370277 0539217176 2931767523 8467481846 7669405132 0005681271 4526356082\n",
    "7785771342 7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235\n",
    "4201995611 2129021960 8640344181 5981362977 4771309960 5187072113 4999999837 2978049951\n",
    "0597317328 1609631859 5024459455 3469083026 4252230825 3344685035 2619311881 7101000313\n",
    "7838752886 5875332083 8142061717 7669147303 5982534904 2875546873 1159562863 8823537875\n",
    "9375195778 1857780532 1712268066 1300192787 6611195909 2164201989\"\"\"\n",
    "pi_1000_places = pi_1000_places.replace(\" \",\"\").replace(\"\\n\",\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r = next(pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str(r)[:20] == pi_1000_places[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Resource*:\n",
    "\n",
    "[Ballad about Ramanujan](https://archive.org/details/Ramanujan) by Mark Engelberg, 2005."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This I-Python Notebook is by Kirby Urner, copyleft MIT License, October 2016."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}