Applications mathématiques en Informatique

Sommaire

Introduction

Voici quelques notes concernant l'algorithmique au travers du langage Pascal, initialement écrites au fil de mes lectures et de mon apprentissage de ce langage. Elles sont rassemblées ici à titre personnel, et éventuellement à titre pédagogique, pour les personnes désireuses d'apprendre ce langage au travers d'exemples choisis, ou pour les personnes connaissant déjà le langage mais recherchant certains détails algorithmiques. Lorsque cela est possible (ou utile), la même version en langage C est fournie, à titre de comparaison. L'ensemble des sources figurant sur cette page ont été testés avec gcc (Gnu C Compiler, version 3.4.1) et fpc (Free Pascal Compiler, version 1.9.8), initialement sous GNU/Linux, puis sous MacOS X, avec gcc version 4.0.1 et fpc version 2.0.4-1.

Le calcul scientifique impose la plupart du temps de travailler avec des réels, et soulève par conséquent le problème (i) de la représentation de ces quantités en machine, et (ii) de la précision des résultats issus des calculs effectués. J'ai rassemblé ici quelques-unes des implications de ce type de calcul, au fur et à mesure de mon "expérimentation" de ce domaine.

Il y a également quelques algorithmes plus généraux que l'on trouve dans des cours de mathématiques pour l'informatique (ou vice-versa...) très ordinaires, mais qui constituent toujours un bon point de départ pour toute personne désireuse de débuter en programmation mathématique. De nombreuses notions présentées sur cette page m'ont été

On commencera par quelques rappels et considérations générales concernant la programmation mathématique, puis on présentera quelques-uns des algorithmes les plus connus, puis on terminera avec des fonctions plus avancées.

Les exemples proposés sont réalisés en Pascal la plupart du temps, quelques fois en C, et quelques implémentations en Scheme ou Maple sont également disponibles. La littérature étant vaste sur ce sujet, je donne ici à titre indicatif les ouvrages que j'ai consultés : [Bers1991], [Knut1997], [Boug1993]. Les autres références ponctuelles (pages internet ou ouvrages divers) sont indiquées directement dans le texte ou dans la bibliographie rassemblée à la fin de ce document.

L'ensemble des sources figurant sur cette page ont été testés avec gcc (Gnu C Compiler, version 3.4.1) et fpc (Free Pascal Compiler, version 2.0.0).

Calcul numérique

Calcul de sommes sans dépassement

Lorsque l'on calcule une somme de manière classique, c'est-à-dire par accumulation itérative des valeurs à sommer, il peut arriver qu'il se produise un dépassemment du résultat du point de vue de la capacité de représentation en machine, c'est-à-dire que la valeur de la somme temporaire calculée à un certain moment n'est pas représentable en machine, comme par exemple une somme codée en entier (integer) dépassant 32765. C'est le cas lorsque l'on travaille avec de grandes bases de données. Il est possible de remédier à cet inconvénient en utilisant l'idée suivante :

Voici le code correspondant

function somme(A : tab): integer;
var
   pos, neg: tab;
   npos, nneg, accu, i: integer;
begin
   npos := 0;
   nneg := 0;
   {-- étape 1 : tri en termes positifs et négatifs --}
   for i:=1 to N do
      if A[i] >= 0 then
      begin
         npos := npos + 1;
         pos[npos] := A[i];
      end
      else
      begin
         nneg := nneg + 1;
         neg[nneg] := A[i];
      end;
   {-- étape 2 : calcul de la somme --}
   accu := 0;
   while (npos > 0) and (nneg > 0) do
   begin
      if accu >= 0 then
      begin
         accu := accu + neg[nneg];
         nneg := nneg - 1;
      end
      else
      begin
         accu := accu + pos[npos];
         npos := npos - 1;
      end
   end;
   {-- étape 3 : ajout des derniers termes de même signe --}
   if npos > 0 then
      for i:=1 to npos do
         accu := accu + pos[i]
      else
         for i:=1 to nneg do
            accu := accu + neg[i];
   somme := accu;
