<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
<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
Step1:  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).
Step2:  If we know the digits, we can directly print them
Step4:  A simple Monte-Carlo method
A simple Monte Carlo method for computing $\pi$ is to draw a circle inscribed in a square, and randomly place dots in the square.
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
Step5:  It is an interesting method, but it is just too limited for approximating digits of $\pi$.
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).
$100$ first digits of $\pi$
$\pi \simeq 3.1415926535 ~ 8979323846 ~ 2643383279 ~ 5028841971 \\ 6939937510 ~ 5820974944 ~ 5923078164 ~ 9862803482 ~ 53421170679$ when computed to the first $100$ digits.
Can you compute up to $1000$ digits? Up to $10000$ digits? Up to $100000$ digits? Up to 1 million digits?
Using mpmath
This is a simple and lazy method, using the mpmath module.
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$.
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).
Step6:  We can arbitrarily set the precision, with the constant mp.dps (digit numbers).
Step7:  Let save it for further comparison of simpler methods.
Step8:  We can solve the initial challenge easily
Step9:  And it will most probably be the quickest method presented here
Step10:  Asking for $10$ times more digits take about $100$ more of time (that's a bad news).
Step12:  The Gauss–Legendre iterative algorithm
More details can be found on this page.
The first method given here is explained in detail.
This algorithm is called the Gauss-Legendre method, and for example it was used to compute the first $206 158 430 000$ decimal digits of $\pi$ on September 18th to 20th, $1999$.
It is a very simply iterative algorithm (ie. step by step, you update the variables, as long as you want)
Step14:  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).
Using decimal.Decimal to improve precision
The limitation of this first algorithm came from using simple float numbers.
Let us now use Decimal numbers to keep as many digits after the decimal as we need.
Step15:  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).
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.
Step17:  Methods based on a convergent series
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$).
Basically, the bigger the $N$, the more digits you get (but the longer the program will run).
Some methods might be really efficient, only needing a few number of steps (a small $N$) for computing many digits.
Try to implement and compare at least two of these methods.
You can use the function sum (or math.fsum) to compute the sum, or a simple while/for loop.
All these partial sums are written up to an integer $N \geq 1$.
A Leibniz formula (easy)
Step19:  This first formula is very inefficient!
Bailey-Borwein-Plouffe series (medium)
Step20:  That's pretty impressive, in only $10$ steps!
But that algorithm is highly dependent on the precision we ask, and on the number of terms in the sum.
Step21:  It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision
Step23:  Bellard's formula (hard)
Step24:  That's pretty impressive, in only $10$ steps!
But that algorithm is again highly dependent on the precision we ask, and on the number of terms in the sum.
Step25:  It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision
Step27:  It is also slower than BBP formula.
One Ramanujan's formula (hard)
Step28:  $1595$ correct digits with $200$ terms, that's quite good!!
Step29:  Let's try to answer my initial question, using this naive implementation.
Step31:  ... It was too slow!
Chudnovsky brothers' formula (hard)
Step32:  It is very efficient, as Ramanujan's formula.
Step33:  It gets $2834$ correct numbers in $200$ steps!
It is more efficient that Ramanujan's formula, but slower to compute.
Step35:  About $2$ seconds to find correctly the first $5671$ digits? That's slow! But hey, it's Python (dynamic typing etc).
Other methods
Trigonometric methods (hard)
Some methods are based on the $\mathrm{arccot}$ or $\arctan$ functions, and use the appropriate Taylor series to approximate these functions.
The best example is Machin's formula
Step37:  Applying Machin's formula
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
Step38:  So we got the first $9995$ digits correctly... in $45$ seconds.
That will still be too slow to have the digits at position $130317$.
Step39:  The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer.
Typical implementation will be in highly efficient compiled language; like C or maybe Julia.
Many Machin-like formulas also exist, like
Step40:  It was too slow too! But at least it worked!
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$.
Conclusion
$\implies$ To the question "What are the $10$ digits of $\pi$ at index $140317..140326$", the answer is $9341076406$.
For more, see http
Step41:  It does not use Decimal numbers.
<ASSISTANT_TASK:>
Python Code:
print(" pi ~= 3.14 (two first digits).")
print(" pi ~= 22/7 = {} (two first digits).".format(22.0 / 7.0))
print(" pi ~= 355/113 = {} (six first digits).".format(355.0 / 113.0))
Explanation: Table of Contents
<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>
# Python Pi Day 2017
> 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.
Today is [Pi Day 2017](http://www.piday.org/), the day celebrating the number $\pi$.
For more details on this number, see [this Wikipedia page](https://en.wikipedia.org/wiki/Pi).
----
Let us use this occasion to showcase a few different approaches to compute the digits of the number $\pi$.
I will use the [Python](https://www.python.org/) programming language, and you are reading a [Jupyter notebook](https://www.jupyter.org/).
[![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/)
# Computing a lot of digits of $\pi$?
- **What to do ?** I will present and implement some methods that can be used to compute the first digits of the irrational number $\pi$.
- **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!
- **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.
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!
(that's the challenge I posted on Facebook)
----
## Two simple methods for finding the first digits of $\pi$
### Fraction approximations, and $\pi$ imported from the `math` module
Three approximations, using fractions, were known from a very long time (Aristote used $\frac{355}{113}$ !).
The first three approximations of pi are:
End of explanation
def mathpi():
    from math import pi
    return pi
print("First method: using math.pi gives pi ~= {:.17f} (17 digits are displayed here).".format(mathpi()))
Explanation: 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).
End of explanation
from decimal import Decimal
bigpi = Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679')
print("The first 100 digits of pi are {}.".format(bigpi))
Explanation: If we know the digits, we can directly print them:
End of explanation
from random import uniform
def montecarlo_pi(nbPoints=10000):
    Returns a probabilist estimate of pi, as a float number.
    nbInside = 0
    # we pick a certain number of points (nbPoints)
    for i in range(nbPoints):
        x = uniform(0, 1)
        y = uniform(0, 1)
        # (x, y) is now a random point in the square [0, 1] × [0, 1]
        if (x**2 + y**2) < 1:
            # This point (x, y) is inside the circle C(0, 1)
            nbInside += 1
    return 4 * float(nbInside) / floor(nbPoints)
 
print("The simple Monte-Carlo method with {} random points gave pi = {}".format(10000, montecarlo_pi(10000)))
Explanation: A simple Monte-Carlo method
A simple Monte Carlo method for computing $\pi$ is to draw a circle inscribed in a square, and randomly place dots in the square.
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:
In Python, such simulation can basically be implemented like this:
End of explanation
import mpmath
# from sympy import mpmath  # on older sympy versions
mp = mpmath.mp
Explanation: It is an interesting method, but it is just too limited for approximating digits of $\pi$.
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).
$100$ first digits of $\pi$
$\pi \simeq 3.1415926535 ~ 8979323846 ~ 2643383279 ~ 5028841971 \\ 6939937510 ~ 5820974944 ~ 5923078164 ~ 9862803482 ~ 53421170679$ when computed to the first $100$ digits.
Can you compute up to $1000$ digits? Up to $10000$ digits? Up to $100000$ digits? Up to 1 million digits?
Using mpmath
This is a simple and lazy method, using the mpmath module.
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$.
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).
End of explanation
mp.dps = 1000  # number of digits
my_pi = mp.pi  # Gives pi to a thousand places
print("A lazy method using the mpmath module:\npi is approximatly {} (with {} digits).".format(my_pi, mp.dps))
Explanation: We can arbitrarily set the precision, with the constant mp.dps (digit numbers).
End of explanation
mp.dps = 100000  # number of digits
len(str(mp.pi))
mpmath_pi = Decimal(str(mp.pi))
Explanation: Let save it for further comparison of simpler methods.
End of explanation
mp.dps = 140330
print(str(mp.pi)[2:][140317:140317+10])
Explanation: We can solve the initial challenge easily:
End of explanation
%timeit mp.dps=140330;print(str(mp.pi)[2:][140317:140317+10])
Explanation: And it will most probably be the quickest method presented here:
End of explanation
%timeit mp.dps=1403230;print(str(mp.pi)[2:][1403217:1403217+10])
Explanation: Asking for $10$ times more digits take about $100$ more of time (that's a bad news).
End of explanation
import math
def gauss_legendre_1(max_step):
    Float number implementation of the Gauss-Legendre algorithm, for max_step steps.
    a = 1.
    b = 1./math.sqrt(2)
    t = 1./4.0
    p = 1.
    for i in range(max_step):
        at = (a + b) / 2.0
        bt = math.sqrt(a*b)
        tt = t - p*(a-at)**2
        pt = 2.0 * p
        a, b, t, p = at, bt, tt, pt
    my_pi = ((a+b)**2)/(4.0*t)
    return my_pi
my_pi = gauss_legendre_1(4)
my_pi
print("pi is approximately: {:.15f} (as a float number, precision is limited).".format(my_pi))
accuracy = 100*abs(math.pi - my_pi)/math.pi
print("Accuracy % with math.pi: {:.4g}".format(accuracy))
accuracy = 100*abs(float(mpmath_pi) - my_pi)/float(mpmath_pi)
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
my_pi = gauss_legendre_1(40)
my_pi
print("pi is approximately: {:.15f} (as a float number, precision is limited).".format(my_pi))
accuracy = 100*abs(math.pi - my_pi)/math.pi
print("Accuracy % with math.pi: {:.4g}".format(accuracy))
accuracy = 100*abs(float(mpmath_pi) - my_pi)/float(mpmath_pi)
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: The Gauss–Legendre iterative algorithm
More details can be found on this page.
The first method given here is explained in detail.
This algorithm is called the Gauss-Legendre method, and for example it was used to compute the first $206 158 430 000$ decimal digits of $\pi$ on September 18th to 20th, $1999$.
It is a very simply iterative algorithm (ie. step by step, you update the variables, as long as you want):
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$.
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:
$a_{n+1} = \frac{a_n + b_n}{2}$,
$b_{n+1} = \sqrt{a_n \times b_n}$,
$t_{n+1} = t_n - p_n (a_n - a_{n+1})^2$,
$p_{n+1} = 2 p_n$.
Finally, the desired approximation of $\pi$ is: $$\pi \simeq \frac{(a_{n+1} + b_{n+1})^2}{4 t_{n+1}}.$$
Using float numbers
The first three iterations give (approximations given up to and including the first incorrect digit):
3.140 …
3.14159264 …
3.1415926535897932382 …
The algorithm has second-order convergent nature, which essentially means that the number of correct digits doubles with each step of the algorithm.
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).
This method is based on computing the limits of the arithmetic–geometric mean of some well-chosen number (see on Wikipédia for more mathematical details).
End of explanation
from decimal import Decimal, getcontext
def gauss_legendre_2(max_step):
    Decimal number implementation of the Gauss-Legendre algorithm, for max_step steps.
    # trick to improve precision
    getcontext().prec = 3 + 2**(max_step + 2)
    cst_2 = Decimal(2.0)
    cst_4 = Decimal(4.0)
    a = Decimal(1.0)
    b = Decimal(0.5).sqrt()
    t = Decimal(0.25)
    p = Decimal(1.0)
    for i in range(max_step):
        new_a = (a+b)/cst_2
        new_b = (a*b).sqrt()
        new_t = Decimal(t - p*(a - new_a)**2)
        new_p = cst_2*p
        a, b, t, p = new_a, new_b, new_t, new_p
    my_pi = Decimal(((a+b)**2)/(cst_4*t))
    return my_pi
my_pi = gauss_legendre_2(5)
print("pi is approximately: {}.".format(my_pi.to_eng_string()[:2**(5+1)]))
accuracy = 100*abs(Decimal(math.pi) - my_pi)/Decimal(math.pi)
print("Accuracy % with math.pi: {:.4g}".format(accuracy))
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: 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).
Using decimal.Decimal to improve precision
The limitation of this first algorithm came from using simple float numbers.
Let us now use Decimal numbers to keep as many digits after the decimal as we need.
End of explanation
%timeit gauss_legendre_1(8)
%timeit gauss_legendre_2(8)
Explanation: 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).
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.
End of explanation
def leibniz(max_step):
     Computing an approximation of pi with Leibniz series.
    my_pi = Decimal(0)
    for k in range(max_step):
        my_pi += Decimal((-1)**k) / Decimal(2*k+1)
    return Decimal(4) * my_pi
