Notebook

Python for Everyone!
Oregon Curriculum Network

Extended Precision with the Native Decimal Type¶

With LaTeX and Generator Functions¶

TAOCP

The Python Standard Library provides a decimal module containing the class Decimal. A Decimal object behaves according to the base 10 algorithms we learn in school.

The precision i.e. number of decimal places to which computations are carried out, is set globally, or within the scope of a context manager.

Note also that Jupyter Notebooks are able to render LaTeX when commanded to do so with the %%latex magic command. As a first example, here is an expression for the mathematical constant e.

In [1]:

%%latex
\begin{align}
e = lim_{n \to \infty} (1 + 1/n)^n
\end{align}

\begin{align} e = lim_{n \to \infty} (1 + 1/n)^n \end{align}

In [1]:

from math import e, pi
print(e)  # as a floating point number
print(pi)

2.718281828459045
3.141592653589793

Lets show setting precision to a thousand places within a scope defined by decimal.localcontext. We set precision internally to the scope. By default, precision is to 28 places. We set n to 1 followed by 102 zeros, so a very large number. The resulting computation matches a published value for e to 100 decimal places.

In [3]:

import decimal

with decimal.localcontext() as ctx:  # context manager
    ctx.prec = 1000
    n = decimal.Decimal(1e102)
    e = (1 + 1/n) ** n

In [4]:

e_1000_places = """2.7182818284590452353602874713526624977572470936999595749669
6762772407663035354759457138217852516642742746639193200305992181741359662904357
2900334295260595630738132328627943490763233829880753195251019011573834187930702
1540891499348841675092447614606680822648001684774118537423454424371075390777449
9206955170276183860626133138458300075204493382656029760673711320070932870912744
3747047230696977209310141692836819025515108657463772111252389784425056953696770
7854499699679468644549059879316368892300987931277361782154249992295763514822082
6989519366803318252886939849646510582093923982948879332036250944311730123819706
8416140397019837679320683282376464804295311802328782509819455815301756717361332
0698112509961818815930416903515988885193458072738667385894228792284998920868058
2574927961048419844436346324496848756023362482704197862320900216099023530436994
1849146314093431738143640546253152096183690888707016768396424378140592714563549
0613031072085103837505101157477041718986106873969655212671546889570350354"""
e_1000_places = e_1000_places.replace("\n",str())

In [5]:

str(e)[2:103] == e_1000_places[2:103]  # skipping "2." and going to 100 decimals

Out[5]:

True

In the context below, the value of n starts at 10 and then gets two more zeros every time around the loop.

The yield keyword is similar to return in handing back an object, however a generator function then pauses to pick up where it left off when nudged by next(), (which triggers __next__ internally).

Generator functions do not forget their internal state as they advance through next values.

Note that when a Decimal type object operates with an integer, that integer is coerced (cast) as a Decimal object.

In [6]:

with decimal.localcontext() as ctx:  # context manager
    ctx.prec = 1000
    def converge():  # generator function
        n = decimal.Decimal('10')
        while True:
            yield (1 + 1/n) ** n
            n = n * 100 # two more zeros
    f = converge()
    for _ in range(9):
        next(f) # f.__next__() <--- not quite like Python 2.x (f.next())

In [7]:

r = next(f)
r

Out[7]:

Decimal('2.718281828459045235224373380')

In [8]:

str(r)[:20] == e_1000_places[:20]

Out[8]:

True

Ramanujan Postage Stamp

The fancier LaTeX below renders a famous equation by Ramanujan, which has been shown to converge to 1/π and therefore π very quickly, relative to many other algorithms. I don't think anyone understands how some random guy could think up such a miraculaous equation.

In [9]:

%%latex
\begin{align}
\frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum^\infty_{k=0} \frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}
\end{align}

\begin{align} \frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum^\infty_{k=0} \frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}} \end{align}

Lets break the above equation into three components, a constant A and a product of Bn and Cn that accumulate as a sum. The latter two may be written as generators, with n increasing by one each time they're fed to next(). An outermost generator ties them all together and takes the reciprocal to get an approximation for π itself.

In [10]:

%%latex
\begin{align*}
\frac{1}{\pi} = \frac{2\sqrt{2}}{9801}\sum_{n = 0}^{\infty}\frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}99^{4n}} = A\sum_{n = 0}^{\infty}B_{n}C_{n},\\
A = \frac{2\sqrt{2}}{9801},\,\\ 
B_{n} = \frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}},\,\\
C_{n} = \frac{1}{99^{4n}}
\end{align*}

\begin{align*} \frac{1}{\pi} = \frac{2\sqrt{2}}{9801}\sum_{n = 0}^{\infty}\frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}99^{4n}} = A\sum_{n = 0}^{\infty}B_{n}C_{n},\\ A = \frac{2\sqrt{2}}{9801},\,\\ B_{n} = \frac{(4n)!(1103 + 26390n)}{4^{4n}(n!)^{4}},\,\\ C_{n} = \frac{1}{99^{4n}} \end{align*}

In [11]:

from math import factorial
from decimal import Decimal as D
decimal.getcontext().prec=100
A = (2 * D('2').sqrt()) / 9801
A

Out[11]:

Decimal('0.0002885855652225477091728780173879600201142070962915923014338699597981292681281720311907739076273118198')

In [12]:

def B():
    n = 0
    while True:
        numerator = factorial(4 * n) * (D(1103) + 26390 * n) 
        denominator = (4 ** (4*n))*(factorial(n))**4
        yield numerator / denominator
        n += 1

def C():
    n = 0
    while True:
        yield 1 / (D('99')**(4*n))
        n += 1
        
def Pi():
    Bn = B()
    Cn = C()
    the_sum = 0
    while True:
        the_sum += next(Bn) * next(Cn) 
        yield 1/(A * the_sum)

In [13]:

pi = Pi()

In [14]:

next(pi)
next(pi)
next(pi)

Out[14]:

Decimal('3.141592653589793238462649065702758898156677480462334781168399595644739794558841580205059234965983146')

In [15]:

pi_1000_places = """3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944
5923078164 0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172
5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461
2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273
7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652
1384146951 9415116094 3305727036 5759591953 0921861173 8193261179 3105118548 0744623799
6274956735 1885752724 8912279381 8301194912 9833673362 4406566430 8602139494 6395224737
1907021798 6094370277 0539217176 2931767523 8467481846 7669405132 0005681271 4526356082
7785771342 7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235
4201995611 2129021960 8640344181 5981362977 4771309960 5187072113 4999999837 2978049951
0597317328 1609631859 5024459455 3469083026 4252230825 3344685035 2619311881 7101000313
7838752886 5875332083 8142061717 7669147303 5982534904 2875546873 1159562863 8823537875
9375195778 1857780532 1712268066 1300192787 6611195909 2164201989"""
pi_1000_places = pi_1000_places.replace(" ","").replace("\n","")

In [16]:

r = next(pi)

In [17]:

str(r)[:20] == pi_1000_places[:20]

Out[17]:

True

Resource:

Ballad about Ramanujan by Mark Engelberg, 2005.

This I-Python Notebook is by Kirby Urner, copyleft MIT License, October 2016.