end;

Une autre solution consiste à travailler avec des entiers "beaucoup" plus grands.

Calcul de la puissance d'un nombre réel

Il n'existe pas de fonction puissance sous Pascal, et plutôt que de répéter x * x * x * ... dans une instruction, il est parfois utile d'en implémenter une rapidement. L'avantage de la méthode récursive (cf. D. Knuth, The Art of Computer Programming, Addison Wesley, 2ème éd., 1981) est qu'elle requiert environ log2(n) multiplications, contrairement à la solution itérative classique qui nécessite n-1 multiplications.

Voici la fonction classique

function puissance(x : real; n : integer): real;
{méthode classique
 N.B.: traite les cas n<0, n=0 et n>0}
var
   i   : integer;

begin
   if n = 0 then
      puissance := 1.0
   else if n > 0 then
   begin
      puissance := x;
      for i:=2 to n do
         puissance := puissance * x
   end
   else
   begin
      puissance := 1/x;
      for i:=2 to abs(n) do
         puissance := puissance * 1/x
   end
end; { puissance }

et une variante

function puissance2(x : real; n : integer): real;
{méthode récursive de la chaîne chinoise
 N.B.: implémentée pour les puissances positives}
var
   moitie : real;

begin
   if n = 0 then
      puissance2 := 1.0
   else
   begin
      moitie := puissance2(x, n div 2);
      puissance2 := sqr(moitie);
      if (n mod 2 = 1) then
         puissance2 := puissance2 * x
   end;
end; { puissance2 }

Le programme puissance.pas permet de tester ces deux fonctions.

Estimation de pi par la méthode de Monte Carlo

Plusieurs méthodes d'estimation de pi existent et sont plus ou moins facilement implémentables en Pascal. La méthode de Monte Carlo est une méthode très simple de simulation qui consiste à

Voici le programme correspondant, pi.pas

program calcul_pi;
{ calcul de pi par la m-béthode de Monte Carlo }

var
   x, y             : real;
   d2       : real;
   pi       : real;
   np, nc, nr : integer; (* nb points tirés, nb points ok, nb répétitions *)
   i, j             : integer;

function aleat: real;
begin
   aleat := random(32766)/32767;
end; { aleat }

{ -- bloc principal --}
begin
   write('Combien de points ? ');
   readln(np);
   write('Combien de répétitions ? ');
   readln(nr);
   randomize;
   for i:=1 to nr do
   begin
      pi := 0.0;
      nc := 0;
      for j:=1 to np do
      begin
         x := aleat;
         y := aleat;
         d2 := (x-0.5)*(x-0.5) + (y-0.5)*(y-0.5);
         if d2 <= 0.25 then
            nc := nc + 1;
      end;
      pi := (4.0*nc)/np;
      writeln('estimation de pi avec ', np, ' points : ', pi);
   end;
end.

Echange de valeurs

L'échange du contenu de deux variables (a prend la valeur de b et réciproquement) peut se faire de plusieurs manières, et dépend du type de langage utilisé. Le plus souvent, on utilise une variable auxiliaire (avec passage par valeur, ou en définissant une fonction d'échange avec passage des paramètres par adresse [1]), mais il est possible d'échanger la valeur de 2 variables sans utiliser une telle variable auxiliaire. Par exemple, les instructions suivantes permettent d'échanger le contenu des variables a et b:

x := x + y;
y := x - y;
x := x - y;

Cependant, cette méthode n'est pas une bonne façon de procéder à la permutation des valeurs de 2 variables, notamment lorsque les variables sont de type réel. En effet, si a = 1E+8 et b=1E-6, a et b auront la même valeur après l'échange. Pourquoi ? Parce que l'écart entre ces deux valeurs dépassent la précision en machine et 1E+8+1E-6 = 1E+8.

Comparaison de deux réels

