Calcul de Dérivées via les nombres duaux

When teaching learns fun topics

Note: This post was originally written for engineering students asking why learning programming matters even if they were not interested by computer science. The C language was imposed, sadly. The attempt here was to explain a concrete example using basic maths (high school) and then use it for solving physical problems from basic mechanics or eletronics. The latter never came to light, yet.

Le travail d'un ingénieur est de résoudre des problèmes et pour cela il puise dans sa boite à outils. L'idée de ce document est donc de montrer comment on construit un code pour répondre à une question issue d'un problème d'ingénierie.

Nous utilisons le langage C pour illustrer même si de notre point de vue ce n'est pas le meilleur langage pour rapidement concrétiser ses idées. La suite logicielle Python/Numpy/Scipy ou le langage Julia nous apparaissent plus pertinents.

L'important dans la résolution d'un problème est de structurer ses idées ; en d'autres termes l'important est de penser en terme d'Algorithmique et de Programmation (voir ici pour une référence complète).


Nous proposons de calculer numériquement la dérivée d'une fonction. Pour cela, nous allons implémenter une différentiation automatique basée sur une accumulation en avant (forward accumulation) utilisant les nombres duaux.

La différentiation automatique est utilisée dans les réseaux neuronaux et autre technologie d'apprentissage semi-automatique (machine learning). Cependant, en général c'est la méthode accumulation en arrière (reverse accumulation ou backpropagation) qui est utilisé ; voir ou .

1 Le problème

Soit une fonction \(f\) allant des nombres réels dans les nombres réels,

\begin{equation*} \begin{array}{lclcl} f &:& \mathbb{R} &\longrightarrow& \mathbb{R}\\ & & x &\mapsto& y=f(x) \end{array} \end{equation*}

et sa dérivée notée \(f^\prime\) est définie comme,

\begin{equation} \label{org813a013} f^\prime(x) = \underset{h\rightarrow 0}{\lim} \frac{f(x+h) - f(x)}{h}. \end{equation}

Pour une valeur donnée, disons \(a\), il est facile d'évaluer la fonction en cette valeur et nous obtenons la nouvelle valeur \(b\) ; explicitement \(b=f(a)\).

Si la fonction \(f\) est analytique, alors nous pouvons calculer symboliquement (avec un papier et un crayon ou SymPy) la fonction dérivée \(f^\prime\) et nous obtenons une formule qui peut être utilisée pour évaluer en la valeur \(a\).

Mais dans la plupart des cas, la fonction \(f\) ne peut pas être dérivée symboliquement.

Comment faire pour évaluer la dérivée en \(a\) ? Comment calculer \(f^\prime(a)\) ?

1.1 Différence finie

En utilisant la définition de la dérivée \eqref{org813a013}, la première stratégie serait de fixer une valeur \(h\) considérée petite comparée à la valeur \(a\). Puis d'approcher le calcul dérivée par,

\begin{equation*} \delta = \frac{f(x+h) - f(x)}{h}. \end{equation*}

La question est alors le choix de cette valeur \(h\).

1.1.1 Exemple

Considérons le cas simple de la fonction \(f(x)=x^2\). La dérivée est trivialement donnée par \(f^\prime(x)=2x\). Choisissons quelques valeurs:

\(a\) \(h\) \(\delta\) \(f^\prime(a)\) err. rel. %
2.1 0.1 4.3 4.2 2.38
2.1 0.01 4.21 4.2 0.23
0.3 0.1 0.7 0.6 16.7
0.3 0.01 0.61 0.6 1.7

On voit clairement que cette approche n'est pas stable et l'erreur varie fortement en fonction des valeurs choisies. Cette erreur s'explique principalement par le développement de Taylor.

Cette approche pourrait être améliorée en utilisant des ordres plus élevés du développement de Taylor. Nous choisissons une autre méthode.

1.2 Méthode des nombres duaux

1.2.1 Maths élémentaires

Revenons aux principes élémentaires de dérivation et aux 4 opérations arithmétiques relatives à la dérivation:

\begin{equation} \label{org06f591c} \begin{array}{lcl} (f+g)^\prime &=& f^\prime + g^\prime\\ (f*g)^\prime &=& f^\prime * g + f* g^\prime\\ (f-g)^\prime &=& f^\prime - g^\prime\\ \displaystyle \left(\frac{f}{g}\right)^\prime &=& \displaystyle \frac{f^\prime * g - f* g^\prime}{g^2}. \end{array} \end{equation}

1.2.2 Nombres duaux

Définissons un nouveau type de nombre composé de 2 réels. Soit \(x\) et \(y\) deux de ces nouveaux nombres qui s'écrivent,

\begin{equation*} \begin{array}{lcl} x &=& (v, d) \\ y &=& (u, p) \end{array} \end{equation*}

et définissons comment ajouter, multiplier, soustraire et diviser ses nombres,

\begin{equation} \label{org7de21a8} \begin{array}{lcll} \texttt{add}(x, y) &=& (v+u ,& d+p) \\ \texttt{mul}(x, y) &=& (u*v ,& d*v ~+~ u*p) \\ \texttt{sub}(x, y) &=& (u-v ,& d-p) \\ \texttt{div}(x, y) &=& (u*v ,& \displaystyle \frac{d*v ~-~ u*p}{p^2}). \end{array} \end{equation}

Nous notons ce nombre \((\texttt{value}, \texttt{deriv})\).

  1. La partie value se comporte comme les nombres réels avec les opérations arithmétiques usuelles.
  2. La partie deriv se comporte avec les opérations arithmétiques de la dérivation \eqref{org06f591c}. Intéressant!

1.2.3 Parallèle avec les nombres complexes

Un nombre complexe peut s'écrire \(z = x + iy\) et le nombre imaginaire \(i\) est défini tel que \(i^2=-1\). À partir de là, en utilisant les règles élémentaires de calculs et les opérations arithmétiques usuelles, il est facile d'obtenir par exemple,

\begin{equation*} z_1 \times z_2=(x_1*x_2 - y_1*y_2) + i(x_1*y_2 + y_1*x_2) \end{equation*}

ce que nous écririons:

\begin{equation*} \texttt{cadd}(z_1, z_2) = (x_1*x_2 - y_1*y_2,~~ x_1*y_2 + y_1*x_2). \end{equation*}

De manière similaire, un nombre dual peut s'écrire \(z = x + \varepsilon y\) et \(\varepsilon\) est défini tel que \(\varepsilon^2=0\). Par conséquent, il est facile d'obtenir \eqref{org7de21a8}, par exemple,

\begin{equation*} (x_1 + \varepsilon y_1)\times(x_2 + \varepsilon y_2) = (x_1*x_2) + \varepsilon(x_1*y_2 + y_1*x_2). \end{equation*}

1.2.4 Exemple

Reprenons l'exemple \(f(x)=x^2\) et définissons le nombre dual \(z=(a, 1.)\). Appliquons les opérations avec \(z\) au lieu de \(x\) dans \(f\), cela donne,

\begin{equation*} \texttt{mul}(z,z) = (a*a, 1.*a + a*1.) = (a*a, 2*a). \end{equation*}

Donc nous avons évalué la fonction en la valeur \(a\) et nous avons aussi obtenu l'évaluation de la dérivée en cette valeur.

Pourquoi le 1. dans \((a, 1.)\) ? Si nous considérons la fonction \(f(x)=x\) alors en substituant \(z\) à la place de \(x\) dans \(f\), nous obtenons bien que la dérivée vaut 1.

Avec le même argument, une constante \(c\) se définit comme \((c, 0.)\) puisque sa dérivée est nulle.

Il faut donc distinguer ce qui représente une variable au sens mathématique comme \(x\) et une constante au sens mathétmatique comme \(c\). La dérivation se fait par rapport à cette variable \(x\).


Maintenant, comment programme-t-on cela ?


2 Implémentations

La version la plus naïve d'un ordinateur ne comprend que les 4 opérations arithmétiques sur les types int et float. À partir de cela nous allons construire tout ce dont nous avons besoin pour répondre à notre problème: comment calculer en une valeur la dérivée d'une fonction.

2.1 Commençons

2.1.1 Représenter un nombre dual   v1

