# 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>?</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></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> 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> 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></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/).

[](http://jupyter.org/)[](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:


```python
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))
```

     pi ~= 3.14 (two first digits).
     pi ~= 22/7 = 3.142857142857143 (two first digits).
     pi ~= 355/113 = 3.1415929203539825 (six first digits).


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).


```python
def mathpi():
    from math import pi
    return pi

print("First method: using math.pi gives pi ~= {:.17f} (17 digits are displayed here).".format(mathpi()))
```

    First method: using math.pi gives pi ~= 3.14159265358979312 (17 digits are displayed here).


If we know the digits, we can directly print them:


```python
from decimal import Decimal
bigpi = Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679')
print("The first 100 digits of pi are {}.".format(bigpi))
```

    The first 100 digits of pi are 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679.


### A simple [Monte-Carlo method](https://en.wikipedia.org/wiki/Pi#Monte_Carlo_methods)
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:


```python
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)
 
```


```python
print("The simple Monte-Carlo method with {} random points gave pi = {}".format(10000, montecarlo_pi(10000)))
```

    The simple Monte-Carlo method with 10000 random points gave pi = 3.138


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`](http://mpmath.org/)
This is a simple and lazy method, using the [`mpmath`](http://mpmath.org/) 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).


```python
import mpmath
# from sympy import mpmath  # on older sympy versions
mp = mpmath.mp
```

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).


```python
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))
```

    A lazy method using the mpmath module:
    pi is approximatly 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198 (with 1000 digits).


Let save it for further comparison of simpler methods.


```python
mp.dps = 100000  # number of digits
len(str(mp.pi))
mpmath_pi = Decimal(str(mp.pi))
```




    100001



We can solve the initial challenge easily:


```python
mp.dps = 140330
print(str(mp.pi)[2:][140317:140317+10])
```

    9341076406


And it will most probably be the quickest method presented here:


```python
%timeit mp.dps=140330;print(str(mp.pi)[2:][140317:140317+10])
```

    9341076406
    9341076406
    9341076406
    9341076406
    1 loop, best of 3: 236 ms per loop


Asking for $10$ times more digits take about $100$ more of time (that's a bad news).


```python
%timeit mp.dps=1403230;print(str(mp.pi)[2:][1403217:1403217+10])
```

    7966916520
    7966916520
    7966916520
    7966916520
    The slowest run took 6.34 times longer than the fastest. This could mean that an intermediate result is being cached.
    1 loop, best of 3: 20.4 s per loop


## The Gauss–Legendre iterative algorithm
> More details can be found on [this page](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm).

The first method given here is explained in detail.
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$.

It is a very simply **iterative algorithm** (ie. step by step, you update the variables, as long as you want):

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$.

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:
   - $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$.

3. 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](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)).


```python
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
```


```python
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))
```




    3.141592653589794



    pi is approximately: 3.141592653589794 (as a float number, precision is limited).
    Accuracy % with math.pi: 2.827e-14
    Accuracy % with mpmath_pi: 2.827e-14



```python
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))
```




    3.141592653589794



    pi is approximately: 3.141592653589794 (as a float number, precision is limited).
    Accuracy % with math.pi: 2.827e-14
    Accuracy % with mpmath_pi: 2.827e-14


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`](https://docs.python.org/3/library/decimal.html) numbers to keep as many digits after the decimal as we need.


```python
from decimal import Decimal, getcontext
```


```python
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
```


```python
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))
```

    pi is approximately: 3.14159265358979323846264338327950288419716939937510582097494459.
    Accuracy % with math.pi: 3.898e-15
    Accuracy % with mpmath_pi: 7.659e-83


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.


```python
%timeit gauss_legendre_1(8)
%timeit gauss_legendre_2(8)
```

    The slowest run took 6.47 times longer than the fastest. This could mean that an intermediate result is being cached.
    100000 loops, best of 3: 4.06 µs per loop
    100 loops, best of 3: 3.9 ms per loop


## 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](https://en.wikipedia.org/wiki/Leibniz_formula_for_pi) (*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}. $$


```python
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
```


```python
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))
```




    Decimal('3.1405926538397929257')



    Accuracy % with mpmath_pi: 0.03183



```python
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))
```




    Decimal('3.1414926535900432389')



    Accuracy % with mpmath_pi: 0.003183


This first formula is very inefficient!

### [Bailey-Borwein-Plouffe series](https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula) (*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). $$


```python
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
```


```python
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))
```




    Decimal('3.1415926535897911465')



    Accuracy % with mpmath_pi: 6.659e-14


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.


```python
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))
```




    Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303817')



    Accuracy % with mpmath_pi: 8.417e-198



```python
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))
```




    Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119494')



    Accuracy % with mpmath_pi: 8.600e-498


It is, of course, slower than the optimized algorithm from `mpmath`.
And it does not scale well with the precision:


```python
getcontext().prec = 10 + 1000  # trick to improve precision
%timeit bbp(1000)

getcontext().prec = 10 + 2000  # trick to improve precision
%timeit bbp(2000)
```

    10 loops, best of 3: 95 ms per loop
    1 loop, best of 3: 661 ms per loop


### [Bellard's formula](https://en.wikipedia.org/wiki/Bellard%27s_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). $$


```python
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))
```


```python
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))
```




    Decimal('3.141592653589793238462643383279490036578')



    Accuracy % with mpmath_pi: 4.090e-31


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.


```python
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))
```




    Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005611526979981607242018469874601593584399987639992795109821175268982496078980620087198769557504210588066903133631260918929384314932145237223424648852944318026933908108945450734428224399103245946394')



    Accuracy % with mpmath_pi: 2.220e-604


It is, of course, slower than the optimized algorithm from `mpmath`.
And it does not scale well with the precision:


```python
getcontext().prec = 10 + 1000  # trick to improve precision
%timeit bellard(1000)

getcontext().prec = 10 + 2000  # trick to improve precision
%timeit bellard(2000)
```

    1 loop, best of 3: 378 ms per loop
    1 loop, best of 3: 2.91 s per loop


It is also slower than BBP formula.

### One [Ramanujan's formula](https://en.wikipedia.org/wiki/Approximations_of_%CF%80#Efficient_methods) (*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)$.


```python
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
```


```python
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))
```




    Decimal('3.141592653589793238462643383279555273157')



    Accuracy % with mpmath_pi: 1.668e-30



```python
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))
```




    Decimal('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588249728368707403112841783926380916831222457885852997711270901024127742895829505936')



    Accuracy % with mpmath_pi: 2.382e-318



```python
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))
```




    Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029553211653449872027559602364806654991198818347977535663698074265425278625518184175746728909777727940092593679577412880114072520662452344498309284389026310583057432459663517241665671357564259409181627161036619398380266804323123256475817186916514676524019320976885611753635255011894984447618075599914209997043044425803288826839697424724950586712681051712136647565220069328838660790206792014356556715545181106432983263815336472301873472006826126207435983354628153181952126467499415321188698738323941189429')



    Accuracy % with mpmath_pi: 6.658e-1596


$1595$ correct digits with $200$ terms, that's quite good!!


```python
getcontext().prec = 10 + 2000  # trick to improve precision
%timeit ramanujan(200)

getcontext().prec = 10 + 5000  # trick to improve precision
%timeit ramanujan(400)
```

    10 loops, best of 3: 50.1 ms per loop
    1 loop, best of 3: 478 ms per loop


Let's try to answer my initial question, using this naive implementation.


```python
%%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])
```

... It was too slow!

### [Chudnovsky brothers' formula](https://en.wikipedia.org/wiki/Chudnovsky_algorithm) (*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](https://en.wikipedia.org/wiki/Agr%C3%A9gation) [website `agreg.org`](http://agreg.org/) :
> 


```python
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
```

It is very efficient, as Ramanujan's formula.


```python
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))
```




    Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009946576407895126946839835259570982582262052248940772671947826848260147699090264013639443745530506820349625245174939965143142980919065925093722169646151570985838741059788595977297549893016175392846813826868386894277415599185592524595395943104997252468084598727364469584865383673622262609912460805124388439045124413654976278079771569143599770012961608944169486855584840635342207222582848864815845602850601684273945226746767889525213852254995466672782398645659611635488623057745649803559363456817432411251507606947945109659609402522887971089314566913686722874894056010150330861792868092087476091782493858900971490967598526136554978189312978482168299894872265880485756401427047755513237964145152374623436454285844479526586782105114135473573952311342716610213596953623144295248493718711014576540359027993440374200731057853906219838744780847852710405603214833731504646600861733152271806431614075966417958570002328166126409832551776784733662072536424251339850953855478119486626262729703543651855423452347092')



    Accuracy % with mpmath_pi: 1.191e-2835


It gets $2834$ correct numbers in $200$ steps!
It is more efficient that Ramanujan's formula, but slower to compute.


```python
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))
```

    Accuracy % with mpmath_pi: 3.947e-5672



```python
getcontext().prec = 3000  # trick to improve precision
%timeit chudnovsky(200)

getcontext().prec = 6000  # trick to improve precision
%timeit chudnovsky(400)
```

    1 loop, best of 3: 216 ms per loop
    1 loop, best of 3: 1.93 s per loop


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](http://en.literateprograms.org/Pi_with_Machin%27s_formula_%28Python%29):
$$\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](http://en.literateprograms.org/Pi_with_Machin%27s_formula_%28Python%29#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:


```python
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
```

#### 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](https://en.wikipedia.org/wiki/John_Machin).


```python
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)
```


```python
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))
```

    Accuracy % with mpmath_pi: 1.501e-9996


So we got the first $9995$ digits correctly... in $45$ seconds.
That will still be too slow to have the digits at position $130317$.


```python
getcontext().prec = 5000  # trick to improve precision
%timeit machin(50)

getcontext().prec = 10000  # trick to improve precision
%timeit machin(100)
```

    1 loop, best of 3: 22.7 s per loop
    1 loop, best of 3: 44.5 s per loop


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](https://en.wikipedia.org/wiki/Machin-like_formula) 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.


```python
%%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])
```

    7126218283
    7126218283
    CPU times: user 1min 3s, sys: 12 ms, total: 1min 3s
    Wall time: 1min 3s



```python
%%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])
```

    9341076406
    9341076406
    CPU times: user 10min 40s, sys: 224 ms, total: 10min 40s
    Wall time: 10min 41s


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](http://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf)
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.


```python
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
```


```python
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:])
```

It does not use `Decimal` numbers.


```python
getcontext().prec = 50  # trick to improve precision
generator_pi(1000)
```




    Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337')




```python
getcontext().prec = 5000  # trick to improve precision
generator_pi(1000)
```




    Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337')



### (*hard*) [Borwein's algorithm](https://en.wikipedia.org/wiki/Borwein%27s_algorithm#Nonic_convergence)
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.
They are not so easy to explain either.
Please check on-line, here [en.WikiPedia.org/wiki/Borwein's_algorithm](https://en.wikipedia.org/wiki/Borwein%27s_algorithm).

The cubic method is similar to the Gauss-Legendre algorithm:

1. Start with $a_0 = 1/3$, and $s_0 = \frac{\sqrt{3}-1}{2}$,
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)$.

Then the numbers $a_k$ will converge to $\frac{1}{\pi}$.

----
## Examples and references
### Links
- [en.WikiPedia.org/wiki/Pi#Modern_quest_for_more_digits](https://en.wikipedia.org/wiki/Pi#Modern_quest_for_more_digits),
- [www.JoyOfPi.com/pi.html](http://www.joyofpi.com/pi.html) and [www.JoyOfPi.com/pilinks.html](http://www.joyofpi.com/pilinks.html),
- [www.EveAndersson.com/pi/digits/](http://www.eveandersson.com/pi/digits/) has great interactive tools,
- 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),
- [one idea with Fibonacci numbers](http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibpi.html#section2),
- [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).
- http://bellard.org/pi/ by Francois Bellard (and http://bellard.org/pi/pi_n2/pi_n2.html)

> That's it for today! Happy Pi Day!

----
### Pie !


----

b82de72c9d52975789d9b7c32c1590d9f3a381e8
ipynb
Jupyter Notebook
Py_Pi_Day_2017.ipynb
IEWbgfnYDwHRoRRSKtkdyMDUzgdwuBYgDKtDJWd/narnt
0eda13a7b8663e218b4fe2e06a974b99db9ff166
MIT
2016-06-25T09:30:00.000Z
2022-03-24T21:02:49.000Z
Py_Pi_Day_2017.ipynb
operade/notebooks
56f97e33e81b5e86905961b09184a41b7616fa90
MIT
2016-06-26T12:21:30.000Z
2021-04-06T09:19:49.000Z
Py_Pi_Day_2017.ipynb
operade/notebooks
56f97e33e81b5e86905961b09184a41b7616fa90
MIT
2017-05-13T23:54:56.000Z
2021-07-17T15:34:24.000Z
Qwen/Qwen-72B
1. YES
2. YES