On évitera de tester directement l'égalité entre deux réels, du type if (a == b), dans la mesure où ceux-ci peuvent différer à cause de la précision machine. On trouve souvent le test suivant (EPS désignant une constante très petite, e.g. 1E-6)

if (fabs (x - y) < EPS

Knuth [Knut1997] recommande cependant de privilégier

if (fabs (x - y) < fabs(x)*EPS

puisque cela permet de s'affranchir du problème lié à une éventuelle différence entre x et y qui dépasserait la précision machine. En effet, si (...) (texte à revoir) [voir également http://www.faqs.org/faqs/C-faq/faq/, sec. 14.5]

Résolution d'équation du 2ème degré

Les solutions classiques d'un polynôme du deuxième degré en x, ax²+bx+c=0, (en supposant a différent de 0, sinon on se ramène à un problème du premier degré), sont obtenues classiquement comme suit :

  • si delta > 0, il existe 2 solutions distinctes :

    x = (-b +/- sqrt(delta))/2a

  • si delta = 0, il existe une solution double :

    x = -b/2a

  • si delta < 0, il n'y a pas de solution réelle.

Une des conséquences de la représentation des réels en machine est que l'addition ou la soustraction de valeurs très proches (ici, -b et sqrt(delta)) peut amener des résultats totalement inexacts.

Voici un exemple en Pascal (trinome.pas)

program trinome;

const EPS = 1E-10;

var
   a, b, c : real;
   delta   : real;
   x1, x2  : real;

function sign(D : real): integer;
begin
   if D > 0.0 then
      sign := +1
   else if D < 0.0 then
      sign := -1
   else
      sign := 0
end; {sign}

begin
   write('a = ');
   readln(a);
   write('b = ');
   readln(b);
   write('c = ');
   readln(c);
   delta := b*b - 4*a*c;
   if delta >= 0 then
   begin
      x1 := (-b - sign(b) * sqr(delta))/(2.0*a);
      x2 := c/(a*x1);
      if delta < EPS then
         writeln('la racine double est : x = ', x1:8:4)
      else
         writeln('les racines réelles sont x1 = ', x1:8:4, ' et x2 = ', x2:8:4)
   end
   else
      writeln('pas de racines réelles.')
end.

On pourrait bien évidemment raffiner le programme pour calculer également les racines complexes, selon le même schéma (voir trinome2.c).

Suites et séries numériques

Calcul de racine

On peut approximer la valeur de la racine de a à l'aide d'une simple suite numérique, définie telle que:

u0 = 10, u(n+1)=un+a/un+1.

En Pascal, cela donne (racine.pas):

program suite1;

const Precision = 1E-9;

var
   n    : integer;
   a    : integer;
   u, v : real;

begin
   n := 1;
   write('a = ');
   readln(a);
   write('u0 = ');
   readln(u);
   repeat
      v := u;
      u := (u+a)/(u+1);
      writeln('u', n, ' = ', u:16:15);
      n := n+1
   until abs(u-v) < Precision
end.

Rien de bien extraordinaire là, mais cela permet de construire, sur le même principe, des calculs sur des suites beaucoup plus complexes, Notons au passage que l'on n'a pas construit une fonction suite qui renverrait le terme demandé, mais qu'on a utilisé une solution itérative. Dans de nombreux cas, cela est suffisant, mais on peut vouloir implémenter des solutions récursives. Les suites se prêtent en effet par excellence à l'utilisation de la récurrence, et on exploitera au mieux la notion de récursivité, comme dans le calcul des termes de la fonction d'Ackerman (cf. infra).

Suite de Fibonacci

Pour la petite histoire, la suite de Fibonacci modélise la croissance des lapins

Voici une solution récursive:

function fibo_rec(n:integer):longint; {Les résultats obtenus sont grands!}
begin
   if n=0 then
      fibo_rec:=1
   else if n=1 then
      fibo_rec:=1
   else
      fibo_rec:=fibo_rec(n-1)+fibo_rec(n-2)
end;

La solution itérative, exposée ci-dessous, est préférable à la procédure récursive car cette dernière entraîne l'évaluation répétée des mêmes expressions : par exemple, pour obtenir F(5), on calcule F(4) et F(3), sachant que F(4) est calculée à partir de F(3) et F(2) -- on évalue donc deux fois F(3), pour le seul calcul de F(5)... ce qui à terme risque de poser de sérieux problèmes en termes de performances.

fibo_rec(5)

Fig. 1 - Illustration des étapes de calcul de fibo_rec(5) (Tiré de [b], fig 1.5)

On peut comparer les temps d'exécution des 2 fonctions afin de s'apercevoir que la solution itérative est de loin la meilleure, dès que n dépasse 30 (j'utilise pour ma part un pentium M à 1.7 GHz). Bien évidemment, l'évaluation de la fonction récursive pour n=500 est hors de propos (en Pascal ou en C). Voici une solution itérative (fibonacci.pas):

function fibo_iter(n:integer):longint;
var
   tab : array[0..MAX] of longint;
   i   : integer;
begin
   tab[0]:=1;
   tab[1]:=1;
   for i:=2 to n do
   begin
      tab[i]:=tab[i-1]+tab[i-2]
   end;
   fibo_iter:=tab[n];
end;
Remarque:

A propos de la solution récursive, lors du calcul de la factorielle d'un nombre, Maple utilise une astuce (option remember) qui permet de conserver les valeurs intermédiaires : ainsi, pour calculer 5000!, on pourra calculer successivement 1000!, 2000!, 3000!, 4000! et 5000! à l'aide d'une fonction du type:

fact:=proc(n) option remember;
if n=1 then 1; else n*fact(n-1) fi;
end;

Le calcul par dichotomie -- selon que n est pair ou impair -- permet d'améliorer sensiblement la rapidité de l'évaluation, lorsque n est grand ; les formules de calcul sont les suivantes

F(2p)= (2×F(p+1)-F(p))×F(p), F(2p+1)= F^2(p+1)+F^2(p)

En utilisant bc, cela donne [a]:

scale=0;
define fib(n) {
  auto a, b;
  if(n<3) {
    if(n==0) return 0;
    return 1;
  }
  a = fib(n/2); b = fib(n/2+1)
  if(n%2) return b^2 + a^2;
  return (2*b-a)*a ;
}
fib(3000);
quit;

On peut également utiliser la propriété suivante : pour 0<= k<= n-2, F(n)=F(k+2)*F(n-k-1)+F(k+1)*F(n-k-2) ; cela permet de ne calculer qu'une petite partie des nombres de Fibonacci précédant F(n) (le gain de temps est d'environ 125 %). Par exemple, en prenant k+1=50, on ne calcule que 2 nombres de Fibonacci consécutifs toutes les 50 positions.

On peut même évaluer cette fonction pour n=3000 en Scheme [a]

1 ]=> (define (F n) (Fs n 1 0))
(define (Fs n s c)
   (if (= n 0)
     c
     (Fs (- n 1) (+ s c) s)))

;Value: f

1 ]=>
;Value: fs