Il apparaît assez clair que nous allons représenter un nombre dual avec struct. Le plus simple est:

typedef struct {
    char name[MAXCHAR];
    float value;
    float deriv;
} dual;

Comme attendu, la structure possède 2 champs: l'un pour représenter les valeurs value et l'autre pour représenter la dérivée deriv. Nous ajoutons un champ name pour associer un nom en espérant faire de plus jolis affichages.

Pour tester, nous écrivons le corps de notre programme et l'affichage que nous obtenons.

/* provide printf/scanf and other */
#include <stdio.h>

/* Length maximum for some "string" (array of char) */
#define MAXCHAR 1000

typedef struct {
    char name[MAXCHAR];
    float value;
    float deriv;
} dual;

int main () {
    printf("Hi");

    /* the value we are interrested in */
    float a = 2.1;

    /* declare a dual number */
    dual z;

    /* fill the fields                                 */
    /*                                                 */
    /* sprintf is similar to printf                    */
    /*  instead of output to the console               */
    /*  sprintf store the output to its first argument */
    /*     (here z.name)                               */
    sprintf(z.name, "%s", "LE nom bien");
    z.value = a;
    z.deriv = 1.;

    /* display the dual number */
    printf("\n");
    printf("\tName : %s\n", z.name);
    printf("\tValue: %f\n", z.value);
    printf("\tDeriv: %f\n", z.deriv);

    printf("Bye.");
    return 0;
  }
Hi
	Name : LE nom bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

2.1.2 Fonction pour définir une variable   v2

Nous souhaitons avoir une fonction qui prend en argument un nombre réel (float) et une chaîne de caractère (char*) et retourne un nombre dual (return z; avec la variable z de type dual). Nous appelons cette fonction newvar et sa signature est donc:

dual newvar(float, char*);

Son implémentation est aussi directe puisque nous avions déjà écrit cela dans le programme précédent.

dual newvar(float val, char* avoile) {
    dual z;
    sprintf(z.name, "%s", avoile);
    z.value = val;
    z.deriv = 1.;
    return z;
}

WARNING: Les quantités val et z sont des variables au sens informatique, et respectivement de type float et dual. Cependant, la variable informatique val ne modélise pas la notion de variable au sens mathématique, c'est-à-dire qu'au nombre val n'est pas associé une dérivée. Nous construisons justement un code pour cela et le nouveau type dual sera une modélisation informatique plus proche de la notion mathématique de variable en considérant la dérivation.

Ainsi nous adaptons notre programme pour obtenir le même affichage que précédemment.

Hi
	Name : LE nome bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

Bien évidemment, nous pouvons aussi définir une fonction similaire à newvar qui crée une constante, disons newcst. Laissé en exercice, sinon regardez dans le fichier source.

Hi
	Name : LE nome bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

2.1.3 Fonction pour afficher un nombre dual   v3

Nous souhaitons avoir une fonction qui prend en argument un nombre dual et ne retourne rien (void). Sa signature est,

void print(dual);

Et nous l'avons déjà implémentée.

void print(dual z) {
    printf("\n");
    printf("\tName : %s\n", z.name);
    printf("\tValue: %f\n", z.value);
    printf("\tDeriv: %f\n", z.deriv);
}
Hi
	Name : LE nome bien
	Value: 2.100000
	Deriv: 1.000000
Bye.

Au lieu d'afficher à chaque fois le programme complet, nous allons maintenant—pour des raisons de lisibilité—afficher uniquement les modifications que nous avons apportées.

--- Supplementary/main-v2b.c    2022-04-19 21:45:57.483116926 +0000
+++ Supplementary/main-v3.c     2022-04-19 21:45:57.484116996 +0000
@@ -14,6 +14,7 @@
       /* defined after the main    */
 dual newvar(float, char*);
 dual newcst(float, char*);
+void print(dual);

 int main () {
     printf("Hi");
@@ -29,10 +30,7 @@
     z = newvar(a, "LE nome bien");

     /* display the dual number */
-    printf("\n");
-    printf("\tName : %s\n", z.name);
-    printf("\tValue: %f\n", z.value);
-    printf("\tDeriv: %f\n", z.deriv);
+    print(z);

     printf("Bye.");
     return 0;
@@ -56,3 +54,10 @@
     z.deriv = 0.;
     return z;
 }
