Le module decimal de Python permet de faire des calculs décimaux de grande précision, avec des nombres comportant un nombre de chiffres significatifs quelconque, aussi grand qu'on le veut.
Mais ce module ne comporte pas les fonctions telles que log, sin, cos, atan, etc… Avec Python v2.5, nous avons cependant la racine carrée (sqrt), l'exponentielle (exp) ainsi que la puissance (power).
L'objectif ici est de définir les fonctions qui manquent.
Nous utiliserons les séries entières classiques.
Apparemment, ce ne sont pas ces séries qui sont utilisées dans les calculatrices, mais l'algorithme CORDIC (http://fr.wikipedia.org/wiki/CORDIC). Cependant, depuis 1959, les ordinateurs ont fait beaucoup de progrès, et les calculs basés sur les séries sont devenus très rapides.
Ces séries comportent une somme de termes de taille décroissante, dont le nombre dépendra du nombre de chiffres significatifs voulu. Il existe de nombreux sites qui parlent de cela sur internet, et en particulier:
Nous calculerons aussi les nombres pi et e avec un nombre quelconque de décimales.
En ce qui concerne pi, il y a aussi de nombreux sites (passionnés) qui en parlent. Par exemple, avec le site suivant, vous avez accès aux premières 450 000 décimales de pi:
Et en voilà 1 million:
J'ai quelques doutes sur l'utilité réelle de calculer des nombres aussi grands, mais c'est très amusant…
Dernier point général, les codes qui suivent supposent:
Nous utiliserons la série suivante:
<m>sin(x) = sum{n=0}{infty}1_n_x_2n_1} / {(2n+1)!}}}</m>
Ce qui nous donne:
<m>sin(x) = 1 - {x^3/{3!}} + {x^5/{5!}} - … +1_n_x_2n_1}/{(2n+1)!}}</m>
Bien entendu, nous n'essaierons pas de calculer les factorielles ni les puissances! Nous utiliserons une formule de récurrence, permettant de calculer chaque nouveau terme à partir du précédent:
Voilà le code:
def sindec(x): """Calcul du sinus de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) xc = x*x t = x s1 = t n = 0 while True: n += 1 t *= -xc/(2*n*(2*n+1)) s2 = s1 + t if s2==s1: break s1 = s2 return s2
Validité: avec x quelconque en radians. Bien entendu, avec x dans l'intervalle ]pi,2pi[ ± 2kpi, la valeur renvoyée sera négative.
Exemples d'utilisation avec 50 chiffres significatifs:
getcontext().prec = 50 x= 0.0842123359644 radians print sindec(x) 0.084112836235638713981716879773177446578104393790879 print math.sin(x) 0.0841128362356 x= 6.02435126247 radians print sindec(x) -0.25595360989556230873783592652804681527158784625035 print math.sin(x) -0.255953609896
Est que les chiffres au delà des 12 premiers sont justes? ⇒ voir plus loin les tests sur les cosinus.
Est-ce que c'est rapide? Oh, oui!
Nous utiliserons la série suivante:
<m>cos(x) = sum{n=0}{infty}1_n_x_2n} / {(2n)!}}}</m>
Ce qui nous donne:
<m>cos(x) = 1 - {x^2/{2!}} + {x^4/{4!}} - … +1_n_x_2n}/{(2n)!}}</m>
Bien entendu, nous n'essaierons pas de calculer les factorielles ni les puissances! Nous utiliserons une formule de récurrence, permettant de calculer chaque nouveau terme à partir du précédent:
Voilà le code:
def cosdec(x): """Calcul du cosinus de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) xc = x*x s1 = 1 t = s1 n = 0 while True: n += 1 t *= -xc/((2*n-1)*2*n) s2 = s1 + t if s2 == s1: break s1 = s2 return s2
Validité: avec x quelconque en radians. Bien entendu, avec x dans l'intervalle ]pi/2,3pi/2[ ± 2kpi, la valeur renvoyée sera négative.
Exemples d'utilisation avec 50 chiffres significatifs:
getcontext().prec = 50 x= 0.412336234282 # radians print cosdec(x) 0.91618707606927757352331972478557308925585912162648 print math.cos(x) 0.916187076069 x= 6.15281253742 # radians print cosdec(x) 0.99151350113535082488508488861516895705624721837360 print math.cos(x) 0.991513501135
Est que les chiffres au delà des 12 premiers sont justes?
Je n'ai pas réussi à trouver sur internet des valeurs assez grandes de sinus et cosinus pour vérifier mes algorithmes. Mais il y a au moins quelque chose de sûr: “sinus carré + cosinus carré” est censé être égal à 1! Et comme mes sinus et cosinus ont été obtenus par des séries différentes, ça vaut le coup de vérifier:
getcontext().prec = 50 x = 6.15281253742 # radians print sindec(x) -0.13000375789306499309100040521235111350540932255885 print cosdec(x) 0.99151350113535082488508488861516895705624721837360 print sindec(x)**2 + cosdec(x)**2 1.0000000000000000000000000000000000000000000000015 x = 5.03247213801 # radians print sindec(x) -0.94920925619221626023725482541844871626870748221895 print cosdec(x) 0.31464549569161093779987720367689854570958208841894 print sindec(x)**2 + cosdec(x)**2 1.0000000000000000000000000000000000000000000000006 x = 5.91654196161 # radians print sindec(x) -0.35848389994454140436127603106340497042737673508754 print cosdec(x) 0.93353590904718390875571747382030888354771755175655 print sindec(x)**2 + cosdec(x)**2 1.0000000000000000000000000000000000000000000000004
Eh! ça ne marche pas trop mal, non?
Est-ce que c'est rapide? oui! Ce sont les même temps que pour le sinus: 3/1000 de seconde en moyenne à chaque appel de fonction pour 28 chiffres significatifs, et env. 6/1000 de seconde pour 56 chiffres significatifs.
Le calcul des tangentes avec les séries est un peu plus complexe que pour les fonctions précédentes. Dans un premier temps, je me contenterai de dire que la tangente est égal au rapport entre le sinus et le cosinus. Je reprendrai le calcul de la tangente lorsque j'aurai trouvé un algorithme rapide pour calculer les nombres de Bernoulli.
Donc:
def tandec(x): """Calcul de la tangente de x""" return sindec(x)/cosdec(x)
Bien entendu, pour x=pi/2 ± kpi, la tangente tendra vers l'infini et donnera donc une exception à gérer.
Nous utiliserons la série suivante:
<m>arcsin(x) = sum{n=0}{infty}a_nx_2n_1} / {2n+1}}}</m>
avec pour <m>a_n</m>:
Ce qui nous donne:
<m>arcsin(x) = x + x_3_2.3 + 1.3.x_5_2.4.5 + 1.3.5.x_7_2.4.6.7 +…+1.3.5..._2n-1_x_2n_1/{2.4.6…(2n)(2n+1)}}</m>
Comme pour les fonctions précédentes, nous utiliserons une formule de récurrence, permettant de calculer chaque nouveau terme à partir du précédent:
Voilà le code. Au départ, la 1ère fonction “_asindec(x)” aurait du suffire, mais elle a des problèmes de convergence lorsqu'on s'approche de x=1. Alors, on restreint ce calcul pour des valeurs de abs(x) inférieures à 0.70710678118654752 (en fait, la valeur exacte importe peu), c'est à dire pour un angle d'environ pi/4. Quand on sort de cette plage, on calcule en fait l'arc cosinus et on corrige ensuite l'angle obtenu avec pi/2. Bien entendu, le 'pi' utilisé doit avoir la même précision que les autres nombres: voir calcul de pi plus loin dans cette page.
Notez qu'on utilise la fonction racine sqrtdec(x) définie par ailleurs sur cette page.
def _asindec(x): """utilitaire pour le calcul de l'arc sinus de x""" xc = x*x k = x s1 = k n = 0 while True: n += 1 k *= xc*(2*n-1)*(2*n-1)/(2*n*(2*n+1)) s2 = s1 + k if s2==s1: break s1 = s2 return s2 def asindec(x): """Calcul de l'arc sinus de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) if abs(x)<Decimal("0.70710678118654752"): return _asindec(x) if x>=0: return pi/2-_asindec(sqrtdec(1-x*x)) return -pi/2+_asindec(sqrtdec(1-x*x))
Validité: pour tout x dans l'intervalle [-1,+1] puisque x est censé être un sinus.
L'arc sinus renvoyé est en radian.
La meilleure façon de vérifier la validité des résultats, c'est de partir d'un angle, de calculer son sinus et de vérifier qu'on retrouve bien l'angle initial:
getcontext().prec = 50 x = 0.863846777147 # radians print sindec(x) 0.76034673109341433401639027976672470571419234482224 print asindec(sindec(x)) 0.86384677714700000000000000000000000000000000000000 x = 0.60455808921 # radians print sindec(x) 0.56839854818359705731026717641031205281762300674039 print asindec(sindec(x)) 0.60455808921000000000000000000000000000000000000000 x = 0.430301700124 # radians print sindec(x) 0.41714501853203110194707214854588640791066844824553 print asindec(sindec(x)) 0.43030170012400000000000000000000000000000000000000
On voit bien qu'on retrouve l'angle initial avec un degré élevé de précision.
def acosdec(x): """Calcul de l'arc sinus de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) return pi/2-asindec(x)
def atandec(x): """Calcul de l'arc tangente de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) return asindec(x/sqrtdec(1+x*x))
Nous utiliserons la série suivante:
<m>log(x) = 2 sum{n=0}{infty}{ {(x-1)^{2n+1}} / {(2n+1)(x+1)^{2n+1}} }</m>
Ce qui nous donne:
<m>log(x) = 2 ( x-1_x_1 + x-1_3_3_x_1_3 + x-1_5_5_x_1_5 + x-1_7_7_x_1_7 + …)</m>
Bien entendu, nous n'essaierons pas de calculer les factorielles ni les puissances! Nous utiliserons une formule de récurrence, permettant de calculer chaque nouveau terme à partir du précédent:
Voilà le code:
def logdec(x): """Calcul du logarithme népérien de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) t = (x-1)/(x+1) tc = t*t s1 = t n = 0 while True: n += 1 t *= tc*(2*n-1)/(2*n+1) s2 = s1 + t if s2 == s1: break s1 = s2 return 2*s2
Validité: pour x>0
Quelques vérifications:
getcontext().prec = 50 x = 1 print logdec(x) 0 print repr(math.log(x)) 0.0 x = 2 print logdec(x) 0.69314718055994530941723212145817656807550013436024 print repr(math.log(x)) 0.69314718055994529 x = 3 print logdec(x) 1.0986122886681096913952452369225257046474905578227 print repr(math.log(x)) 1.0986122886681098 x = 4 print logdec(x) 1.3862943611198906188344642429163531361510002687205 print repr(math.log(x)) 1.3862943611198906
Et, bien sûr, en utilisant en plus la fonction exponentielle, on devrait avoir log(e)=1
getcontext().prec = 50 x = 1 print logdec(expdec(x)) 1.0000000000000000000000000000000000000000000000000 print repr(math.log(math.exp(x))) 1.0
Pas mal, non?
def log10dec(x): """Calcul du logarithme décimal de x""" return logdec(x)/logdec(10)
Nous utiliserons la série suivante:
<m>exp(x) = sum{n=0}{infty}x_n_n</m>
Ce qui nous donne:
<m>exp(x) = 1 + {x}/{1!} + {x^2}/{2!} + {x^3}/{3!} + … + {x^n}/{n!} + …</m>
Bien entendu, nous n'essaierons pas de calculer les factorielles ni les puissances! Nous utiliserons une formule de récurrence, permettant de calculer chaque nouveau terme à partir du précédent:
Voilà le code:
def expdec(x): """Calcul de e à la puissance x""" if type(x)!=type(Decimal): x = Decimal(str(x)) s1 = 1 k = 1 n = 0 while True: n += 1 k *= x/n s2 = s1 + k if s2==s1: break else: s1 = s2 return s2
Validité: pour x quelconque.
Quelques vérifications:
getcontext().prec = 50 # Calcul de exponentielle de 3 x = 3 print expdec(x) 20.085536923187667740928529654581717896987907838557 print repr(math.exp(x)) 20.085536923187668
Comme les fonctions log et exp ont été calculés avec des séries différentes, on peut vérifier que log(exp(x)) redonne bien x:
getcontext().prec = 50 print logdec(expdec(1)) 1.0000000000000000000000000000000000000000000000000
On peut aussi calculer le nombre “e”:
getcontext().prec = 50 # Calcul du nombre e par calcul de l'exponentielle de 1 x = 1 print expdec(x) 2.7182818284590452353602874713526624977572470936998 print repr(math.exp(x)) 2.7182818284590451
La fonction existe dans le module décimal sous forme de méthode (x.sqrt()), mais je préfère la mettre sous forme de fonction:
def sqrtdec(x): """Calcul de la racine carrée de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) return x.sqrt()
Petit exercice amusant, puisque <m>sqrt{2} = {e^{0.5log(2)}}</m>:
getcontext().prec = 50 x = 2 print sqrtdec(x) 1.4142135623730950488016887242096980785696718753769 print expdec(Decimal("0.5")*logdec(x)) 1.4142135623730950488016887242096980785696718753769 print 2**(Decimal("1")/2) 1.4142135623730950488016887242096980785696718753769 print math.sqrt(x) 1.41421356237
Pas mal, non, pour des fonctions calculées à partir d'algorithmes différents?
Puisqu'on a la fonction “x à la puissance y” dans le module decimal, on va l'utiliser: racine n_ième de x est égal à x à la puissance (1/n).
def nrtdec(x, n): """Calcul de la racine n_ième de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) if type(n)!=type(Decimal): n = Decimal(str(n)) return x**(Decimal("1")/n)
Vérification:
getcontext().prec = 50 print nrtdec(27, 3) 3.0000000000000000000000000000000000000000000000000
Puisque la fonction “x à la puissance y” existe dans le module decimal, on va l'utiliser:
def powdec(x, y): """Calcul de x à la puissance y""" if type(x)!=type(Decimal): x = Decimal(str(x)) if type(y)!=type(Decimal): y = Decimal(str(y)) return x**y
Vérification en reprenant les variables pi et e définis sur cette page:
getcontext().prec = 50 print powdec(pi,e) 22.459157718361045473427152204543735027589315133997 print pi**e 22.459157718361045473427152204543735027589315133997 print expdec(e*logdec(pi)) 22.459157718361045473427152204543735027589315133997
Comme on a un bon algorithme pour calculer l'arc sinus, on va s'en servir!
Comme sin(pi/6) = 1/2, alors:
def pidec(): """Calcul de pi avec l'arc sinus""" return 6*asindec(0.5)
Simple, non?
On peut, bien entendu faire une fonction spéciale et autonome intégrant l'algorithme de l'arc sinus:
def pidec(): """Calcul de pi avec la précision courante (méthode avec arc sinus)""" getcontext().prec += 2 x = Decimal("0.5") xc = x*x k = x s1 = k n = 0 while True: n += 1 k *= xc*(2*n-1)*(2*n-1)/(2*n*(2*n+1)) s2 = s1 + k if s2 == s1: break s1 = s2 getcontext().prec -= 2 return 6*s2
Exemple:
getcontext().prec = 50 print pidec() 3.1415926535897932384626433832795028841971693993751
Vous pouvez vérifier avec ce site: http://www.brouty.fr/Maths/pi.html
Voilà la valeur de pi avec 1000 chiffres après la virgule que nous pouvons utiliser dans les calculs.
Si nécessaire, on peut toujours réduire le nombre de chiffres utiles en multipliant par le nombre entier 1: la valeur résultante est alors réduite à la valeur en-cours de getcontext().prec.
pi = Decimal("3.\ 14159265358979323846264338327950288419716939937510\ 58209749445923078164062862089986280348253421170679\ 82148086513282306647093844609550582231725359408128\ 48111745028410270193852110555964462294895493038196\ 44288109756659334461284756482337867831652712019091\ 45648566923460348610454326648213393607260249141273\ 72458700660631558817488152092096282925409171536436\ 78925903600113305305488204665213841469519415116094\ 33057270365759591953092186117381932611793105118548\ 07446237996274956735188575272489122793818301194912\ 98336733624406566430860213949463952247371907021798\ 60943702770539217176293176752384674818467669405132\ 00056812714526356082778577134275778960917363717872\ 14684409012249534301465495853710507922796892589235\ 42019956112129021960864034418159813629774771309960\ 51870721134999999837297804995105973173281609631859\ 50244594553469083026425223082533446850352619311881\ 71010003137838752886587533208381420617177669147303\ 59825349042875546873115956286388235378759375195778\ 18577805321712268066130019278766111959092164201989")
Si on veut calculer Pi plus rapidement et/ou avec plus de chiffres, on peut utiliser l'un des nombreux algorithmes que les fanas de Pi ont élaboré depuis toujours.
Par exemple, avec l'algorithme de Salamin et Brent:
def pidec(): """Calcul de pi avec la précision courante (Salamin et Brent)""" getcontext().prec += 2 a1 = Decimal("1.0") b1 = Decimal("1.0")/Decimal("2.0").sqrt() p1 = Decimal("1.0") t1 = Decimal("0.25") eps = Decimal("1.0e-" + str(getcontext().prec-1)) # condition d'arrêt while True: a2 = (a1 + b1)/2 b2 = (a1*b1).sqrt() p2 = 2*p1 a12 = a1-a2 t2 = t1-p1*a12*a12 if abs(a2-b2)<eps: break else: a1, b1, p1, t1 = a2, b2, p2, t2 ab = a2+b2 res = (ab*ab)/(4*t2) getcontext().prec -= 2 return res*1 # le *1 est pour revenir à la précision courante avec arrondi
Avec cet algorithme, on peut calculer:
On retrouve, bien entendu les mêmes résultats que ici: http://www.brouty.fr/Maths/pi.html
Facile, puisqu'on a la fonction exponentielle expdec(x) et que e = exponentielle de 1:
getcontext().prec = 50 print expdec(1) 2.7182818284590452353602874713526624977572470936998
Voilà la valeur de e avec 1000 chiffres après la virgule que nous pouvons utiliser dans les calculs.
Si nécessaire, on peut toujours réduire le nombre de chiffres utiles en multipliant par le nombre entier 1: la valeur résultante est alors réduite à la valeur en-cours de getcontext().prec.
e = Decimal("\ 2.718281828459045235360287471352662497757247093699\ 95957496696762772407663035354759457138217852516642\ 74274663919320030599218174135966290435729003342952\ 60595630738132328627943490763233829880753195251019\ 01157383418793070215408914993488416750924476146066\ 80822648001684774118537423454424371075390777449920\ 69551702761838606261331384583000752044933826560297\ 60673711320070932870912744374704723069697720931014\ 16928368190255151086574637721112523897844250569536\ 96770785449969967946864454905987931636889230098793\ 12773617821542499922957635148220826989519366803318\ 25288693984964651058209392398294887933203625094431\ 17301238197068416140397019837679320683282376464804\ 29531180232878250981945581530175671736133206981125\ 09961818815930416903515988885193458072738667385894\ 22879228499892086805825749279610484198444363463244\ 96848756023362482704197862320900216099023530436994\ 18491463140934317381436405462531520961836908887070\ 16768396424378140592714563549061303107208510383750\ 51011574770417189861068739696552126715468895703503\ 63")
(à vérifier)
def shdec(x): """Calcul du sinus hyperbolique de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) xc = x*x t = x s1 = t n = 0 while True: n += 1 t *= xc/(2*n*(2*n+1)) s2 = s1 + t if s2 == s1: break s1 = s2 return s2
(à vérifier)
def chdec(x): """Calcul du cosinus hyperbolique de x""" if type(x)!=type(Decimal): x = Decimal(str(x)) xc = x*x t = 1 s1 = t n = 0 while True: n += 1 t *= xc/((2*n-1)*2*n) s2 = s1 + t if s2 == s1: break s1 = s2 return s2
(à vérifier)
def thdec(x): """Calcul de la tangente hyperbolique de x""" return shdec(x)/chdec(x)
(à vérifier)
def ashdec(x): """Calcul de l'argument du sinus hyperbolique de x""" return logdec(x + sqrtdec(1+x*x))
(à terminer)
(à terminer)