1 ]=> (F 3000)

qui donne comme résultat

41061588630797126033356837871926710522012510863736925240888543092690
55842741134037313304916608500445608300368357069422745885693621454765
02674373045446852160486606292497360503469773453733196887405847255290
08204908690751262205905454219588975803110922267084927479385953913331
83712447955431476110732762400667379340851917318109932017067768389347
66764778739502174470268627820918553842225858306408301661862900358266
85723821023580250435195147299791967652400478423637645334726836415264
83462458405732142414199379172429186026398100978669423920154046201538
18671425739835074851396421139982713640679581178458198658692285968043
243656709796000
[a](1, 2) http://perso.wanadoo.fr/jean-paul.davalan/divers/fibonacci/f02.html#PROG
[b](1, 2) http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html#%_sec_1.2.1

Fonction d'Ackerman

On appelle fonction d'Ackerman la fonction définie de N² dans N par :

F(x,y) = y+1, si x=0

= F(x-1,1), si y=0

= F(x-1,F(x,y-1)), si x<>0 et y<>0

Cette fonction est en apparence banale, mais peut se révèler très complexe lorsqu'il s'agit d'énumérer les étapes intermédiaires de calcul (à la différence de l'algorithme récursif de la suite de Fibonacci, pour lequel il s'agit simplement de répéter ou réécrire des calculs déjà effectués au préalable). Par exemple, pour calculer F(0,2), la réponse est immédiate : F(0,2)=3. En fait, dès qu'on tombe dans le cas x=0, c'est rapide ! Maintenant, calculons juste F(2,0) :