getcontext().prec = 20  # trick to improve precision
my_pi = leibniz(1000)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
getcontext().prec = 20  # trick to improve precision
my_pi = leibniz(10000)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: Methods based on a convergent series
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$).
Basically, the bigger the $N$, the more digits you get (but the longer the program will run).
Some methods might be really efficient, only needing a few number of steps (a small $N$) for computing many digits.
Try to implement and compare at least two of these methods.
You can use the function sum (or math.fsum) to compute the sum, or a simple while/for loop.
All these partial sums are written up to an integer $N \geq 1$.
A Leibniz formula (easy):
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.
$$\pi \simeq 4\sum_{n=0}^{N} \cfrac {(-1)^n}{2n+1}. $$
End of explanation
def bbp(max_step):
     Computing an approximation of pi with Bailey-Borwein-Plouffe series.
    my_pi = Decimal(0)
    for k in range(max_step):
        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)))
    return my_pi
getcontext().prec = 20  # trick to improve precision
my_pi = bbp(10)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: This first formula is very inefficient!
Bailey-Borwein-Plouffe series (medium):
It also has a number of digits linear in the number $N$ of terms in the sum.
$$\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). $$
End of explanation
getcontext().prec = 200  # trick to improve precision
my_pi = bbp(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
getcontext().prec = 500  # trick to improve precision
my_pi = bbp(500)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: That's pretty impressive, in only $10$ steps!
But that algorithm is highly dependent on the precision we ask, and on the number of terms in the sum.
End of explanation
getcontext().prec = 10 + 1000  # trick to improve precision
%timeit bbp(1000)
getcontext().prec = 10 + 2000  # trick to improve precision
%timeit bbp(2000)
Explanation: It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision:
End of explanation
def bellard(max_step):
     Computing an approximation of pi with Bellard series.
    my_pi = Decimal(0)
    for k in range(max_step):
        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))
    return my_pi * Decimal(1.0/(2**6))