+
+void print(dual z) {
+    printf("\n");
+    printf("\tName : %s\n", z.name);
+    printf("\tValue: %f\n", z.value);
+    printf("\tDeriv: %f\n", z.deriv);
+}

2.2 Opérons sur 2 nombres duaux

2.2.1 Multiplication   v4a

La signature de la fonction qui permet la multiplication entre deux nombres duaux est assez naturelle. Cette fonction a 2 arguments qui sont des nombres duaux (dual) et elle retourne un nombre dual (return z; avec z une variable de type dual).

dual mul(dual, dual);

Et l'implémentation est donnée par \eqref{org7de21a8}.

dual mul(dual z1, dual z2) {
    dual z;
    sprintf(z.name, "(%s*%s)", z1.name, z2.name);

    z.value = z1.value * z2.value;
    z.deriv = (z1.deriv * z2.value) + (z1.value * z2.deriv);

    return z;
}

Pour tester cette nouvelle fonction, nous écrivons le programme principal ci-dessous, que nous exécutons et nous obtenons bien le résultat attendu.

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : (x*x)
	Value: 4.409999
	Deriv: 4.200000
Bye.

Nous obtenons bien pour la fonction \(x^2\) évaluée en 2.1 la valeur \(2.1^2=\) 4.41 (à l'arrondi près) et sa dérivée \(2\times 2.1=\) 4.2.

--- Supplementary/main-v3.c     2022-04-19 21:45:57.484116996 +0000
+++ Supplementary/main-v4a.c    2022-04-19 21:45:57.486117136 +0000
@@ -16,6 +16,8 @@
 dual newcst(float, char*);
 void print(dual);

+dual mul(dual, dual);
+
 int main () {
     printf("Hi");

@@ -27,11 +29,14 @@


     /* fill the fields */
-    z = newvar(a, "LE nome bien");
+    z = newvar(a, "x");

     /* display the dual number */
     print(z);

+    dual zz = mul(z, z);
+    print(zz);
+
     printf("Bye.");
     return 0;
   }
@@ -61,3 +66,13 @@
     printf("\tValue: %f\n", z.value);
     printf("\tDeriv: %f\n", z.deriv);
 }
+
+dual mul(dual z1, dual z2) {
+    dual z;
+    sprintf(z.name, "(%s*%s)", z1.name, z2.name);
+
+    z.value = z1.value * z2.value;
+    z.deriv = (z1.deriv * z2.value) + (z1.value * z2.deriv);
+
+    return z;
+}

Sur le même principe, il est facile d'implémenter les opérations manquantes add, sub et div. Laissé en exercice, sinon regardez dans le fichier source.

2.2.2 Vérifications des 4 opérations   v4b

Nous commençons pour définir une valeur variable (newvar). C'est une variable informatique arbitrairement nommée z et d'autre part nous lui attribuons le nom "x" dans notre modélisation des nombres duaux.

Nous définissons aussi une valeur constante (newcst). C'est une variable informatique arbitrairement nommée c et d'autre part nous lui attribuons le nom "c" dans notre modélisation des nombres duaux.

Puis nous définissons un programme principal qui affiche:

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : c
	Value: 1.200000
	Deriv: 0.000000

	Name : (x*c)
	Value: 2.520000
	Deriv: 1.200000

	Name : (x+x)
	Value: 4.200000
	Deriv: 2.000000

	Name : (c-x)
	Value: -0.900000
	Deriv: -1.000000

	Name : ((x+c)/(c-x))
	Value: -3.666667
	Deriv: 2.962964
Bye.

À partir de ces exemples, nous pouvons tester si notre implémentation est correcte. Il faut faire plusieurs tests pour s'assurer que tout fonctionne comme espéré.

Nous commençons à voir l'intérêt d'avoir un champ name dans notre modélisation des nombres duaux.

2.3 Test avec une fonction plus compliquée

Nous souhaitons vérifier que notre implémentation fait bien les calculs escomptés. Pour cela nous voulons définir une fonction suffisamment compliquée et en même temp nous voulons pouvoir facilement vérifier le résultat.

