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:

• http://fr.wikipedia.org/wiki/S%C3%A9rie_enti%C3%A8re

• http://serge.mehl.free.fr/anx/dev_lim.html

• etc…

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:

• http://www.brouty.fr/Maths/pi.html

Et en voilà 1 million:

• http://www.piday.org/million.php

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:

• l'importation du module “decimal” avec “from decimal import *”

• que les calculs se feront avec la précision donnée par “getcontext().prec=nbc” (nbc=nb de chiffres significatifs, 28 par défaut)

• que les divisions avec '/' seront décimales, ce qui, avec Python 2.6, nécessite la ligne "from __future__ import division" juste après le shebang.

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:

• premier terme pour n=0: <m>terme (0) = x</m>

• termes suivants: <m>terme (n) = {terme (n-1)}*x_2_2n_2n_1</m>

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!

• Test réalisé: 50000 tirages au hasard de valeurs de x pris de 0 à 2pi avec calcul sur 28 chiffres: env. 3/1000 de secondes par appel de la fonction!

• Et avec les mêmes calculs sur 56 chiffres (soit le double): env. 6/1000 de secondes

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:

• premier terme pour n=0: <m>terme(0) = 1</m>

• termes suivants: <m>terme(n) = {terme(n-1)}*x_2_2n-1_2n</m>

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>:

• si n=0 ⇒ <m>a_n = 1</m>

• sinon ⇒ <m>a_n = {prod{k=1}{n}{(2k-1)}}/{prod{k=1}{n}{2x}}</m> soit <m>a_n = {1.3.5…(2n-1)}/{2.4.6…(2n)}</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:

• premier terme pour n=0: <m>terme(0) = x</m>

• termes suivants: <m>terme(n) = {terme(n-1)}*x_2_2n-1_2n-1_2n_2n_1</m>

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:

• premier terme pour n=0: <m>terme(0) = x-1_x_1</m>

• termes suivants: <m>terme(n) = {terme(n-1)}*{ {(x-1)^2(2n-1)}/{(x+1)^2(2n+1)} }</m>

• et ne pas oublier de multiplier par 2 à la fin!

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:

• premier terme pour n=0: <m>terme(0) = 1</m>

• termes suivants: <m>terme(n) = {terme(n-1)}*x_n</m>

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:

• Pi avec 1000 décimales en 0.14 secondes (au lieu de 7 secondes avec l'arc sinus)

• Pi avec 5000 décimales en 4 secondes

• Pi avec 50000 (cinquante mille) décimales en 10 minutes environ

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)