<jupyter_start><jupyter_text>A Brief Investigation of Fraction Approximations to Pi *Xiuwen Tu*_2021.10_ Pi is an important and interesting number. Because of the sustained interest, the number of decimal digits computed has been going through the rough, as shown in the next graph from Wikipedia:  ![Record approximations of pi](https://upload.wikimedia.org/wikipedia/commons/thumb/2/25/Record_pi_approximations.svg/1920px-Record_pi_approximations.svg.png) The above graph was created by Nageh in 2012. In recent years, there were at least these 3 new records: - 2019/3/14: 31.4E12 digits- 2020/1/29: 50E12 digits- 2021/8/5: 62.8E12 digits   Best fraction approximations of pi  I got interested in pi and I was learning py. So I asked myself these two questions:   What are good fraction approximations of pi?  What python features are useful in this investigation?<jupyter_code>import math
from decimal import *
from fractions import Fraction

import numpy as np
import matplotlib.pyplot as plt
# set the precision to a couple of digits more than the pi with 50 digits after the decimal  point
DefaultContext.prec = 52
setcontext(DefaultContext)
math.pi
np.pi
math.pi == np.pi
math.pi.as_integer_ratio()<jupyter_output><empty_output><jupyter_text>First, use Decimal to define pi_50 with 50 correct decimal digits<jupyter_code># First 50 digits of pi from https://en.wikipedia.org/wiki/Pi
pi_50_str = '3.14159265358979323846264338327950288419716939937510'    
len(pi_50_str)
getcontext().prec
pi_50 = Decimal(pi_50_str)
pi_50
def rel_pi_error(pi_approx): 
    return (Decimal(pi_approx) - pi_50) / pi_50
rel_pi_error(math.pi)
# this is exact, but not much insight
pi_50.as_integer_ratio()<jupyter_output><empty_output><jupyter_text>Next, use Fraction to obtain fraction approximations of pi_50<jupyter_code>Fraction(pi_50)
# this is a useful one 
Fraction(pi_50).limit_denominator()    #  the default max_denominator is 1000000)
data_list = []
for i in range(2,90): 
    max_denom = 2**i
    t = Fraction(pi_50).limit_denominator(max_denom).as_integer_ratio()
    error = (Decimal(t[0])/Decimal(t[1]) - pi_50) / pi_50
    data_list.append((max_denom, error))
    print(f'max denom: 2**{i} or {2**i},\t\t approx ratio: {t},\t\t error: {error:.1e}')<jupyter_output>max denom: 2**2 or 4,		 approx ratio: (13, 4),		 error: 3.5e-2
max denom: 2**3 or 8,		 approx ratio: (22, 7),		 error: 4.0e-4
max denom: 2**4 or 16,		 approx ratio: (22, 7),		 error: 4.0e-4
max denom: 2**5 or 32,		 approx ratio: (22, 7),		 error: 4.0e-4
max denom: 2**6 or 64,		 approx ratio: (201, 64),		 error: -3.1e-4
max denom: 2**7 or 128,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**8 or 256,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**9 or 512,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**10 or 1024,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**11 or 2048,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**12 or 4096,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**13 or 8192,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**14 or 16384,		 approx ratio: (355, 113),		 error: 8.5e-8
max denom: 2**15 or 32768,		 approx ratio: (102928, 32763),		 error: -1.1e-9
max denom: 2**16 or 65536,		 approx ratio: (104[...]<jupyter_text>A couple of observations in the above results: - The allowed max denominator exceeded 2**64 without problem. - The relative error exhausted the 50 decimal-digit precision we had in pi_50<jupyter_code># print the max_denom in scientific notation to have a little tidier output 
data_list = []
for i in range(2,90): 
    max_denom = 2**i
    t = Fraction(pi_50).limit_denominator(max_denom).as_integer_ratio()
    error = (Decimal(t[0])/Decimal(t[1]) - pi_50) / pi_50
    data_list.append((max_denom, error))
    print(f'max denom: 2**{i} or ~{2**i:.1e},\t\t approx ratio: {t},\t\t error: {error:.1e}')
len(data_list)
data_list
# drop the last 7 with questionable precisions
data_list = data_list[:-7]
data_list<jupyter_output><empty_output><jupyter_text>Now, visualize the results by plotting relative error vs max denominator<jupyter_code># separate the ones with positive errors from those with negative errors  
pos_data_list = [t for t in data_list if t[1].compare(Decimal(0))==Decimal('1')]
pos_denom_list = [t[0] for t in pos_data_list]
pos_error_list = [t[1] for t in pos_data_list]

# for negative errors, multiple by Decimal(-1) so that they show up in log plot
neg_data_list = [t for t in data_list if t[1].compare(Decimal(0))==Decimal('-1')]
neg_denom_list = [t[0] for t in neg_data_list]
neg_error_list = [t[1]*Decimal(-1) for t in neg_data_list]
fig=plt.figure(figsize=(8,6), dpi= 100, facecolor='w', edgecolor='k')
plt.scatter(pos_denom_list, pos_error_list, marker='+')
plt.scatter(neg_denom_list, neg_error_list, marker='x')
ax = plt.gca()
ax.set_ylabel('Relative error')
ax.set_yscale('log')
ax.set_xlabel('Max denominator')
ax.set_xscale('log')<jupyter_output><empty_output><jupyter_text>A few fraction approximations are good over a sizable range of max denominator Two notable ones are 22/7 and 355/113.   Let's check out the accuracy of 355/113 and 22/7 BTW, 355/113 is known as Zu's ratio, named after Zu Chongzhi (祖冲之), a Chinese mathematician in the 5th century<jupyter_code># In case we need a simple approximation to pi (but don't want to use math.pi or np.pi), we can use 355/113
pi_zu = 355/113    
rel_pi_error(pi_zu)<jupyter_output><empty_output><jupyter_text>355/113 works out to have accuracy of 85 ppb, impressive<jupyter_code>rel_pi_error(22/7)<jupyter_output><empty_output><jupyter_text>22/7 works out to have accuracy of 402 ppm, not bad  Finally, let's look into the continued fraction representation of pi Here is an image of the formula for the simple continued fraction from Wikipedia:  ![Image of Continued Fraction](https://wikimedia.org/api/rest_v1/media/math/render/svg/e294255007892e9424e6d3f49183cad6b133bbcf) Let's define a recursive function to compute the truncated first 20 terms.<jupyter_code>#  The fist 20 terms from https://mathworld.wolfram.com/PiContinuedFraction.html 
pi_cf_20 = (3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2)

def cont_fract(cf_tuple): 
    if len(cf_tuple)==1: 
        return cf_tuple[0]
    elif len(cf_tuple)==2: 
        return cf_tuple[0] + Fraction(1,cf_tuple[1])
    elif len(cf_tuple)>2: 
        x = cf_tuple[-2] + Fraction(1,cf_tuple[-1])
        new_cf_list = list(cf_tuple[:-2])
        new_cf_list.append(x)
        return cont_fract(tuple(new_cf_list))        
    else: 
        print(f'Something is wrong. Unexpected len(cf_tuple): \t{len(cf_tuple)}')
        return null 

for i in range(1,len(pi_cf_20)): 
    print(cont_fract(pi_cf_20[:i]))<jupyter_output>3
22/7
333/106
355/113
103993/33102
104348/33215
208341/66317
312689/99532
833719/265381
1146408/364913
4272943/1360120
5419351/1725033
80143857/25510582
165707065/52746197
245850922/78256779
411557987/131002976
1068966896/340262731
2549491779/811528438
6167950454/1963319607
{'fraction_non_alphanumeric': 0.11557073726159978, 'fraction_numerical': 0.1098775974950185, 'mean_word_length': 4.666935483870968, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 32, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '234556', 'license': 'MIT', 'n_tokens_mistral': 3102, 'n_tokens_neox': 2586, 'n_words': 779, 'path': 'pi_fraction_approximations_with_py.ipynb', 'repo_name': 'xiuwen-tu/pi_approx_with_py', 'size': '7026'}