Une fonction toute indiquée est la fonction exponentielle car sa dérivée est elle-même,

\begin{equation*} \left(e^x\right)^\prime = e^x \end{equation*}

Cependant, l'ordinateur ne sait pas implicitement comment calculer la fonction exponentielle. Qu'à cela ne tienne!

Nous avons besoin d'être capable d'évaluer la fonction exponentielle à partir uniquement des 4 opérations arithmétiques. Or la fonction exponentielle s'écrit aussi sous la forme d'une série convergente (et partout!),

\begin{equation*} e^x = \sum_{k=0}^\infty \frac{x^k}{k!} \end{equation*}

Par simplicité ici, nous tronquerons arbitrairement la somme à 20 termes.

2.3.1 Fonction exponentielle   v5

Il est assez clair que notre fonction exponential aura comme argument un nombre dual et retournera un nombre dual, donc sa signature est,

#define NEXP 20
dual exponential(dual);

Et son implémentation est simple mais mérite un peu d'attention. La somme \(\sum\) va être naturellement transformée en boucle for. Par ailleurs, cette somme commence par k=0 mais comme nous allons accumuler son résultat, nous allons initialiser l'accumulateur (variable informatique expx) à la valeur k=0 et faire commencer la boucle à k=1.

dual exponential(dual x) {

    /* temporary variables */
    dual expx  = newcst(1., "");
    dual xk    = newcst(1., "");
    dual kfac  = newcst(1., "");

    float k;
    for (k=1; k<=NEXP; k++) {
        dual kk = newcst(k, "");

        xk   = mul(xk, x);      /* x, x*x, x*x*x, etc. => eval x^k */
        kfac = mul(kfac, kk);   /* 1, 1*2, 1*2*3, etc. => eval k!  */

        dual term_k = div(xk, kfac);

        /* update the Sum with another term */
        expx = add(expx, term_k);

        /* to avoid segmentation fault (because of ugly C)   */
        /* instead of accumalating all the ( and ) and +,*,/ */
        sprintf(expx.name, "");
        sprintf(xk.name, "");
        sprintf(kfac.name, "");
        sprintf(kk.name, "");
    }

    /* pretty print the name */
    sprintf(expx.name, "~(e^(%s))", x.name);
    return expx;
}

Avec un papier et un crayon, et en écrivant les premiers termes, il est rapidement clair pourquoi expx, xk et kfrac sont initialisés avec la fonction newcst et non la fonction newvar.

Finalement, nous écrivons un programme principal pour tester.

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : c
	Value: 1.200000
	Deriv: 0.000000

	Name : ~(e^(x))
	Value: 8.166168
	Deriv: 8.166168

	Name : ~(e^((((c*x)*x)+c)))
	Value: 659.838440
	Deriv: 3325.548828
Bye.

Il est clair que nous obtenons bien le bon résultat, même pour une fonction non-triviale. Notre implémentation n'utilise que les 4 opérations arithmétiques.

2.3.2 Retour des nombres complexes

Nous pourrions facilement étendre cela à toutes les fonctions (sinus, cosinus, tangente, etc.). Par exemple, nous pourrions calculer la fonction sinus en remarquant que le sinus est la partie imaginaire d'une fonction exponentielle complexe. Par conséquent, nous aurions besoin d'implémenter une structure représentant les nombres complexes puis d'implémenter les 4 opérations arithmétiques pour ce nouveau type complex. Et enfin modifier notre structure dual en remplaçant le type float par ce type complex.

Cette amélioration utilisant les nombres complexes est un bon exercice!

2.3.3 Autre fonction

Nous avons choisi la fonction exponentielle mais nous aurions aussi pu choisir la fonction racine carrée comme illustration. Elle se calcule par la méthode de Héron comme,

