from sympy import *
import sympy as sp
x = sp.symbols('x')
def fu(n,x):
if n==0: return x
else: return (fu(n-1,x)+fv(n-1,x))/2
def fv(n,x):
if n==0: return 1
else: return sqrt(fu(n-1,x)*fv(n-1,x))
def fy(n,y):
return (fu(n,x)/fv(n,x)).subs(x,y)
def fz(n,y):
return (diff(fv(n,x),x)/diff(fu(n,x),x)).subs(x,y)
def fpi(n):
if n==0: return (2+sqrt(2))
else: return fpi(n-1) * (1 + fy(n,1/sqrt(2))) / (1 + fz(n,1/sqrt(2)))
for n in range(10):
p = 1000
print('n='+ str(n)
+ ': erreur 10**('
+ str(N(ln(N(fpi(n)-pi,p))/ln(10),5))
+ ')')
en faisant les calculs approchés avec $1000$ chiffres significatifs, l'exécution de ce code donne $\pi$ avec 700 décimales au bout de seulement 8 étapes (la précision double à chaque étape, en fait):
n=0: erreur 10**(-0.56444)
n=1: erreur 10**(-2.9939)
n=2: erreur 10**(-8.1322)
n=3: erreur 10**(-18.737)
n=4: erreur 10**(-40.262)
n=5: erreur 10**(-83.619)
n=6: erreur 10**(-170.64)
n=7: erreur 10**(-344.98)
n=8: erreur 10**(-693.95)
