Py_Pi_Day_2017.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#Python-Pi-Day-2017\" data-toc-modified-id=\"Python-Pi-Day-2017-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Python Pi Day 2017</a></div><div class=\"lev1 toc-item\"><a href=\"#Computing-a-lot-of-digits-of-$\\pi$?\" data-toc-modified-id=\"Computing-a-lot-of-digits-of-$\\pi$?-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Computing a lot of digits of <span class=\"MathJax_Preview\" style=\"color: inherit;\"></span><span class=\"MathJax\" id=\"MathJax-Element-388-Frame\" tabindex=\"0\" style=\"position: relative;\" data-mathml=\"<math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;><mi>&amp;#x03C0;</mi></math>\" role=\"presentation\"><nobr aria-hidden=\"true\"><span class=\"math\" id=\"MathJax-Span-2745\" role=\"math\" style=\"width: 0.702em; display: inline-block;\"><span style=\"display: inline-block; position: relative; width: 0.554em; height: 0px; font-size: 125%;\"><span style=\"position: absolute; clip: rect(1.941em, 1000.55em, 2.572em, -1000em); top: -2.462em; left: 0em;\"><span class=\"mrow\" id=\"MathJax-Span-2746\"><span class=\"mi\" id=\"MathJax-Span-2747\" style=\"font-family: STIXMathJax_Main; font-style: italic;\">π<span style=\"display: inline-block; overflow: hidden; height: 1px; width: 0.032em;\"></span></span></span><span style=\"display: inline-block; width: 0px; height: 2.462em;\"></span></span></span><span style=\"display: inline-block; overflow: hidden; vertical-align: -0.061em; border-left: 0px solid; width: 0px; height: 0.634em;\"></span></span></nobr><span class=\"MJX_Assistive_MathML\" role=\"presentation\"><math xmlns=\"http://www.w3.org/1998/Math/MathML\"><mi>π</mi></math></span></span><script type=\"math/tex\" id=\"MathJax-Element-388\">\\pi</script>?</a></div><div class=\"lev2 toc-item\"><a href=\"#Two-simple-methods-for-finding-the-first-digits-of-$\\pi$\" data-toc-modified-id=\"Two-simple-methods-for-finding-the-first-digits-of-$\\pi$-21\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Two simple methods for finding the first digits of <span class=\"MathJax_Preview\" style=\"color: inherit;\"></span><span class=\"MathJax\" id=\"MathJax-Element-391-Frame\" tabindex=\"0\" style=\"position: relative;\" data-mathml=\"<math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;><mi>&amp;#x03C0;</mi></math>\" role=\"presentation\"><nobr aria-hidden=\"true\"><span class=\"math\" id=\"MathJax-Span-2754\" role=\"math\" style=\"width: 0.693em; display: inline-block;\"><span style=\"display: inline-block; position: relative; width: 0.541em; height: 0px; font-size: 126%;\"><span style=\"position: absolute; clip: rect(1.953em, 1000.54em, 2.615em, -1000em); top: -2.489em; left: 0em;\"><span class=\"mrow\" id=\"MathJax-Span-2755\"><span class=\"mi\" id=\"MathJax-Span-2756\" style=\"font-family: STIXMathJax_Main; font-style: italic;\">π<span style=\"display: inline-block; overflow: hidden; height: 1px; width: 0.032em;\"></span></span></span><span style=\"display: inline-block; width: 0px; height: 2.489em;\"></span></span></span><span style=\"display: inline-block; overflow: hidden; vertical-align: -0.068em; border-left: 0px solid; width: 0px; height: 0.653em;\"></span></span></nobr><span class=\"MJX_Assistive_MathML\" role=\"presentation\"><math xmlns=\"http://www.w3.org/1998/Math/MathML\"><mi>π</mi></math></span></span><script type=\"math/tex\" id=\"MathJax-Element-391\">\\pi</script></a></div><div class=\"lev3 toc-item\"><a href=\"#Fraction-approximations,-and-$\\pi$-imported-from-the-math-module\" data-toc-modified-id=\"Fraction-approximations,-and-$\\pi$-imported-from-the-math-module-211\"><span class=\"toc-item-num\">2.1.1&nbsp;&nbsp;</span>Fraction approximations, and <span class=\"MathJax_Preview\" style=\"color: inherit;\"></span><span class=\"MathJax\" id=\"MathJax-Element-392-Frame\" tabindex=\"0\" style=\"position: relative;\" data-mathml=\"<math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;><mi>&amp;#x03C0;</mi></math>\" role=\"presentation\"><nobr aria-hidden=\"true\"><span class=\"math\" id=\"MathJax-Span-2757\" role=\"math\" style=\"width: 0.736em; display: inline-block;\"><span style=\"display: inline-block; position: relative; width: 0.56em; height: 0px; font-size: 129%;\"><span style=\"position: absolute; clip: rect(1.898em, 1000.56em, 2.602em, -1000em); top: -2.455em; left: 0em;\"><span class=\"mrow\" id=\"MathJax-Span-2758\"><span class=\"mi\" id=\"MathJax-Span-2759\" style=\"font-family: STIXMathJax_Main; font-style: italic;\">π<span style=\"display: inline-block; overflow: hidden; height: 1px; width: 0.032em;\"></span></span></span><span style=\"display: inline-block; width: 0px; height: 2.455em;\"></span></span></span><span style=\"display: inline-block; overflow: hidden; vertical-align: -0.079em; border-left: 0px solid; width: 0px; height: 0.686em;\"></span></span></nobr><span class=\"MJX_Assistive_MathML\" role=\"presentation\"><math xmlns=\"http://www.w3.org/1998/Math/MathML\"><mi>π</mi></math></span></span><script type=\"math/tex\" id=\"MathJax-Element-392\">\\pi</script> imported from the <code>math</code> module</a></div><div class=\"lev3 toc-item\"><a href=\"#A-simple-Monte-Carlo-method\" data-toc-modified-id=\"A-simple-Monte-Carlo-method-212\"><span class=\"toc-item-num\">2.1.2&nbsp;&nbsp;</span>A simple <a href=\"https://en.wikipedia.org/wiki/Pi#Monte_Carlo_methods\" target=\"_blank\">Monte-Carlo method</a></a></div><div class=\"lev2 toc-item\"><a href=\"#$100$-first-digits-of-$\\pi$\" data-toc-modified-id=\"$100$-first-digits-of-$\\pi$-22\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span><span class=\"MathJax_Preview\" style=\"color: inherit;\"></span><span class=\"MathJax\" id=\"MathJax-Element-397-Frame\" tabindex=\"0\" style=\"position: relative;\" data-mathml=\"<math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;><mn>100</mn></math>\" role=\"presentation\"><nobr aria-hidden=\"true\"><span class=\"math\" id=\"MathJax-Span-2778\" role=\"math\" style=\"width: 1.92em; display: inline-block;\"><span style=\"display: inline-block; position: relative; width: 1.515em; height: 0px; font-size: 126%;\"><span style=\"position: absolute; clip: rect(1.885em, 1001.49em, 2.792em, -1000em); top: -2.67em; left: 0em;\"><span class=\"mrow\" id=\"MathJax-Span-2779\"><span class=\"mn\" id=\"MathJax-Span-2780\" style=\"font-family: STIXMathJax_Main;\">100</span></span><span style=\"display: inline-block; width: 0px; height: 2.67em;\"></span></span></span><span style=\"display: inline-block; overflow: hidden; vertical-align: -0.063em; border-left: 0px solid; width: 0px; height: 0.96em;\"></span></span></nobr><span class=\"MJX_Assistive_MathML\" role=\"presentation\"><math xmlns=\"http://www.w3.org/1998/Math/MathML\"><mn>100</mn></math></span></span><script type=\"math/tex\" id=\"MathJax-Element-397\">100</script> first digits of <span class=\"MathJax_Preview\" style=\"color: inherit;\"></span><span class=\"MathJax\" id=\"MathJax-Element-398-Frame\" tabindex=\"0\" style=\"position: relative;\" data-mathml=\"<math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;><mi>&amp;#x03C0;</mi></math>\" role=\"presentation\"><nobr aria-hidden=\"true\"><span class=\"math\" id=\"MathJax-Span-2781\" role=\"math\" style=\"width: 0.693em; display: inline-block;\"><span style=\"display: inline-block; position: relative; width: 0.541em; height: 0px; font-size: 126%;\"><span style=\"position: absolute; clip: rect(1.953em, 1000.54em, 2.615em, -1000em); top: -2.489em; left: 0em;\"><span class=\"mrow\" id=\"MathJax-Span-2782\"><span class=\"mi\" id=\"MathJax-Span-2783\" style=\"font-family: STIXMathJax_Main; font-style: italic;\">π<span style=\"display: inline-block; overflow: hidden; height: 1px; width: 0.032em;\"></span></span></span><span style=\"display: inline-block; width: 0px; height: 2.489em;\"></span></span></span><span style=\"display: inline-block; overflow: hidden; vertical-align: -0.068em; border-left: 0px solid; width: 0px; height: 0.653em;\"></span></span></nobr><span class=\"MJX_Assistive_MathML\" role=\"presentation\"><math xmlns=\"http://www.w3.org/1998/Math/MathML\"><mi>π</mi></math></span></span><script type=\"math/tex\" id=\"MathJax-Element-398\">\\pi</script></a></div><div class=\"lev2 toc-item\"><a href=\"#Using-mpmath\" data-toc-modified-id=\"Using-mpmath-23\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Using <a href=\"http://mpmath.org/\" target=\"_blank\"><code>mpmath</code></a></a></div><div class=\"lev2 toc-item\"><a href=\"#The-Gauss–Legendre-iterative-algorithm\" data-toc-modified-id=\"The-Gauss–Legendre-iterative-algorithm-24\"><span class=\"toc-item-num\">2.4&nbsp;&nbsp;</span>The Gauss–Legendre iterative algorithm</a></div><div class=\"lev3 toc-item\"><a href=\"#Using-float-numbers\" data-toc-modified-id=\"Using-float-numbers-241\"><span class=\"toc-item-num\">2.4.1&nbsp;&nbsp;</span>Using float numbers</a></div><div class=\"lev3 toc-item\"><a href=\"#Using-decimal.Decimal-to-improve-precision\" data-toc-modified-id=\"Using-decimal.Decimal-to-improve-precision-242\"><span class=\"toc-item-num\">2.4.2&nbsp;&nbsp;</span>Using <code>decimal.Decimal</code> to improve precision</a></div><div class=\"lev2 toc-item\"><a href=\"#Methods-based-on-a-convergent-series\" data-toc-modified-id=\"Methods-based-on-a-convergent-series-25\"><span class=\"toc-item-num\">2.5&nbsp;&nbsp;</span>Methods based on a convergent series</a></div><div class=\"lev3 toc-item\"><a href=\"#A-Leibniz-formula-(easy):\" data-toc-modified-id=\"A-Leibniz-formula-(easy):-251\"><span class=\"toc-item-num\">2.5.1&nbsp;&nbsp;</span><a href=\"https://en.wikipedia.org/wiki/Leibniz_formula_for_pi\" target=\"_blank\">A Leibniz formula</a> (<em>easy</em>):</a></div><div class=\"lev3 toc-item\"><a href=\"#Bailey-Borwein-Plouffe-series-(medium):\" data-toc-modified-id=\"Bailey-Borwein-Plouffe-series-(medium):-252\"><span class=\"toc-item-num\">2.5.2&nbsp;&nbsp;</span><a href=\"https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula\" target=\"_blank\">Bailey-Borwein-Plouffe series</a> (<em>medium</em>):</a></div><div class=\"lev3 toc-item\"><a href=\"#Bellard's-formula-(hard):\" data-toc-modified-id=\"Bellard's-formula-(hard):-253\"><span class=\"toc-item-num\">2.5.3&nbsp;&nbsp;</span><a href=\"https://en.wikipedia.org/wiki/Bellard%27s_formula\" target=\"_blank\">Bellard's formula</a> (<em>hard</em>):</a></div><div class=\"lev3 toc-item\"><a href=\"#One-Ramanujan's-formula-(hard):\" data-toc-modified-id=\"One-Ramanujan's-formula-(hard):-254\"><span class=\"toc-item-num\">2.5.4&nbsp;&nbsp;</span>One <a href=\"https://en.wikipedia.org/wiki/Approximations_of_%CF%80#Efficient_methods\" target=\"_blank\">Ramanujan's formula</a> (<em>hard</em>):</a></div><div class=\"lev3 toc-item\"><a href=\"#Chudnovsky-brothers'-formula-(hard):\" data-toc-modified-id=\"Chudnovsky-brothers'-formula-(hard):-255\"><span class=\"toc-item-num\">2.5.5&nbsp;&nbsp;</span><a href=\"https://en.wikipedia.org/wiki/Chudnovsky_algorithm\" target=\"_blank\">Chudnovsky brothers' formula</a> (<em>hard</em>):</a></div><div class=\"lev2 toc-item\"><a href=\"#Other-methods\" data-toc-modified-id=\"Other-methods-26\"><span class=\"toc-item-num\">2.6&nbsp;&nbsp;</span>Other methods</a></div><div class=\"lev3 toc-item\"><a href=\"#Trigonometric-methods-(hard)\" data-toc-modified-id=\"Trigonometric-methods-(hard)-261\"><span class=\"toc-item-num\">2.6.1&nbsp;&nbsp;</span>Trigonometric methods (<em>hard</em>)</a></div><div class=\"lev4 toc-item\"><a href=\"#High-precision-arccot-computation\" data-toc-modified-id=\"High-precision-arccot-computation-2611\"><span class=\"toc-item-num\">2.6.1.1&nbsp;&nbsp;</span><a href=\"http://en.literateprograms.org/Pi_with_Machin%27s_formula_%28Python%29#High-precision_arccot_computation\" target=\"_blank\">High-precision arccot computation</a></a></div><div class=\"lev4 toc-item\"><a href=\"#Applying-Machin's-formula\" data-toc-modified-id=\"Applying-Machin's-formula-2612\"><span class=\"toc-item-num\">2.6.1.2&nbsp;&nbsp;</span>Applying Machin's formula</a></div><div class=\"lev4 toc-item\"><a href=\"#Trying-to-solve-my-question!\" data-toc-modified-id=\"Trying-to-solve-my-question!-2613\"><span class=\"toc-item-num\">2.6.1.3&nbsp;&nbsp;</span>Trying to solve my question!</a></div><div class=\"lev4 toc-item\"><a href=\"#Conclusion\" data-toc-modified-id=\"Conclusion-2614\"><span class=\"toc-item-num\">2.6.1.4&nbsp;&nbsp;</span>Conclusion</a></div><div class=\"lev3 toc-item\"><a href=\"#(hard)-Unbounded-Spigot-Algorithm\" data-toc-modified-id=\"(hard)-Unbounded-Spigot-Algorithm-262\"><span class=\"toc-item-num\">2.6.2&nbsp;&nbsp;</span>(<em>hard</em>) <a href=\"http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf\" target=\"_blank\">Unbounded Spigot Algorithm</a></a></div><div class=\"lev3 toc-item\"><a href=\"#(hard)-Borwein's-algorithm\" data-toc-modified-id=\"(hard)-Borwein's-algorithm-263\"><span class=\"toc-item-num\">2.6.3&nbsp;&nbsp;</span>(<em>hard</em>) <a href=\"https://en.wikipedia.org/wiki/Borwein%27s_algorithm#Nonic_convergence\" target=\"_blank\">Borwein's algorithm</a></a></div><div class=\"lev2 toc-item\"><a href=\"#Examples-and-references\" data-toc-modified-id=\"Examples-and-references-27\"><span class=\"toc-item-num\">2.7&nbsp;&nbsp;</span>Examples and references</a></div><div class=\"lev3 toc-item\"><a href=\"#Links\" data-toc-modified-id=\"Links-271\"><span class=\"toc-item-num\">2.7.1&nbsp;&nbsp;</span>Links</a></div><div class=\"lev3 toc-item\"><a href=\"#Pie-!\" data-toc-modified-id=\"Pie-!-272\"><span class=\"toc-item-num\">2.7.2&nbsp;&nbsp;</span>Pie !</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Python Pi Day 2017\n",
    "> This is heavily inspired by what I did two years ago, see this page [cs101/hackhaton/14_03/2015](http://perso.crans.org/besson/cs101/hackathon/14_03_2015) on my website.\n",
    "\n",
    "Today is [Pi Day 2017](http://www.piday.org/), the day celebrating the number $\\pi$.\n",
    "For more details on this number, see [this Wikipedia page](https://en.wikipedia.org/wiki/Pi).\n",
    "\n",
    "----\n",
    "\n",
    "Let us use this occasion to showcase a few different approaches to compute the digits of the number $\\pi$.\n",
    "I will use the [Python](https://www.python.org/) programming language, and you are reading a [Jupyter notebook](https://www.jupyter.org/).\n",
    "\n",
    "[![made-with-jupyter](https://img.shields.io/badge/Made%20with-Jupyter-1f425f.svg)](http://jupyter.org/)[![made-with-python](https://img.shields.io/badge/Made%20with-Python-1f425f.svg)](https://www.python.org/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing a lot of digits of $\\pi$?\n",
    "\n",
    "- **What to do ?** I will present and implement some methods that can be used to compute the first digits of the irrational number $\\pi$.\n",
    "- **How ?** One method is based on random numbers, but all the other are simple (or not so simple) iterative algorithm: the more steps you compute, the more digits you will have!\n",
    "- **How to compare / assess the result ?** It is simple: the more digits you got, the better. We will also test the different methods implemented, and for the most efficient, see how fast it is to go up-to $140000$ digits.\n",
    "\n",
    "The simple goal is to write a *small* function that computes digits of pi, as fast as possible, and find the 10 digits from position 140317 to 140327!\n",
    "(that's the challenge I posted on Facebook)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## Two simple methods for finding the first digits of $\\pi$\n",
    "### Fraction approximations, and $\\pi$ imported from the `math` module"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Three approximations, using fractions, were known from a very long time (Aristote used $\\frac{355}{113}$ !).\n",
    "The first three approximations of pi are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " pi ~= 3.14 (two first digits).\n",
      " pi ~= 22/7 = 3.142857142857143 (two first digits).\n",
      " pi ~= 355/113 = 3.1415929203539825 (six first digits).\n"
     ]
    }
   ],
   "source": [
    "print(\" pi ~= 3.14 (two first digits).\")\n",
    "print(\" pi ~= 22/7 = {} (two first digits).\".format(22.0 / 7.0))\n",
    "print(\" pi ~= 355/113 = {} (six first digits).\".format(355.0 / 113.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method is extremely limited, and will not give you more than 13 correct digits, as `math.pi` is stored as a *float* number (limited precision)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First method: using math.pi gives pi ~= 3.14159265358979312 (17 digits are displayed here).\n"
     ]
    }
   ],
   "source": [
    "def mathpi():\n",
    "    from math import pi\n",
    "    return pi\n",
    "\n",
    "print(\"First method: using math.pi gives pi ~= {:.17f} (17 digits are displayed here).\".format(mathpi()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we know the digits, we can directly print them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first 100 digits of pi are 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679.\n"
     ]
    }
   ],
   "source": [
    "from decimal import Decimal\n",
    "bigpi = Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679')\n",
    "print(\"The first 100 digits of pi are {}.\".format(bigpi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simple [Monte-Carlo method](https://en.wikipedia.org/wiki/Pi#Monte_Carlo_methods)\n",
    "A simple Monte Carlo method for computing $\\pi$ is to draw a circle inscribed in a square, and randomly place dots in the square.\n",
    "The ratio of dots inside the circle to the total number of dots will approximately equal $\\pi / 4$, which is also the ratio of the area of the circle by the area of the square:\n",
    "\n",
    "![Example of a random simulation of this Monte-Carlo method](https://upload.wikimedia.org/wikipedia/commons/8/84/Pi_30K.gif \"Example of a random simulation of this Monte-Carlo method\")\n",
    "\n",
    "In Python, such simulation can basically be implemented like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import uniform\n",
    "\n",
    "def montecarlo_pi(nbPoints=10000):\n",
    "    \"\"\"Returns a probabilist estimate of pi, as a float number.\"\"\"\n",
    "    nbInside = 0\n",
    "    # we pick a certain number of points (nbPoints)\n",
    "    for i in range(nbPoints):\n",
    "        x = uniform(0, 1)\n",
    "        y = uniform(0, 1)\n",
    "        # (x, y) is now a random point in the square [0, 1] × [0, 1]\n",
    "        if (x**2 + y**2) < 1:\n",
    "            # This point (x, y) is inside the circle C(0, 1)\n",
    "            nbInside += 1\n",
    "\n",
    "    return 4 * float(nbInside) / floor(nbPoints)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The simple Monte-Carlo method with 10000 random points gave pi = 3.138\n"
     ]
    }
   ],
   "source": [
    "print(\"The simple Monte-Carlo method with {} random points gave pi = {}\".format(10000, montecarlo_pi(10000)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is an interesting method, but it is just too limited for approximating digits of $\\pi$.\n",
    "\n",
    "\n",
    "The previous two methods are quite limited, and not efficient enough if you want more than 13 digits (resp. 4 digits for the Monte-Carlo method)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $100$ first digits of $\\pi$\n",
    "$\\pi \\simeq 3.1415926535 ~ 8979323846 ~ 2643383279 ~ 5028841971 \\\\\\\\ 6939937510 ~ 5820974944 ~ 5923078164 ~ 9862803482 ~ 53421170679$ when computed to the first $100$ digits.\n",
    "\n",
    "Can you compute up to $1000$ digits? Up to $10000$ digits? Up to $100000$ digits? **Up to 1 million digits?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using [`mpmath`](http://mpmath.org/)\n",
    "This is a simple and lazy method, using the [`mpmath`](http://mpmath.org/) module.\n",
    "MPmath is a Python library for arbitrary-precision floating-point arithmetic (Multi-Precision), and it has a builtin highly-optimized algorithm to compute digits of $\\pi$.\n",
    "\n",
    "It works really fine up-to 1000000 digits (56 ms), from 1 million digits to be printed, printing them starts to get too time consuming (the IDE or the system might freeze)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import mpmath\n",
    "# from sympy import mpmath  # on older sympy versions\n",
    "mp = mpmath.mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can [arbitrarily set the precision](http://docs.sympy.org/dev/modules/mpmath/basics.html#setting-the-precision), with the constant `mp.dps` (digit numbers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A lazy method using the mpmath module:\n",
      "pi is approximatly 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198 (with 1000 digits).\n"
     ]
    }
   ],
   "source": [
    "mp.dps = 1000  # number of digits\n",
    "my_pi = mp.pi  # Gives pi to a thousand places\n",
    "print(\"A lazy method using the mpmath module:\\npi is approximatly {} (with {} digits).\".format(my_pi, mp.dps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let save it for further comparison of simpler methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100001"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.dps = 100000  # number of digits\n",
    "len(str(mp.pi))\n",
    "mpmath_pi = Decimal(str(mp.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve the initial challenge easily:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9341076406\n"
     ]
    }
   ],
   "source": [
    "mp.dps = 140330\n",
    "print(str(mp.pi)[2:][140317:140317+10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And it will most probably be the quickest method presented here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9341076406\n",
      "9341076406\n",
      "9341076406\n",
      "9341076406\n",
      "1 loop, best of 3: 236 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit mp.dps=140330;print(str(mp.pi)[2:][140317:140317+10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Asking for $10$ times more digits take about $100$ more of time (that's a bad news)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7966916520\n",
      "7966916520\n",
      "7966916520\n",
      "7966916520\n",
      "The slowest run took 6.34 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "1 loop, best of 3: 20.4 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit mp.dps=1403230;print(str(mp.pi)[2:][1403217:1403217+10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Gauss–Legendre iterative algorithm\n",
    "> More details can be found on [this page](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm).\n",
    "\n",
    "The first method given here is explained in detail.\n",
    "This algorithm is called the [Gauss-Legendre method](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm), and for example it was used to compute the first $206 158 430 000$ decimal digits of $\\pi$ on September 18th to 20th, $1999$.\n",
    "\n",
    "It is a very simply **iterative algorithm** (ie. step by step, you update the variables, as long as you want):\n",
    "\n",
    "1. First, start with 4 variables $a_0, b_0, t_0, p_0$, and their initial values are $a_0 = 1, b_0 = 1/\\sqrt{2}, t_0 = 1/4, p_0 = 1$.\n",
    "\n",
    "2. Then, update the variables (or create 4 new ones $a_{n+1}, b_{n+1}, t_{n+1}, p_{n+1}$) by repeating the following instructions (in this order) until the difference of $a_{n}$ and  $b_{n}$, is within the desired accuracy:\n",
    "   - $a_{n+1} = \\frac{a_n + b_n}{2}$,\n",
    "   - $b_{n+1} = \\sqrt{a_n \\times b_n}$,\n",
    "   - $t_{n+1} = t_n - p_n (a_n - a_{n+1})^2$,\n",
    "   - $p_{n+1} = 2 p_n$.\n",
    "\n",
    "3. Finally, the desired approximation of $\\pi$ is: $$\\pi \\simeq \\frac{(a_{n+1} + b_{n+1})^2}{4 t_{n+1}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using float numbers\n",
    "The first three iterations give (approximations given up to and including the first incorrect digit):\n",
    "\n",
    "    3.140 …\n",
    "    3.14159264 …\n",
    "    3.1415926535897932382 …\n",
    "\n",
    "The algorithm has **second-order convergent nature**, which essentially means that the number of correct digits doubles with each step of the algorithm.\n",
    "Try to see how far it can go in less than 120 seconds (print the approximation of $\\pi$ after every step, and stop/kill the program after 2 minutes).\n",
    "\n",
    "> This method is based on [computing the limits of the arithmetic–geometric mean](https://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean) of some well-chosen number ([see on Wikipédia for more mathematical details](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm#Mathematical_background))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def gauss_legendre_1(max_step):\n",
    "    \"\"\"Float number implementation of the Gauss-Legendre algorithm, for max_step steps.\"\"\"\n",
    "    a = 1.\n",
    "    b = 1./math.sqrt(2)\n",
    "    t = 1./4.0\n",
    "    p = 1.\n",
    "    for i in range(max_step):\n",
    "        at = (a + b) / 2.0\n",
    "        bt = math.sqrt(a*b)\n",
    "        tt = t - p*(a-at)**2\n",
    "        pt = 2.0 * p\n",
    "        a, b, t, p = at, bt, tt, pt\n",
    "\n",
    "    my_pi = ((a+b)**2)/(4.0*t)\n",
    "    return my_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589794"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi is approximately: 3.141592653589794 (as a float number, precision is limited).\n",
      "Accuracy % with math.pi: 2.827e-14\n",
      "Accuracy % with mpmath_pi: 2.827e-14\n"
     ]
    }
   ],
   "source": [
    "my_pi = gauss_legendre_1(4)\n",
    "my_pi\n",
    "print(\"pi is approximately: {:.15f} (as a float number, precision is limited).\".format(my_pi))\n",
    "accuracy = 100*abs(math.pi - my_pi)/math.pi\n",
    "print(\"Accuracy % with math.pi: {:.4g}\".format(accuracy))\n",
    "accuracy = 100*abs(float(mpmath_pi) - my_pi)/float(mpmath_pi)\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.141592653589794"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi is approximately: 3.141592653589794 (as a float number, precision is limited).\n",
      "Accuracy % with math.pi: 2.827e-14\n",
      "Accuracy % with mpmath_pi: 2.827e-14\n"
     ]
    }
   ],
   "source": [
    "my_pi = gauss_legendre_1(40)\n",
    "my_pi\n",
    "print(\"pi is approximately: {:.15f} (as a float number, precision is limited).\".format(my_pi))\n",
    "accuracy = 100*abs(math.pi - my_pi)/math.pi\n",
    "print(\"Accuracy % with math.pi: {:.4g}\".format(accuracy))\n",
    "accuracy = 100*abs(float(mpmath_pi) - my_pi)/float(mpmath_pi)\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This first implementation of the Gauss-Legendre algorithm is limited to a precision of 13 or 14 digits. But it converges quickly ! (4 steps here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `decimal.Decimal` to improve precision\n",
    "The limitation of this first algorithm came from using simple *float* numbers.\n",
    "Let us now use [`Decimal`](https://docs.python.org/3/library/decimal.html) numbers to keep as many digits after the decimal as we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from decimal import Decimal, getcontext"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gauss_legendre_2(max_step):\n",
    "    \"\"\"Decimal number implementation of the Gauss-Legendre algorithm, for max_step steps.\"\"\"\n",
    "    # trick to improve precision\n",
    "    getcontext().prec = 3 + 2**(max_step + 2)\n",
    "\n",
    "    cst_2 = Decimal(2.0)\n",
    "    cst_4 = Decimal(4.0)\n",
    "    a = Decimal(1.0)\n",
    "    b = Decimal(0.5).sqrt()\n",
    "    t = Decimal(0.25)\n",
    "    p = Decimal(1.0)\n",
    "\n",
    "    for i in range(max_step):\n",
    "        new_a = (a+b)/cst_2\n",
    "        new_b = (a*b).sqrt()\n",
    "        new_t = Decimal(t - p*(a - new_a)**2)\n",
    "        new_p = cst_2*p\n",
    "\n",
    "        a, b, t, p = new_a, new_b, new_t, new_p\n",
    "\n",
    "    my_pi = Decimal(((a+b)**2)/(cst_4*t))\n",
    "    return my_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pi is approximately: 3.14159265358979323846264338327950288419716939937510582097494459.\n",
      "Accuracy % with math.pi: 3.898e-15\n",
      "Accuracy % with mpmath_pi: 7.659e-83\n"
     ]
    }
   ],
   "source": [
    "my_pi = gauss_legendre_2(5)\n",
    "print(\"pi is approximately: {}.\".format(my_pi.to_eng_string()[:2**(5+1)]))\n",
    "\n",
    "accuracy = 100*abs(Decimal(math.pi) - my_pi)/Decimal(math.pi)\n",
    "print(\"Accuracy % with math.pi: {:.4g}\".format(accuracy))\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second implementation of the Gauss-Legendre algorithm is now working better (when we adapt the precision). And it converges quickly ! (8 steps give a precision upto the 697th digits).\n",
    "\n",
    "How much did we lost in term of performance by using decimal numbers? About a factor $1000$, but that's normal as the second solution stores a lot of digits. They don't even compute the same response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 6.47 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "100000 loops, best of 3: 4.06 µs per loop\n",
      "100 loops, best of 3: 3.9 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit gauss_legendre_1(8)\n",
    "%timeit gauss_legendre_2(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Methods based on a convergent series\n",
    "For the following formulae, you can try to write a program that computes the partial sum of the series, up to a certain number of term ($N \\geq 1$).\n",
    "Basically, the bigger the $N$, the more digits you get (but the longer the program will run).\n",
    "\n",
    "Some methods might be really efficient, only needing a few number of steps (a small $N$) for computing many digits.\n",
    "Try to implement and compare at least two of these methods.\n",
    "You can use the function `sum` (or `math.fsum`) to compute the sum, or a simple `while`/`for` loop.\n",
    "\n",
    "All these partial sums are written up to an integer $N \\geq 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [A Leibniz formula](https://en.wikipedia.org/wiki/Leibniz_formula_for_pi) (*easy*):\n",
    "It has a number of digits sub-linear in the number $N$ of terms in the sum: we need $10$ times more terms to win just one digit.\n",
    "$$\\pi \\simeq 4\\sum_{n=0}^{N} \\cfrac {(-1)^n}{2n+1}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def leibniz(max_step):\n",
    "    \"\"\" Computing an approximation of pi with Leibniz series.\"\"\"\n",
    "    my_pi = Decimal(0)\n",
    "    for k in range(max_step):\n",
    "        my_pi += Decimal((-1)**k) / Decimal(2*k+1)\n",
    "    return Decimal(4) * my_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1405926538397929257')"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 0.03183\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 20  # trick to improve precision\n",
    "my_pi = leibniz(1000)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1414926535900432389')"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 0.003183\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 20  # trick to improve precision\n",
    "my_pi = leibniz(10000)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This first formula is very inefficient!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Bailey-Borwein-Plouffe series](https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula) (*medium*):\n",
    "It also has a number of digits linear in the number $N$ of terms in the sum.\n",
    "$$\\pi \\simeq \\sum_{n = 1}^{N} \\left( \\frac{1}{16^{n}} \\left( \\frac{4}{8n+1} - \\frac{2}{8n+4} - \\frac{1}{8n+5} - \\frac{1}{8n+6} \\right) \\right). $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bbp(max_step):\n",
    "    \"\"\" Computing an approximation of pi with Bailey-Borwein-Plouffe series.\"\"\"\n",
    "    my_pi = Decimal(0)\n",
    "    for k in range(max_step):\n",
    "        my_pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6)))\n",
    "    return my_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897911465')"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 6.659e-14\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 20  # trick to improve precision\n",
    "my_pi = bbp(10)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's pretty impressive, in only $10$ steps!\n",
    "But that algorithm is highly dependent on the precision we ask, and on the number of terms in the sum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303817')"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 8.417e-198\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 200  # trick to improve precision\n",
    "my_pi = bbp(200)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119494')"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 8.600e-498\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 500  # trick to improve precision\n",
    "my_pi = bbp(500)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is, of course, slower than the optimized algorithm from `mpmath`.\n",
    "And it does not scale well with the precision:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 95 ms per loop\n",
      "1 loop, best of 3: 661 ms per loop\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 10 + 1000  # trick to improve precision\n",
    "%timeit bbp(1000)\n",
    "\n",
    "getcontext().prec = 10 + 2000  # trick to improve precision\n",
    "%timeit bbp(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Bellard's formula](https://en.wikipedia.org/wiki/Bellard%27s_formula) (*hard*):\n",
    "It is a more efficient formula.\n",
    "$$\\pi \\simeq \\frac1{2^6} \\sum_{n=0}^N \\frac{(-1)^n}{2^{10n}} \\, \\left(-\\frac{2^5}{4n+1} - \\frac1{4n+3} + \\frac{2^8}{10n+1} - \\frac{2^6}{10n+3} - \\frac{2^2}{10n+5} - \\frac{2^2}{10n+7} + \\frac1{10n+9} \\right). $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bellard(max_step):\n",
    "    \"\"\" Computing an approximation of pi with Bellard series.\"\"\"\n",
    "    my_pi = Decimal(0)\n",
    "    for k in range(max_step):\n",
    "        my_pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3))\n",
    "    return my_pi * Decimal(1.0/(2**6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.141592653589793238462643383279490036578')"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 4.090e-31\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 40  # trick to improve precision\n",
    "my_pi = bellard(10)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's pretty impressive, in only $10$ steps!\n",
    "But that algorithm is again highly dependent on the precision we ask, and on the number of terms in the sum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005611526979981607242018469874601593584399987639992795109821175268982496078980620087198769557504210588066903133631260918929384314932145237223424648852944318026933908108945450734428224399103245946394')"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 2.220e-604\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 800  # trick to improve precision\n",
    "my_pi = bellard(200)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is, of course, slower than the optimized algorithm from `mpmath`.\n",
    "And it does not scale well with the precision:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 378 ms per loop\n",
      "1 loop, best of 3: 2.91 s per loop\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 10 + 1000  # trick to improve precision\n",
    "%timeit bellard(1000)\n",
    "\n",
    "getcontext().prec = 10 + 2000  # trick to improve precision\n",
    "%timeit bellard(2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also slower than BBP formula."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One [Ramanujan's formula](https://en.wikipedia.org/wiki/Approximations_of_%CF%80#Efficient_methods) (*hard*):\n",
    "It is efficient but harder to compute easily with high precision, due to the factorial.\n",
    "But hopefully, the function `math.factorial` works with Python *integers*, of arbitrary size, so this won't be an issue.\n",
    "\n",
    "$$\\frac{1}{\\pi} \\simeq \\frac{2\\sqrt{2}}{9801} \\sum_{n=0}^N \\frac{(4n)!(1103+26390n)}{(n!)^4 396^{4n}}. $$\n",
    "\n",
    "*Remark:* This method and the next one compute approximation of $\\frac{1}{\\pi}$, not $\\pi$.\n",
    "\n",
    "This method has a quadratic precision (the number of correct digits is of the order $\\mathcal{O}(N^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import factorial\n",
    "\n",
    "def ramanujan(max_step):\n",
    "    \"\"\" Computing an approximation of pi with a Ramanujan's formula.\"\"\"\n",
    "    my_pi = Decimal(0)\n",
    "    d_1103 = Decimal(1103)\n",
    "    d_26390 = Decimal(26390)\n",
    "    d_396 = Decimal(396)\n",
    "    for k in range(max_step):\n",
    "        my_pi += ((Decimal(factorial(4 * k))) * (d_1103 + d_26390 * Decimal(k))) / ( (Decimal(factorial(k)))**4 * (d_396**(4*k)))\n",
    "    my_pi = my_pi * 2 * Decimal(2).sqrt() / Decimal(9801)\n",
    "    my_pi = my_pi**(-1)\n",
    "    return my_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.141592653589793238462643383279555273157')"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 1.668e-30\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 40  # trick to improve precision\n",
    "my_pi = ramanujan(4)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588249728368707403112841783926380916831222457885852997711270901024127742895829505936')"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 2.382e-318\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 400  # trick to improve precision\n",
    "my_pi = ramanujan(40)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029553211653449872027559602364806654991198818347977535663698074265425278625518184175746728909777727940092593679577412880114072520662452344498309284389026310583057432459663517241665671357564259409181627161036619398380266804323123256475817186916514676524019320976885611753635255011894984447618075599914209997043044425803288826839697424724950586712681051712136647565220069328838660790206792014356556715545181106432983263815336472301873472006826126207435983354628153181952126467499415321188698738323941189429')"
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 6.658e-1596\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 2000  # trick to improve precision\n",
    "my_pi = ramanujan(200)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1595$ correct digits with $200$ terms, that's quite good!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 50.1 ms per loop\n",
      "1 loop, best of 3: 478 ms per loop\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 10 + 2000  # trick to improve precision\n",
    "%timeit ramanujan(200)\n",
    "\n",
    "getcontext().prec = 10 + 5000  # trick to improve precision\n",
    "%timeit ramanujan(400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to answer my initial question, using this naive implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-128-f3c69c467c77>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'time'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'getcontext().prec = 140350  # trick to improve precision\\ni = 140317\\nmy_pi = ramanujan(10000)\\nprint(str(my_pi)[2:][i:i+10])\\n\\nmp.dps=140330\\nprint(str(mp.pi)[2:][i:i+10])'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2113\u001b[0m             \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_expand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2114\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2115\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2116\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2117\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-59>\u001b[0m in \u001b[0;36mtime\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    187\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 188\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    189\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    190\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/IPython/core/magics/execution.py\u001b[0m in \u001b[0;36mtime\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n\u001b[1;32m   1183\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1184\u001b[0m             \u001b[0mst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mclock2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1185\u001b[0;31m             \u001b[0mexec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mglob\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlocal_ns\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1186\u001b[0m             \u001b[0mend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mclock2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1187\u001b[0m             \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m<ipython-input-123-e46be6bb0995>\u001b[0m in \u001b[0;36mramanujan\u001b[0;34m(max_step)\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0md_396\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mDecimal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m396\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax_step\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m         \u001b[0mmy_pi\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mDecimal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfactorial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0md_1103\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0md_26390\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mDecimal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mDecimal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfactorial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m4\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0md_396\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     11\u001b[0m     \u001b[0mmy_pi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmy_pi\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mDecimal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mDecimal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m9801\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \u001b[0mmy_pi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmy_pi\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "%%time\n",
    "getcontext().prec = 140350  # trick to improve precision\n",
    "i = 140317\n",
    "my_pi = ramanujan(10000)\n",
    "print(str(my_pi)[2:][i:i+10])\n",
    "\n",
    "mp.dps=140330\n",
    "print(str(mp.pi)[2:][i:i+10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... It was too slow!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Chudnovsky brothers' formula](https://en.wikipedia.org/wiki/Chudnovsky_algorithm) (*hard*):\n",
    "$$\\frac{1}{\\pi} \\simeq 12 \\sum_{n=0}^N \\frac {(-1)^{n}(6n)!(545140134n+13591409)}{(3n)!(n!)^{3}640320^{{3n+3/2}}}. $$\n",
    "In 2015, the best method is still the Chudnovsky brothers's formula.\n",
    "\n",
    "> Be careful when you use these formulae, *check carefully* the constants you wrote (545140134 will work well, 545140135 might not work as well!).\n",
    "\n",
    "> This formula is used as the logo of the [French agrégation of Mathematics](https://en.wikipedia.org/wiki/Agr%C3%A9gation) [website `agreg.org`](http://agreg.org/) :\n",
    "> ![http://agreg.org/LogoAgreg.gif](http://agreg.org/LogoAgreg.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from math import factorial\n",
    "\n",
    "def chudnovsky(max_step):\n",
    "    \"\"\" Computing an approximation of pi with Bellard series.\"\"\"\n",
    "    my_pi = Decimal(0)\n",
    "    for k in range(max_step):\n",
    "        my_pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))\n",
    "    my_pi = my_pi * Decimal(10005).sqrt()/4270934400\n",
    "    my_pi = my_pi**(-1)\n",
    "    return my_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is very efficient, as Ramanujan's formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009946576407895126946839835259570982582262052248940772671947826848260147699090264013639443745530506820349625245174939965143142980919065925093722169646151570985838741059788595977297549893016175392846813826868386894277415599185592524595395943104997252468084598727364469584865383673622262609912460805124388439045124413654976278079771569143599770012961608944169486855584840635342207222582848864815845602850601684273945226746767889525213852254995466672782398645659611635488623057745649803559363456817432411251507606947945109659609402522887971089314566913686722874894056010150330861792868092087476091782493858900971490967598526136554978189312978482168299894872265880485756401427047755513237964145152374623436454285844479526586782105114135473573952311342716610213596953623144295248493718711014576540359027993440374200731057853906219838744780847852710405603214833731504646600861733152271806431614075966417958570002328166126409832551776784733662072536424251339850953855478119486626262729703543651855423452347092')"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 1.191e-2835\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 3000  # trick to improve precision\n",
    "my_pi = chudnovsky(200)\n",
    "my_pi\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It gets $2834$ correct numbers in $200$ steps!\n",
    "It is more efficient that Ramanujan's formula, but slower to compute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 3.947e-5672\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 6000  # trick to improve precision\n",
    "my_pi = chudnovsky(400)\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 216 ms per loop\n",
      "1 loop, best of 3: 1.93 s per loop\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 3000  # trick to improve precision\n",
    "%timeit chudnovsky(200)\n",
    "\n",
    "getcontext().prec = 6000  # trick to improve precision\n",
    "%timeit chudnovsky(400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "About $2$ seconds to find correctly the first $5671$ digits? That's slow! But hey, it's Python (dynamic typing etc)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## Other methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trigonometric methods (*hard*)\n",
    "Some methods are based on the $\\mathrm{arccot}$ or $\\arctan$ functions, and use the appropriate Taylor series to approximate these functions.\n",
    "The best example is [Machin's formula](http://en.literateprograms.org/Pi_with_Machin%27s_formula_%28Python%29):\n",
    "$$\\pi = 16 \\;\\mathrm{arccot}(5) - 4 \\;\\mathrm{arccot}(239).$$\n",
    "\n",
    "And we use the Taylor series to approximate $\\mathrm{arccot}(x)$:\n",
    "$$\\mathrm{arccot}(x) = \\frac{1}{x} - \\frac{1}{3x^3} + \\frac{1}{5x^5} - \\frac{1}{7x^7} + \\dots = \\sum_{n=0}^{+\\infty} \\frac{(-1)^n}{(2n+1) x^{2n+1}} .$$\n",
    "\n",
    "This method is also explained here with some details.\n",
    "In order to obtain $n$ digits, we will use *fixed-point* arithmetic to compute $\\pi \\times 10^n$ as a Python `long` integer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### [High-precision arccot computation](http://en.literateprograms.org/Pi_with_Machin%27s_formula_%28Python%29#High-precision_arccot_computation)\n",
    "To calculate $\\mathrm{arccot}$ of an argument $x$, we start by dividing the number $1$ (represented by $10^n$, which we provide as the argument `unity`) by $x$ to obtain the first term.\n",
    "\n",
    "We then repeatedly divide by $x^2$ and a counter value that runs over $3$, $5$, $7$ etc, to obtain each next term.\n",
    "The summation is stopped at the first zero `term`, which in this *fixed-point* representation corresponds to a real value less than $10^{-n}$:\n",
    "\n",
    "```python\n",
    "def arccot(x, unity):\n",
    "    xpower = unity / x\n",
    "    sum = xpower\n",
    "    n = 3\n",
    "    sign = -1\n",
    "    while True:\n",
    "        xpower = xpower / (x*x)\n",
    "        term = xpower / n\n",
    "        if term == 0:\n",
    "            break  # we are done\n",
    "        sum += sign * term\n",
    "        sign = -sign\n",
    "        n += 2\n",
    "    return sum\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapting it to use Decimal numbers is easy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def arccot(x, unity):\n",
    "    \"\"\"Compute arccot(x) with a certain level of precision.\"\"\"\n",
    "    x = Decimal(x)\n",
    "    unity = Decimal(unity)\n",
    "    mysum = xpower = unity / x\n",
    "    n = 3\n",
    "    sign = -1\n",
    "    while True:\n",
    "        xpower = xpower / (x*x)\n",
    "        term = xpower / n\n",
    "        if not term:\n",
    "            break\n",
    "        mysum += sign * term\n",
    "        sign = -sign  # we alternate the sign\n",
    "        n += 2\n",
    "    return mysum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Applying Machin's formula\n",
    "Finally, the main function uses Machin's formula to compute $\\pi$ using the necessary level of precision, by using this high precision $\\mathrm{arccot}$ function:\n",
    "$$\\pi = 16 \\;\\mathrm{arccot}(5) - 4 \\;\\mathrm{arccot}(239).$$\n",
    "\n",
    "```python\n",
    "def machin(digits):\n",
    "    unity = 10**(digits + 10)\n",
    "    pi = 4 * (4*arccot(5, unity) - arccot(239, unity))\n",
    "    return pi / unity\n",
    "```\n",
    "\n",
    "To avoid rounding errors in the result, we use 10 guard digits internally during the calculation.\n",
    "We may now reproduce the historical result obtained by [Machin in 1706](https://en.wikipedia.org/wiki/John_Machin)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def machin(digits):\n",
    "    \"\"\"Compute pi with Machin's formula, with precision at least digits.\"\"\"\n",
    "    unity = 10**(digits + 10)\n",
    "    my_pi = Decimal(4) * (Decimal(4)*arccot(5, unity) - arccot(239, unity))\n",
    "    return my_pi / Decimal(unity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy % with mpmath_pi: 1.501e-9996\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 10000  # trick to improve precision\n",
    "my_pi = machin(100)\n",
    "\n",
    "accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi\n",
    "print(\"Accuracy % with mpmath_pi: {:.4g}\".format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we got the first $9995$ digits correctly... in $45$ seconds.\n",
    "That will still be too slow to have the digits at position $130317$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 22.7 s per loop\n",
      "1 loop, best of 3: 44.5 s per loop\n"
     ]
    }
   ],
   "source": [
    "getcontext().prec = 5000  # trick to improve precision\n",
    "%timeit machin(50)\n",
    "\n",
    "getcontext().prec = 10000  # trick to improve precision\n",
    "%timeit machin(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer.\n",
    "Typical implementation will be in highly efficient compiled language; like C or maybe Julia.\n",
    "\n",
    "Many [Machin-like formulas](https://en.wikipedia.org/wiki/Machin-like_formula) also exist, like:\n",
    "$$\\pi = 4\\arctan\\left(\\frac{1}{2}\\right) + 4 \\arctan\\left(\\frac{1}{3}\\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Trying to solve my question!\n",
    "The important parameter to tune is not the precision given to `machin()` but the `Decimal.prec` precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7126218283\n",
      "7126218283\n",
      "CPU times: user 1min 3s, sys: 12 ms, total: 1min 3s\n",
      "Wall time: 1min 3s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "i = 14031\n",
    "getcontext().prec = i + 20  # trick to improve precision\n",
    "mp.dps = i + 20\n",
    "print(str(mp.pi)[2:][i:i+10])\n",
    "\n",
    "my_pi = machin(11)\n",
    "print(str(my_pi)[2:][i:i+10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9341076406\n",
      "9341076406\n",
      "CPU times: user 10min 40s, sys: 224 ms, total: 10min 40s\n",
      "Wall time: 10min 41s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "i = 140317\n",
    "getcontext().prec = i + 20  # trick to improve precision\n",
    "mp.dps = i + 20\n",
    "print(str(mp.pi)[2:][i:i+10])\n",
    "\n",
    "my_pi = machin(50)\n",
    "print(str(my_pi)[2:][i:i+10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It was too slow too! But at least it worked!\n",
    "\n",
    "My manual implementation of Machin's formula took about $10$ minutes, on an old laptop with $1$ core running Python 3.5, to find the $10$ digits of $\\pi$ at index $140317$.\n",
    "\n",
    "#### Conclusion\n",
    "$\\implies$ To the question \"What are the $10$ digits of $\\pi$ at index $140317..140326$\", the answer is $9341076406$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more, see http://fredrikj.net/blog/2011/03/100-mpmath-one-liners-for-pi/."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (*hard*) [Unbounded Spigot Algorithm](http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf)\n",
    "This algorithm is quite efficient, but not easy to explain. Go check on-line if you want.\n",
    "\n",
    "See this page (http://codepad.org/3yDnw0n9) for a 1-line C program that uses a simpler Spigot algorithm for computing the first 15000 digits\n",
    "\n",
    "A nice method, with a generator yielding the next digit each time, comes from http://stackoverflow.com/a/9005163.\n",
    "It uses only integer, so no need to do any `Decimal` trick here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def next_pi_digit(max_step):\n",
    "    q, r, t, k, m, x = 1, 0, 1, 1, 3, 3\n",
    "    for j in range(max_step):\n",
    "        if 4 * q + r - t < m * t:\n",
    "            yield m\n",
    "            # More details on Python generators can be found here http://stackoverflow.com/a/231855\n",
    "            q, r, t, k, m, x = 10*q, 10*(r-m*t), t, k, (10*(3*q+r))//t - 10*m, x\n",
    "        else:\n",
    "            q, r, t, k, m, x = q*k, (2*q+r)*x, t*x, k+1, (q*(7*k+2)+r*x)//(t*x), x+2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generator_pi(max_step):\n",
    "    big_str = ''.join(str(d) for d in next_pi_digit(max_step))\n",
    "    return Decimal(big_str[0] + '.' + big_str[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It does not use `Decimal` numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337')"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getcontext().prec = 50  # trick to improve precision\n",
    "generator_pi(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337')"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getcontext().prec = 5000  # trick to improve precision\n",
    "generator_pi(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (*hard*) [Borwein's algorithm](https://en.wikipedia.org/wiki/Borwein%27s_algorithm#Nonic_convergence)\n",
    "It has several versions, one with a cubic convergence (each new step multiplies by $3$ the number of digits), one with a nonic convergence (of order $9$) etc.\n",
    "They are not so easy to explain either.\n",
    "Please check on-line, here [en.WikiPedia.org/wiki/Borwein's_algorithm](https://en.wikipedia.org/wiki/Borwein%27s_algorithm).\n",
    "\n",
    "The cubic method is similar to the Gauss-Legendre algorithm:\n",
    "\n",
    "1. Start with $a_0 = 1/3$, and $s_0 = \\frac{\\sqrt{3}-1}{2}$,\n",
    "2. And then iterate, as much as you want, by defining $r_{k+1} = \\frac{3}{1+2(1-s_k^3)^{1/3}}$, and updating $s_{k+1} = \\frac{r_{k+1}-1}{2}$ and $a_{k+1} = r_{k+1}^2 a_k - 3^k (r_{k+1}^2 - 1)$.\n",
    "\n",
    "Then the numbers $a_k$ will converge to $\\frac{1}{\\pi}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## Examples and references\n",
    "### Links\n",
    "- [en.WikiPedia.org/wiki/Pi#Modern_quest_for_more_digits](https://en.wikipedia.org/wiki/Pi#Modern_quest_for_more_digits),\n",
    "- [www.JoyOfPi.com/pi.html](http://www.joyofpi.com/pi.html) and [www.JoyOfPi.com/pilinks.html](http://www.joyofpi.com/pilinks.html),\n",
    "- [www.EveAndersson.com/pi/digits/](http://www.eveandersson.com/pi/digits/) has great interactive tools,\n",
    "- more crazy stuff [MathWorld.Wolfram.com/PiDigits.html](http://mathworld.wolfram.com/PiDigits.html), or [MathWorld.Wolfram.com/Pi.html](http://mathworld.wolfram.com/Pi.html),\n",
    "- [one idea with Fibonacci numbers](http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibpi.html#section2),\n",
    "- [and this incredibly long list of digits](http://piworld.calico.jp/estart.html) at [PiWorld.calico.jp/estart.html](http://piworld.calico.jp/estart.html).\n",
    "- http://bellard.org/pi/ by Francois Bellard (and http://bellard.org/pi/pi_n2/pi_n2.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> That's it for today! Happy Pi Day!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "### Pie !\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/4/4b/Apple_pie.jpg/600px-Apple_pie.jpg)\n",
    "\n",
    "----"
   ]
  }
 ],
 "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.5.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "11px",
    "width": "251px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": true,
   "toc_position": {
    "height": "672px",
    "left": "0px",
    "right": "1388px",
    "top": "124px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}