F(2,0)=F(1,1)=F(0,F(1,0))=F(0,F(0,1))=F(0,2)=3

ou (autre méthode) :

F(2,0)=F(1,1)=F(0,F(1,0))=F(1,0)+1=F(0,1)+1=3

En fait, on peut montrer que pour tout y (entier naturel) :

F(0,y) = y+1

F(1,y) = y+2

F(2,y) = 2y+3

F(3,y) = 8*(2y-3)

Mais, pour x > 3, il n'existe plus de procédé systématique de calcul. Lorsque l'on songe que F(3,10)=8189, il est évident que son calcul manuel relève du défi...

La fonction d'Ackerman peut donc être évaluée comme suit (ackerman.pas)

function func(x, y : integer) : integer;
begin
   if x = 0 then
      func := y+1
   else
      if y = 0 then
         func := func(x-1, 1)
      else
         func := func(x-1, func(x, y-1))
end; { func }

On notera qu'il existe d'autres définitions de la fonction d'ACkerman (cf. http://perso.wanadoo.fr/jean-paul.davalan/mots/suites/ack/), mais qu'elles sont toutes basées sur la fonction définie ci-dessus, qui est celle proposée par Ackerman.

Dans l'ouvrage de Abelson et al. [b], on trouve une procédure de calcul de la fonction d'Ackerman, écrite en Scheme

(define (A x y)
  (cond ((= y 0) 0)
        ((= x 0) (* 2 y))
        ((= y 1) 2)
        (else (A (- x 1)
                 (A x (- y 1))))))

On voit ici la concision avec laquelle ce langage permet de représenter des fonctions mathématiques. C'est ici bien évidemment une solution récursive, puisque le langage fonctionnel Scheme repose exclusivement sur ce mode d'évaluation, mais...

Arithmétique

Calcul de pgcd

Le type integer ne convient pas à la manipulation de très grands entiers, en raison de la capacité de représentation de tels nombres (limitée à 32765). La représentation des entiers en machine peut se faire de différentes manières :

  • utiliser un type constitué de tableaux de nombres de type integer classique mais écrits en base 10 (ou 100)
  • modifier le type réel (e.g. real en Pascal) existant

La première méthode (dite de calcul en précision multiple) impose de recréer également l'ensemble des opérations arithmétiques courantes : addition, multiplication et division euclidienne. La deuxième méthode est plus rapide, mais exige de limiter l'étendue des nombres utilisés, de façon à obtenir des calculs exacts.

Notons que l'exercice est un exercice de forme car la plupart des compilateurs autorise l'usage de type prédéfini pour les entiers (e.g. ??). D'autre part, en Pascal, le type longint permet également de travailler avec des nombres entiers plus grands.

Pour calculer le pgcd de 2 nombres entiers, on a besoin de 2 fonctions arithmétiques spéciales : le calcul de la partie entière et la division euclidienne (avec le calcul du reste entier). On suppose qu'on a défini le type entier à partir du type prédéfini réel, e.g. type entier = real;. Une solution possible pour ces trois fonctions est la suivante

function ent(x: entier): entier;
var s : entier;

begin
   s := int(x);
   if (x < 0) and (s <> x) then
      s := s-1;
   ent := s
end; { ent }

function div_e(x, y : entier) : entier;
var q : entier;

begin
   q := x/y;
   div_e := ent(q)