getcontext().prec = 40  # trick to improve precision
my_pi = bellard(10)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: Bellard's formula (hard):
It is a more efficient formula.
$$\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). $$
End of explanation
getcontext().prec = 800  # trick to improve precision
my_pi = bellard(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: That's pretty impressive, in only $10$ steps!
But that algorithm is again highly dependent on the precision we ask, and on the number of terms in the sum.
End of explanation
getcontext().prec = 10 + 1000  # trick to improve precision
%timeit bellard(1000)
getcontext().prec = 10 + 2000  # trick to improve precision
%timeit bellard(2000)
Explanation: It is, of course, slower than the optimized algorithm from mpmath.
And it does not scale well with the precision:
End of explanation
from math import factorial
def ramanujan(max_step):
     Computing an approximation of pi with a Ramanujan's formula.
    my_pi = Decimal(0)
    d_1103 = Decimal(1103)
    d_26390 = Decimal(26390)
    d_396 = Decimal(396)
    for k in range(max_step):
        my_pi += ((Decimal(factorial(4 * k))) * (d_1103 + d_26390 * Decimal(k))) / ( (Decimal(factorial(k)))**4 * (d_396**(4*k)))
    my_pi = my_pi * 2 * Decimal(2).sqrt() / Decimal(9801)
    my_pi = my_pi**(-1)
    return my_pi
getcontext().prec = 40  # trick to improve precision
my_pi = ramanujan(4)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
getcontext().prec = 400  # trick to improve precision
my_pi = ramanujan(40)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
getcontext().prec = 2000  # trick to improve precision
my_pi = ramanujan(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: It is also slower than BBP formula.
One Ramanujan's formula (hard):
It is efficient but harder to compute easily with high precision, due to the factorial.
But hopefully, the function math.factorial works with Python integers, of arbitrary size, so this won't be an issue.
$$\frac{1}{\pi} \simeq \frac{2\sqrt{2}}{9801} \sum_{n=0}^N \frac{(4n)!(1103+26390n)}{(n!)^4 396^{4n}}. $$
Remark: This method and the next one compute approximation of $\frac{1}{\pi}$, not $\pi$.
This method has a quadratic precision (the number of correct digits is of the order $\mathcal{O}(N^2)$.
End of explanation
getcontext().prec = 10 + 2000  # trick to improve precision
%timeit ramanujan(200)
getcontext().prec = 10 + 5000  # trick to improve precision
%timeit ramanujan(400)
Explanation: $1595$ correct digits with $200$ terms, that's quite good!!
End of explanation
%%time
getcontext().prec = 140350  # trick to improve precision
i = 140317
my_pi = ramanujan(10000)
print(str(my_pi)[2:][i:i+10])
mp.dps=140330
print(str(mp.pi)[2:][i:i+10])
Explanation: Let's try to answer my initial question, using this naive implementation.
End of explanation
from math import factorial
def chudnovsky(max_step):
     Computing an approximation of pi with Bellard series.
    my_pi = Decimal(0)
    for k in range(max_step):
        my_pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
    my_pi = my_pi * Decimal(10005).sqrt()/4270934400
    my_pi = my_pi**(-1)
    return my_pi
Explanation: ... It was too slow!
Chudnovsky brothers' formula (hard):
$$\frac{1}{\pi} \simeq 12 \sum_{n=0}^N \frac {(-1)^{n}(6n)!(545140134n+13591409)}{(3n)!(n!)^{3}640320^{{3n+3/2}}}. $$
In 2015, the best method is still the Chudnovsky brothers's formula.
Be careful when you use these formulae, check carefully the constants you wrote (545140134 will work well, 545140135 might not work as well!).
This formula is used as the logo of the French agrégation of Mathematics website agreg.org :
End of explanation
getcontext().prec = 3000  # trick to improve precision
my_pi = chudnovsky(200)
my_pi
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: It is very efficient, as Ramanujan's formula.
End of explanation
getcontext().prec = 6000  # trick to improve precision
my_pi = chudnovsky(400)
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
getcontext().prec = 3000  # trick to improve precision
%timeit chudnovsky(200)
getcontext().prec = 6000  # trick to improve precision
%timeit chudnovsky(400)
Explanation: It gets $2834$ correct numbers in $200$ steps!
It is more efficient that Ramanujan's formula, but slower to compute.
End of explanation
def arccot(x, unity):
    Compute arccot(x) with a certain level of precision.
    x = Decimal(x)
    unity = Decimal(unity)
    mysum = xpower = unity / x
    n = 3
    sign = -1
    while True:
        xpower = xpower / (x*x)
        term = xpower / n
        if not term:
            break
        mysum += sign * term
        sign = -sign  # we alternate the sign
        n += 2
    return mysum
Explanation: About $2$ seconds to find correctly the first $5671$ digits? That's slow! But hey, it's Python (dynamic typing etc).
Other methods
Trigonometric methods (hard)
Some methods are based on the $\mathrm{arccot}$ or $\arctan$ functions, and use the appropriate Taylor series to approximate these functions.
The best example is Machin's formula:
$$\pi = 16 \;\mathrm{arccot}(5) - 4 \;\mathrm{arccot}(239).$$
And we use the Taylor series to approximate $\mathrm{arccot}(x)$:
$$\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}} .$$
This method is also explained here with some details.
In order to obtain $n$ digits, we will use fixed-point arithmetic to compute $\pi \times 10^n$ as a Python long integer.
High-precision arccot computation
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.
We then repeatedly divide by $x^2$ and a counter value that runs over $3$, $5$, $7$ etc, to obtain each next term.
The summation is stopped at the first zero term, which in this fixed-point representation corresponds to a real value less than $10^{-n}$:
python
def arccot(x, unity):
    xpower = unity / x
    sum = xpower
    n = 3
    sign = -1
    while True:
        xpower = xpower / (x*x)
        term = xpower / n
        if term == 0:
            break  # we are done
        sum += sign * term
        sign = -sign
        n += 2
    return sum
Adapting it to use Decimal numbers is easy:
End of explanation
def machin(digits):
    Compute pi with Machin's formula, with precision at least digits.
    unity = 10**(digits + 10)
    my_pi = Decimal(4) * (Decimal(4)*arccot(5, unity) - arccot(239, unity))
    return my_pi / Decimal(unity)
getcontext().prec = 10000  # trick to improve precision
my_pi = machin(100)
accuracy = 100*abs(mpmath_pi - my_pi)/mpmath_pi
print("Accuracy % with mpmath_pi: {:.4g}".format(accuracy))
Explanation: Applying Machin's formula
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:
$$\pi = 16 \;\mathrm{arccot}(5) - 4 \;\mathrm{arccot}(239).$$
python
def machin(digits):
    unity = 10**(digits + 10)
    pi = 4 * (4*arccot(5, unity) - arccot(239, unity))
    return pi / unity
To avoid rounding errors in the result, we use 10 guard digits internally during the calculation.
We may now reproduce the historical result obtained by Machin in 1706.
End of explanation
getcontext().prec = 5000  # trick to improve precision
%timeit machin(50)
getcontext().prec = 10000  # trick to improve precision
%timeit machin(100)
Explanation: So we got the first $9995$ digits correctly... in $45$ seconds.
That will still be too slow to have the digits at position $130317$.
End of explanation
%%time
i = 14031
getcontext().prec = i + 20  # trick to improve precision
mp.dps = i + 20
print(str(mp.pi)[2:][i:i+10])
my_pi = machin(11)
print(str(my_pi)[2:][i:i+10])
%%time
i = 140317
getcontext().prec = i + 20  # trick to improve precision
mp.dps = i + 20
print(str(mp.pi)[2:][i:i+10])
my_pi = machin(50)
print(str(my_pi)[2:][i:i+10])
Explanation: The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer.
Typical implementation will be in highly efficient compiled language; like C or maybe Julia.
Many Machin-like formulas also exist, like:
$$\pi = 4\arctan\left(\frac{1}{2}\right) + 4 \arctan\left(\frac{1}{3}\right).$$
Trying to solve my question!
The important parameter to tune is not the precision given to machin() but the Decimal.prec precision.
End of explanation
def next_pi_digit(max_step):
    q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
    for j in range(max_step):
        if 4 * q + r - t < m * t:
            yield m
            # More details on Python generators can be found here http://stackoverflow.com/a/231855
            q, r, t, k, m, x = 10*q, 10*(r-m*t), t, k, (10*(3*q+r))//t - 10*m, x
        else:
            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
def generator_pi(max_step):
    big_str = ''.join(str(d) for d in next_pi_digit(max_step))
    return Decimal(big_str[0] + '.' + big_str[1:])
Explanation: It was too slow too! But at least it worked!
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$.
Conclusion
$\implies$ To the question "What are the $10$ digits of $\pi$ at index $140317..140326$", the answer is $9341076406$.
For more, see http://fredrikj.net/blog/2011/03/100-mpmath-one-liners-for-pi/.
(hard) Unbounded Spigot Algorithm
This algorithm is quite efficient, but not easy to explain. Go check on-line if you want.
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
A nice method, with a generator yielding the next digit each time, comes from http://stackoverflow.com/a/9005163.
It uses only integer, so no need to do any Decimal trick here.
End of explanation
getcontext().prec = 50  # trick to improve precision
generator_pi(1000)
getcontext().prec = 5000  # trick to improve precision
generator_pi(1000)
Explanation: It does not use Decimal numbers.
End of explanation
<END_TASK>