Une histoire des logarithmes

Où l'on essaye de calculer les logarithmes…

Méthodes numériquesIPT Sup

Introduction — Prostaphérèse et intérêts composés

Les problèmes de complexité existent depuis bien avant l'invention de l'informatique. Lorsque l'on voulait faire des calculs compliqués à la main, cela pouvait valoir la peine de réfléchir un peu pour trouver des méthodes de calculs raisonnablement complexes.

Un domaine important d'amélioration des calculs concerne les multiplications, surtout lorsque l'on manipule des grands nombres (ou des nombres avec beaucoup de chiffres après la virgule). L'idée de transformer les multiplications en additions (beaucoup moins coûteuses) est en particulier une piste très intéressante. Une première méthode pour faire cela se nomme prostaphérèse et date de la fin du XVIème siècle. Elle repose sur l'utilisation d'identités trigonométriques comme :

cosacosb=12(cos(ab)+cos(a+b)){`\\cos a \\cos b = \\frac 1 2 \\bigl(\\cos(a - b) + \\cos(a + b)\\bigr)`}

Par exemple, on veut multiplier 1,4243{`1,4243`} par 8,2001{`8,2001`}. Notons tout d'abord que

1,4243×8,2001=100×0,14243×0,82001{`1,4243 \\times 8,2001 = 100 \\times 0,14243 \\times 0,82001`}

En cherchant dans les tables trigonométriques (on en disposait déjà à la fin de l'antiquité), on trouve que

0,14243=cos(81,8115...)et0,82001=cos(34,9142...){`0,14243 = \\cos (81,8115...) \\qquad \\mathrm{et} \\qquad 0,82001 = \\cos(34,9142...)`}

où les angles sont exprimés en degrés, les radians n'existant pas encore. On en déduit que :

On en déduit finalement que 1,4243×8,2001=11,6787{`1,4243 \\times 8,2001 = 11,6787`}, à comparer au résultat exact 11,67940243{`11,67940243`}. Ainsi, avec cette méthode, on trouve raisonnablement le produit de deux nombres en n'effectuant que des additions et des lectures de tables trigonométriques. Cette méthode sera encore en usage au milieu du XVIIème siècle, soit bien après l'invention des logarithmes.

Elle souffre cependant d'un certain nombre de faiblesses. En particulier, pour élever une nombre à puissance entière n{`n`}, on doit effectuer n{`n`} « aller-retours » avec les tables trigonométriques. Or, élever à une puissance n{`n`} est précisément le type d'opérations dont se servent les banquiers pour calculer les intérêts composés de leurs prêts et emprunts. Un autre problème de cette méthode est que l'on ne peut l'utiliser pour calculer des racines.

Les mathématiciens vont alors développer les fonctions logarithmes que vous connaissez bien. Dans ce TP, nous allons présenter un certain nombre de méthodes utilisées pour calculer les logarithmes, et nous intéresser en particulier aux questions de leur efficacité (on dirait maintenant complexité, mais il faut avoir à l'esprit que l'on compte ici des opérations effectuées à la main).

Pour visualiser la précision des calculs, on se basera sur une première fonction à programmer.

Question Écrire une fonction plpc (pour « plus long préfixe commun ») qui, étant donné deux chaînes de caractères, renvoie, comme son nom l'indique, le plus long préfixe commun aux deux chaînes. Par exemple, on veut le comportement suivant :

>>> plpc('Bonjour', 'Bonsoir')
'Bon'
>>> plpc('2.718281828', '2.718253968')
'2.7182'

Une fois cette fonction obtenue, on définira la fonction suivante :

def affiche(i, a, b) :
  stra = "{}".format(a)
  strb = "{}".format(b)
  print("{:2d} : {}".format(i, plpc(stra, strb)))

Un exemple d'utilisation est le suivant :

>>> np.pi
3.141592653589793
>>> 355/113
3.1415929203539825
>>> affiche(0, np.pi, 355/113)
 0 : 3.141592

Dans la suite, on s'intéresse donc à des méthodes permettant de calculer de manière approchée des réels (en l'occurence, ln2{`\\ln 2`}), et l'appel de la fonction affiche permettre de mettre en évidence les chiffres corrects obtenus. Dans l'exemple précédent, on compare la valeur de π{`\\pi`} à celle de 355/113{`355/113`}. On remarque donc que la fraction a en commun avec π{`\\pi`} les 6 premiers chiffres après la virgule.