\begin{equation*} \left\{ \begin{array}{lcl} x_{0} &=& a \\ x_{n+1} &=& \displaystyle 0.5 \left( x_n + \frac{a}{x_n} \right) \end{array} \right. \end{equation*}

Son implémentation et la vérification est un bon exercice!

2.4 Pour aller plus loin

Le langage C possède une biobliothèque standard définisssant des fonctions mathématiques: math.h. Nous souhaiterions utiliser cette bibliotèque pour calculer des fonctions connues plutôt que de tout re-implémenter par nous-même.

2.4.1 Exemple: fonction exponentielle (encore)   v6

Pour illustrer notre propos, étendons la fonction exp de math.h. Ou pour le dire autrement, nous allons utiliser la fonction exp de math.h avec notre type dual.

En sachant que la fonction exp de la bibliothèque math.h a la signature:

double exp(double);

or double et float sont des types compatibles. Nous voulons une signature similaire pour le type dual ; en fait la même signature que la fonction exponential:

#include <math.h>
dual dexp(dual);

Nous choisissons arbitraiment le nom dexp, et il nous semble parlant: l'extension de exp au type dual.

Pour l'implémentation, il faut commencer par se rappeler la définition mathématique:

\begin{equation*} \left(e^u\right)^\prime = u^\prime e^u \end{equation*}

Ce qui se traduit naturellement par:

dual dexp(dual u) {
    dual z;
    sprintf(z.name, "(e^(%s))", u.name);

    z.value = exp(u.value);
    z.deriv = u.deriv * z.value;

    return z;
}

Il faut noter que la multiplication u.deriv * z.value se fait dans le type float. Finalement, nous comparons nos 2 versions.

Hi
	Name : x
	Value: 2.100000
	Deriv: 1.000000

	Name : c
	Value: 1.200000
	Deriv: 0.000000

# by +,*,-,/

	Name : ~(e^((((c*x)*x)+c)))
	Value: 659.838440
	Deriv: 3325.548828

# by math.h

	Name : (e^((((c*x)*x)+c)))
	Value: 659.841492
	Deriv: 3325.601074
Bye.

Sur ce même principe, nous pouvons étendre toutes les fonctions qui sont dans la bibliothèque math.h et dont nous connaissons les formules de dérivations, comme sinus, cosinus, tangente, etc.

2.4.2 Dérivée de toute fonction

Avec notre code, nous sommes maintenant capable de calculer la dérivée de n'importe quelque fonction. Même des fonctions qui sont impossibles à dériver à la main. Par exemple,

dual fweird(dual x) {
    dual y;
    sprintf(y.name, "%s", x.name);
    y.deriv = x.deriv;
    if (x.value < 0) {
        return x;
    } else if (((int)x.value) % 2 == 0) { /* even */
        y.value = x.value - 1.;
        return add(div(sub(x, newcst(1, "1")),
                       add(x, newcst(2, "2"))),
                   fweird(y));
    } else {                              /* odd   */
        y.value = x.value - 1.;
        return mul(x, fweird(y));
    }
}
Hi
	Name : x
	Value: 12.300000
	Deriv: 1.000000

	Name : (((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+(x*(((x-1)/(x+2))+x)))))))))))))
	Value: -10758.750000
	Deriv: 6566.896484
Bye.

2.4.3 Intervalle de définitions des nombres

Pour connaître la taille mémoire utilisée pour représenter un nombre entier en mémoire, le plus simple est d'utiliser sizeof(int) (retourne en octet/byte). En général, ceci est de 4 octets soit 32 bits. Seul les entiers dans l'intervalle \(\left[-2^{31} ~;~ 2^{31} - 1\right]\) sont représentables dans l'ordinateur. Le plus grand entier est donc 2147483647. Il est possible de représenter un nombre 2 fois plus grand en considérant le type unsigned c'est-à-dire uniquement les entiers positifs. Pour des nombres encore plus grands, il faut utiliser le type long et unsigned long.

Mais qu'en est-il du type float ?

Ceci dépend ! La virgule est flottante et donc il y a une précision variable.

