#!/usr/bin/env python # coding: utf-8 # # Pi Day Fun # March 14, 2016 # Updated since. # In[1]: from IPython.display import YouTubeVideo YouTubeVideo("HrRMnzANHHs") # Reference Pi: #
# 3.14159265358979323846264338327950288419716939937510 # 58209749445923078164062862089986280348253421170679 # 82148086513282306647093844609550582231725359408128 # 48111745028410270193852110555964462294895493038196 # 44288109756659334461284756482337867831652712019091 # 45648566923460348610454326648213393607260249141273 # 72458700660631558817488152092096282925409171536436 # 78925903600113305305488204665213841469519415116094 # 33057270365759591953092186117381932611793105118548 # 07446237996274956735188575272489122793818301194912 # 98336733624406566430860213949463952247371907021798 # 60943702770539217176293176752384674818467669405132 # 00056812714526356082778577134275778960917363717872 # 14684409012249534301465495853710507922796892589235 # 42019956112129021960864034418159813629774771309960 # 51870721134999999837297804995105973173281609631859 # 50244594553469083026425223082533446850352619311881 # 71010003137838752886587533208381420617177669147303 # 59825349042875546873115956286388235378759375195778 # 18577805321712268066130019278766111959092164201989 ## The sequence of odd fractions, as a running total, converges to pi/4, albeit slowly... # # In[2]: from fractions import Fraction from itertools import count, islice from decimal import Decimal, localcontext def convert(f): """get a Decimal from a Fraction (and multiply by 4)""" return (Decimal(f.numerator) / Decimal(f.denominator)) * Decimal(4) def pi_series(): "...converges very slowly" denoms = count(1,2) # odd numbers from 1 total = Fraction(1,next(denoms)) # 1/1 while True: yield total total -= Fraction(1, next(denoms)) # - 1/3 total += Fraction(1, next(denoms)) # + 1/5 and so on def nth(iterable, n, default=None): "Returns the nth item or a default value" return next(islice(iterable, n, None), default) with localcontext() as ctx: # <-- context manager object ctx.prec = 3000 pi = pi_series() print("{0}".format(convert(nth(pi, 1000)))[:10]) # The Youtube above describes how to use successive primes in successive terms to build a running product that converges to 2/pi. # In[3]: def Primes(): """generate successive prime numbers (trial by division)""" candidate = 1 _primes_so_far = [2] # first prime, only even prime yield _primes_so_far[-1] while True: candidate += 2 # check odds only from now on for prev in _primes_so_far: if prev**2 > candidate: yield candidate _primes_so_far.append(candidate) break if not divmod(candidate, prev)[1]: # no remainder! break # done looping p = Primes() print([next(p) for _ in range(100)]) # next 30 primes please! # In[4]: def convert(f): """get a Decimal from a Fraction (and multiply by 4)""" return (Decimal(f.denominator) / Decimal(f.numerator)) def Pi(): primes = Primes() result = Fraction(1,1) while True: p = next(primes) if divmod(p, 4)[1] == 1: term = (1 + Fraction(1,p)) else: term = (1 - Fraction(1,p)) result *= term yield result with localcontext() as ctx: # <-- context manager object ctx.prec = 300 # feel free to boost pi = Pi() print("{0}".format(convert(nth(pi, 333)))[:10]) # print("{0}".format(convert(nth(pi, 3000)))[:20]) # [Ramanujan](http://mathforum.org/kb/thread.jspa?threadID=2246748&tstart=0)