La fonction prends trois arguments : le premier est un entier qui permet d'afficher une information comme le numéro de l'itération courante, et les deux arguments suivants sont des réels dont on affichera le plus grand préfixe commun de leur écriture décimale.

Remarques sur la fonction de formatage On utilise la fonction format (ou, plus précisement, la méthode format de la classe str) dans la fonction affiche.

Tout d'abord, pour convertir une flottant en chaîne de caractères. Cette façon de procéder est meilleure que la convertion directe à l'aide de la fonction str, car elle permet d'avoir plus de décimales :

>>> np.sqrt(2)
1.4142135623730951
>>> str(np.sqrt(2))
'1.41421356237'
>>> "{}".format(np.sqrt(2))
'1.4142135623730951'

Les accolades permettent d'indiquer où afficher les arguments de la fonction.

>>> "{} ou {} ?".format("Fromage", "dessert")
'Fromage ou dessert ?'

Enfin, dans la dernière ligne de la fonction affiche, les caractères entre les premières accolades donnent des informations supplémentaires sur l'affichage demandé (à droite de « : »). Ici, le « d » indique que l'on affiche un entier, et le « 2 » qui précède indique la largeur minimale à utiliser.

>>> "{}".format(3)
'3'
>>> "{:10d}".format(3)
'         3'

La Méthode de Briggs

Au début du XVIIème siècle, Henry Briggs publia une table de logarithmes. Ayant remarqué qu'une fonction logarithmique log{`\\log`} (quelqu'en soit la base) vérifie les deux propriétés suivantes :

  • loga=12loga{`\\log \\sqrt{a\\,} = \\frac 1 2 \\log a`},
  • si a{`a`} est très proche de 0{`0`}, alors log(1+a){`\\log (1 + a)`} est proportionnel à a{`a`} (le coefficient dépend de la base du logarithme),

il procéda ainsi :

  1. on pose u0=2{`u_0 = 2`} et v0=10{`v_0 = 10`} ;
  2. on en calcule les racines carrées successives, soit en écriture moderne :
    nN, un+1=un et vn+1=vn{`\\forall \\, n \\in \\mathbf N, \\ u_{n+1} = \\sqrt{u_n} \\ \\mathrm{et} \\ v_{n+1} = \\sqrt{v_n}`}
  3. le rapport un1vn1{`\\dfrac {u_n - 1}{v_n - 1}`} converge alors vers ln2ln10=log102{`\\dfrac {\\ln 2}{\\ln 10} = \\log_{10} 2`}.

Briggs itérera le processus 54 fois (autrement dit, 54 fois deux calculs de racines carrées faits à la main).

Petite explication mathématique On a un=212n1{`u_n = 2^{\\frac 1 {2^n}} \\rightarrow 1`}, donc

12nln2=lnunnun1{`\\frac 1 {2^n} \\ln 2 = \\ln u_n \\sim_{n \\rightarrow \\infty} u_n - 1`}

et de même,

vn1n12nln10{`v_n - 1 \\sim_{n \\rightarrow \\infty} {\\frac 1 {2^n}} \\ln 10`}

On en déduit que l'on a bien

un1vn1nln2ln10{`\\frac {u_n - 1} {v_n - 1} \\rightarrow_{n \\rightarrow \\infty} \\frac {\\ln 2}{\\ln 10}`}

Question Écrire une fonction qui calcule une valeur approchée de log10(2){`\\log_{10}(2)`} à l'aide de cette méthode. Quelques remarques supplémentaires :

  • La fonction utilise le module numpy renommé en np. Il faudra donc avoir entré auparavant la commande import numpy as np ;
  • La valeur « exacte » de log10(2){`\\log_{10}(2)`}, calculée par l'expression np.log(2) / np.log(10) ne sert qu'à afficher les décimales correctes des résultats obtenus ;
  • La racine carrée d'un flottant s'obtient avec la fonction np.sqrt.
  • Briggs a itéré le processus 54 fois. Essayez de retrouver le résultat qu'il a obtenu.