#include <stdio.h>
#include <math.h>
#include <limits.h>
int main () {
    float x=1., f=1., ff;
    while (!isinf(f)) {
        ff = f;

        f = f * x;

        x = x + 1.;
        printf("%.0f!  \t= %f\n", x-1, f);
    }
    return 0;
}
1!  	= 1.000000
2!  	= 2.000000
3!  	= 6.000000
4!  	= 24.000000
5!  	= 120.000000
6!  	= 720.000000
7!  	= 5040.000000
8!  	= 40320.000000
9!  	= 362880.000000
10!  	= 3628800.000000
11!  	= 39916800.000000
12!  	= 479001600.000000
13!  	= 6227020800.000000
14!  	= 87178289152.000000
15!  	= 1307674279936.000000
16!  	= 20922788478976.000000
17!  	= 355687414628352.000000
18!  	= 6402373530419200.000000
19!  	= 121645096004222976.000000
20!  	= 2432902023163674624.000000
21!  	= 51090940837169725440.000000
22!  	= 1124000724806013026304.000000
23!  	= 25852017444594485559296.000000
24!  	= 620448454699064672387072.000000
25!  	= 15511211079246240657965056.000000
26!  	= 403291499589617303175561216.000000
27!  	= 10888870415132690890901946368.000000
28!  	= 304888371623715344945254498304.000000
29!  	= 8841763079319199907069674127360.000000
30!  	= 265252889961724357982831874408448.000000
31!  	= 8222839685527520666638122083155968.000000
32!  	= 263130869936880661332419906660990976.000000
33!  	= 8683318509846655538309012935952826368.000000
34!  	= 295232822996533287161359432338880069632.000000
35!  	= inf

À comparer avec la séquence exacte:

1!      = 1
2!      = 2
3!      = 6
4!      = 24
5!      = 120
6!      = 720
7!      = 5040
8!      = 40320
9!      = 362880
10!     = 3628800
11!     = 39916800
12!     = 479001600
13!     = 6227020800
14!     = 87178291200
15!     = 1307674368000
16!     = 20922789888000
17!     = 355687428096000
18!     = 6402373705728000
19!     = 121645100408832000
20!     = 2432902008176640000
21!     = 51090942171709440000
22!     = 1124000727777607680000
23!     = 25852016738884976640000
24!     = 620448401733239439360000
25!     = 15511210043330985984000000
26!     = 403291461126605635584000000
27!     = 10888869450418352160768000000
28!     = 304888344611713860501504000000
29!     = 8841761993739701954543616000000
30!     = 265252859812191058636308480000000
31!     = 8222838654177922817725562880000000
32!     = 263130836933693530167218012160000000
33!     = 8683317618811886495518194401280000000
34!     = 295232799039604140847618609643520000000
35!     = 10333147966386144929666651337523200000000
36!     = 371993326789901217467999448150835200000000
37!     = 13763753091226345046315979581580902400000000

Donc nous voyons que la représentation en nombre flottant fait un calcul faux pour factorielle 14 (14!), et ensuite les erreurs se cumulent. Ceci est attendu et la précision dépend de ce que l'on cherche à calculer. Pour les calculs standards, nous sommes très rarement dans ce cas. Ici, nous testons les limites.

Par ailleurs, nous voyons que nous ne pouvons pas utiliser plus de 36 termes dans notre calcul de l'exponentielle utilisant la série.


3 Outils utilisés pour générer ce document

Le document est écrit avec GNU Emacs et Org-mode. À partir d'un unique fichier source example.org tout est généré avec la commande :

emacs -batch example.org -f org-babel-execute-buffer

ou en ouvrant le document example.org avec Emacs puis en pressant Control c Control v b.

Le dépôt contenant toutes les sources se trouve ici :

3.1 Erreur de coloration du code

Avec la configuration d'Emacs par défaut, vous risquez d'avoir ce message d'erreur:

Please install htmlize from https://github.com/hniksic/emacs-htmlize

dans ce cas, nous vous suggérons de lancer la commande---récuperation de la dépendance manquante:

git clone https://github.com/hniksic/emacs-htmlize.git /tmp/emacs-htmlize

et voilà. Sinon ouvrez et lisez le fichier example.org ou ouvrez un boggue.

3.2 Erreur d'exportation

Loading /home/simon/.guix-profile/share/emacs/site-lisp/gettext-autoloads.el (source)...
Symbol’s function definition is void: org-outline-overlay-data

Changements incompatibles dans la version 9.2. Voir ici.


© 2014-2022 Simon Tournier <simon (at) tournier.info>

(last update: 2022-04-19 Tue 21:45)