end; { div_e }

function mod_e(x, y : entier) : entier;
begin
   mod_e := x - div_e(x, y)*y
end; { mod_e }

Cela permet d'écrire une fonction de calcul du pgcd (pgcd.pas) sous la forme

function pgcd(x, y: entier) : entier;
var r : entier;

begin
   x := abs(x);
   y := abs(y);
   if x < y then
   begin
      r := x;
      x := y;
      y := r
   end;
   while y <> 0 do
   begin
      r := mod_e(x, y);
      x := y;
      y := r
   end;
   pgcd := x
end; { pgcd }

Combinatoire

Permutations

On peut tester le programme à l'aide de test_perms.pas, en donnant comme série {1,3,5,4,2}. Le programme indique que le rang de cette permutation de 5 éléments est 12.

On peut vérifier que le rang indiqué (12) correspond à celui que donne Mathematica

ln[5]:= Permutations[{1,2,3,4,5}]
Out[5]= {{1, 2, 3, 4, 5}, {1, 2, 3, 5, 4}, {1, 2, 4, 3, 5}, {1, 2, 4, 5, 3},

>    {1, 2, 5, 3, 4}, {1, 2, 5, 4, 3}, {1, 3, 2, 4, 5}, {1, 3, 2, 5, 4},

>    {1, 3, 4, 2, 5}, {1, 3, 4, 5, 2}, {1, 3, 5, 2, 4}, {1, 3, 5, 4, 2},
                                                        ^^^^^^^^^^^^^^^
>    {1, 4, 2, 3, 5}, {1, 4, 2, 5, 3}, {1, 4, 3, 2, 5}, {1, 4, 3, 5, 2},

>    {1, 4, 5, 2, 3}, {1, 4, 5, 3, 2}, {1, 5, 2, 3, 4}, {1, 5, 2, 4, 3},

> ...

Notons que la bibliothèque STL de C++ définit des fonctions pour travailler avec les permutations, dans stl_algo.h (généralement dans /usr/include/c++/4.0/bits). Il suffit d'inclure la bibliothèque algorithm pour pouvoir utiliser ces fonctions

$ grep permutation stl_algo.h
// next_permutation and prev_permutation, with and without an explicitly
 *  @return  False if wrapped to first permutation, true otherwise.
 *  Treats all permutations of the range as a set of "dictionary" sorted
  next_permutation(_BidirectionalIterator __first,
 *  @return  False if wrapped to first permutation, true otherwise.
 *  Treats all permutations of the range [first,last) as a set of
  next_permutation(_BidirectionalIterator __first,
 *  @return  False if wrapped to last permutation, true otherwise.
 *  Treats all permutations of the range as a set of "dictionary" sorted
  prev_permutation(_BidirectionalIterator __first,
 *  @return  False if wrapped to last permutation, true otherwise.
 *  Treats all permutations of the range [first,last) as a set of
  prev_permutation(_BidirectionalIterator __first,

Le programme permutation.cc permet de tester très brièvement la génération de permutations en utilisant C++.

Calcul matriciel

Références

[Bers1991]Berstel, J., Pin, J.-E. et Pocchiola, M. (1991). Mathématiques et Informatique, 1. Algèbre. McGraw-Hill.
[Knut1997](1, 2) Knuth, D.E. (1997). The Art of Computer Programming, Volume 1: Fundamental Algorithms. Addison-Wesley.
[Boug1993]Bougé, L., Kenyon, C., Muller, J.-M., Robert, Y. (1993). Algorithmique - Exercices corrigés, Oral du concours d'entrée à l'Ecole Normale Supérieure de Lyon. (Ellipses)
[Dela1996]Delannoy, C. (1996). Exercices en Turbo-Pascal. (Eyrolles)
[web]http://www.cs.sunysb.edu/~algorith/ ; http://www2.toki.or.id/book/AlgDesignManual/

Notes

[1]

ou encore en définissant une macro, par exemple en C:

#define echange(a,b)