Le programme à compléter est :

...
for i in range(55) :
  affiche(i, np.log(2) / np.log(10), ...)
  ...

Les points de suspension sont, bien sûr, ce qu'il faut modifier.

Le début de l'affichage devra ressembler à ceci :

 0 : 0.
 1 : 0.
 2 : 0.
 3 : 0.
 4 : 0.
 5 : 0.
 6 : 0.
 7 : 0.
 8 : 0.
 9 : 0.30
10 : 0.30
11 : 0.30
12 : 0.30
13 : 0.30
14 : 0.3010

Normalement, lors de l'exécution du programme (a priori correct) précédent, nous n'arrivez pas à faire aussi bien que Briggs. Le nombre de décimales correctes augmente jusqu'à la 27ème itération, puis diminue.

Question Pourquoi ?

De plus, à partir de la 54ème itération, même la partie 0. n'est plus affichée. Il est même possible d'avoir un message d'erreur. La raison est, à nouveau, la taille finie de la mantisse. Pour illustrer cela, on a représenté dans le tableau suivant les racines itérées de 1.4{`1.4`}, avec la représentation en binaire. On voit clairement qu'entre le 1 de gauche et le 1 suivant, il y a des 0 dont le nombre augmente d'un à chaque étape.

DécimalBinaire
1.400001.011001100110
1.183221.001011101110
1.087761.000101100111
1.042961.000010101111
1.021251.000001010111
1.010571.000000101011
1.005271.000000010101
1.002641.000000001010
1.001321.000000000101

Question Expliquer l'origine du problème évoqué.

Question

  1. Montrer que si a1{`a - 1`} est compris entre 0{`0`} et 2n{`2^{-n}`}, alors a1{`\\sqrt {a\\,} - 1`} est compris entre 0{`0`} et 2(n+1){`2^{-(n+1)}`} (on pourra utiliser la contraposée, soit si a12n{`a - 1 \\geq 2^{-n}`}, alors a2121n{`a^2 - 1 \\geq 2^{1-n}`}).
  2. En déduire (approximativement) la taille de la mantisse utilisée pour représenter les nombres flottants.

Remarque La taille de la mantisse s'obtient avec la commande suivante :

>>> np.finfo(float).nmant
52

L'aire sous une courbe

Les calculs précédents faisant intervenir les logarithmes en base 10{`10`}. Cependant, on découvrit un peu plus tard que l'aire sous la courbe de la fonction x1x{`x \\mapsto \\frac 1 x`} est aussi une fonction logarithme, que l'on appelle maintenant logarithme népérien. Cette fonction logarithmique est la plus pertinente mathématiquement.

Nous allons programmer une méthode pour calculer ln2{`\\ln 2`} de façon approchée à l'aide de l'intégrale de la fonction inverse entre 1{`1`} et 2{`2`}, et en utilisant la bibliothèque numpy.

Question Écrire une fonction rect qui calcule l'intégrale demandée à l'aide de la méthode des rectangles « à gauche. » Elle prendra comme argument un entier indiquant le nombre de rectangles à utiliser.

C'est bien connu, la vitesse de convergence pour la méthode des rectangles est de l'ordre de 1n{`\\frac 1 n`}. Pour en voir une illustration, exécutez la commande suivante.

for i in range(8) :
  affiche(i, np.log(2), rect(10 ** i))

Question À quoi sert le 10 ** i ?

Ainsi, à la dernière itération, l'intervalle [1,2]{`[1, 2]`} est découpé en 107=10000000{`10^7 = 10000000`} rectangles. On a approximativement une nouvelle décimale correcte à chaque fois que l'on multiplie le nombre de subdivisions par 10{`10`}.

On connaît une méthode plus efficace pour calculer des intégrales, puisqu'elle converge en 1n2{`\\frac 1 {n^2}`}. Tâchons de l'utiliser.

Question Écrire une fonction trap qui calcule la même intégrale à l'aide de la méthode des trapèzes.

Expérimentons cette fonction à l'aide de la commande :

for i in range(8) :
  affiche(i, np.log(2), trap(10 ** i))

Deux remarques s'imposent :

  • Difficile d'avoir plus d'une dizaine de décimales, alors que 52 bits de mantisse correspondent plutôt à une quinzaine.
  • Le nombre d'opérations pour obtenir une précision décente reste vraiment très élevé.

Calculs en séries

Les progrés mathématiques du XVIIème vont permettre de développer de nouvelles techniques bien plus efficaces, notamment à l'aide des séries numériques.

Série de Mercator

En 1668, dans son traîté «Logarithmotechnia», Nicolas Mercator (connu pour ses projections cartographiques) démontre la formule suivante (écrite de façon moderne) :

x]1,1], ln(1+x)=n=1(1)n+1xnn{`\\forall \\, x \\in \\mathopen]-1,1], \\ \\ln (1 + x) = \\sum_{n = 1}^\\infty \\frac {(-1)^{n+1} x^n} n`}

Avec des points de suspension, cela donne :

ln(1+x)=xx22+x33x44+x55x66+ {`\\ln (1 + x) = x - \\frac {x^2} 2 + \\frac {x^3} 3 - \\frac {x^4} 4 + \\frac {x^5} 5 - \\frac {x^6} 6 + \\ \\cdots`}

Pour x=1{`x = 1`}, on reconnaît la série harmonique alternée :

ln2=112+1314+1516+  =n=0(1)nn+1{`\\ln 2 = 1 - \\frac 1 2 + \\frac 1 3 - \\frac 1 4 + \\frac 1 5 - \\frac 1 6 + \\ \\cdots \\ = \\sum_{n = 0}^\\infty \\frac {(-1)^n} {n + 1}`}

Question Écrire une fonction mercator1 qui prends en argument un entier n et retourne la somme partielle précédente évaluée jusqu'à ce rang n{`n`}.

Question Pour évaluer son (in-)efficacité, modifier la fonction précédente pour qu'elle affiche (à l'aide de la bien-nommée fonction affiche que vous devez bien connaître maintenant) les valeurs approchées de ln2{`\\ln 2`} obtenues à chaque itération.

On peut être un peu déçu par cette méthode de calcul. Mais la formule précédente s'applique aussi à d'autres valeurs de x{`x`}. Pour le calcul de ln2{`\\ln 2`}, on peut remarquer que l'on a :

ln2=ln12=ln(112).{`\\ln 2 = - \\ln \\frac 1 2 = - \\ln \\Bigl(1 - \\frac 1 2\\Bigr).`}

Question Écrire une fonction mercator2 basée sur le calcul de ln2{`\\ln 2`} vu comme ln(112){`- \\ln \\bigl(1 - \\frac 1 2\\bigr)`}, et visualiser l'efficacité de cette nouvelle méthode à l'aide d'un judicieux appel à la fonction affiche.

Normalement on a, enfin, une méthode de calcul raisonnablement efficace.

Série de Gregory

Une amélioration simple pour l'utilisation de la série de Mercator est dûe à James Gregory, qui utilise intelligemment les propriétés du logarithme. Il soutrait les séries correspondant à ln(1+x){`\\ln (1 + x)`} à ln(1x){`\\ln (1 - x)`}, pour obtenir :

x]1,1[, ln(1+x1x)=2x+2x33+2x55+2x77+  =2n=0x2n+12n+1{`\\forall \\, x \\in \\mathopen]-1, 1\\mathclose[, \\ \\ln \\Bigl(\\frac {1+x}{1-x}\\Bigr) = 2 x + 2 \\frac {x^3} 3 + 2 \\frac {x^5} 5 + 2 \\frac {x^7} 7 + \\ \\cdots \\ = 2 \\sum_{n = 0}^\\infty \\frac {x^{2n + 1}}{2n + 1}`}

Question Pour quelle valeur de x{`x`} peut-on calculer ln2{`\\ln 2`} ?

Question Écrire une fonction gregory correspondant à cette méthode, et visualiser son efficacité.