Algorithmes pour l’arithmétique des entiers

Le but est d’initier à l’écriture d’algorithmes récursifs terminaux, sur des objets mathématiquement simples. Les bouts de codes en python sont rassemblés dans le fichier arith_int.py.

Factorielles et récursivité

définition

On note n! (lire factorielle n) le produit des nombres entiers de 1 à n : n! = \prod_{j=1}^{n}j.

Quand n vaut 0, la liste des facteurs est vide, le résultat doit être l’élément neutre de la multiplication : 0!=1.

Le code python qui correspond à cette définition est :

def fact1(n):
  f = 1
  for j in range(1,n+1):
    f = f*j
  return f

Mais en python, on peut plus directement utiliser reduce :

def fact1bis(n):
  return reduce(n.__mul__,range(1,n+1),1)

Il y a aussi une définition récursive de n! : le cas de base est 0!=1, et le cas de récursion est (n+1)!=(n+1) \times n!, ce qui conduit au code :

def fact2a(n):
  if n==0:
    return 1
  else:
    return n*fact2a(n-1)

Là encore python dispose d’une syntaxe souple et riche pour les expressions conditionnelles :

def fact2abis(n):
  return 1 if n==1 else n*fact2abis(n-1)

Le problème avec cette version récursive est la limitation de la taille de la pile qui permet à python de gérer les appels récursifs.

On peut chercher une version récursive terminale :

def fact2b(n,f=1):
    # return f if n==1 else fact3(n-1,n*f)
    return f if n==0 else fact2b(n-1, n*f)

(noter l’argument optionnel f qui sert d’accumulateur et prend par défaut la valeur 1) puis la dérécursifier, ce qui donne peu ou prou la première version :

def fact2(n, f=1):
    while n>0:
        n, f = n-1, n*f
    else:
        return f

On a gardé l’argument optionnel : le résultat de fact2(n,f) est n! \times f. C’est aussi l’invariant de la boucle while.

Pour s’entainer, généralisons aux nombres d’arrangements.

définition

On note n^{\underline{p}} (lire n puissance p sousligné) le produit dégressif de p facteurs commençant à n : n^{\underline{p}} = \prod_{j=0}^{p-1} (n-j) = n(n-1)\dots(n-p+1).

La définition récursive est limpide : si p vaut 1, le résultat est n ; sinon, le résultat est n fois le résultat obtenu pour n-1 et p-1 à la place de n et p :

n^{\underline{p+1}} &= n (n-1)^{\underline{p}}, \qquad [p\in\mathbb{N}] \\ n^{\underline{0}} &= 1

Il en découle que n^{\underline{n}} = n! et que n^{\underline{p}} = 0 dès que p>n. On obtient le code :

def puis_sous(n,p):
  if p==0:
    return 1
  elif p==1:
    return n
  else:
    return n*puis_sous(n-1,p-1)

C’est bien mais c’est récursif, pas terminal.

exercice

Dérécursifier la fonction puis_sous comme on l’a fait pour fact2.

Puissances et diviser pour régner

On veut maintenant calculer a^n, qui est le produit de n facteurs égaux à a. Bien sûr on peut utiliser la définition récursive évidente : si n vaut 0, le résultat est 1, et sinon le résultat est a \times a^{n-1}. On a alors n-1 multiplications à effectuer (car on ne compte évidemment pas la première, de la forme 1 \times a).

On peut décomposer le problème en deux : si n=p+q, on calcule a^p \times a^q. Par exemple si n=2p le résultat est (a^p)^2 ou encore (a^2)^p. Si on note C(n) le nombre de multiplications éffectuées pour élever à la puissance n le nombre a, cette méthode donne C(p+q) = 1 + C(p)+C(q) si p \neq q, mais seulement C(2p) = 1 + C(p) quand p = q. Pour les petites valeurs on obtient :

n 1 2 4 8
C(n) 0 1 2 3

On devine, et on peut démontrer par récurrence, que C(2^j)=j, ie C(n) = \ln_2 n pour n une puissance entière de deux.

Quand n est impair, il y a une mutliplication par a de plus : C(2p+1) = 2 + C(p). En notant n \div 2 et n \mod 2 le quotient et le reste dans la division euclidienne de n par 2, la formule générale est C(n) = 1 + n \mod 2 + C(n \div 2).

Voici trois versions de l’algorithme correspondant. Les bits de n (ie les chiffres de son écriture en base 2) sont lus et utilisés dans l’ordre des poids croissants : d’abord 1=2^0, puis 2=2^1, puis 4= 2^2, puis 8=2^3, etc. C’est en rapport avec l’écriture sous forme de polynôme, ou de série :

n = \eps_0 2^0 + \eps_1 2^1 + \eps_3 2^3 + \dots +\eps_j^j

Nous avons décidé de ne pas prendre en compte dans la suite le cas trivial a^0 = 1. Les fonctions suivantes sont correctes quand leur argument n est un entier naturel non nul.

  • version récursive :
def puis_rap1a(a,n):
    if n==1:
        return a
    else:
        b = puis_rap1a(a,n/2)
        return (1 if n%2==0 else a)*b*b
  • version récursive terminale :
def puis_rap1b(a,n,b=1):
    return a*b if n==1 else puis_rap1b(a*a, n/2, b if n%2==0 else a*b)
  • version dérécursifié :
def puis_rap1(a,n,b=1):
  while n>1:
    a, n, b = a*a, n/2, b if n%2==0 else a*b
  else:
    return a*b
  

exercice

Démontrer que dans cet algotihtme \lfloor\ln_2 n \leq C(n) \leq 2\lfloor\ln_2 n. On a noté \lfloor x la partie entière d’un réel x, c’est-à-dire le plus grand entier j tel que j\leq x.

Les multiplications que l’on effectue dans cet algorithme concernent des nombres dont la taille augmente au cours du calcul.

Voici maintenant un algorithme où les multiplications sont soit des mises au carré soit des multiplication par a. Il est relié au shéma de Hörner pour l’évaluation de l’écriture polynômiale de n :

n = (\dots((\eps_j \times 2 + \eps_{j-1})\times 2 + \eps_{j-2})\times 2 + \dots \eps_1)\times 2 + \eps_0

En revanche, puisqu’il utilise les bits de n dans l’ordre des poids décroissants, il faut commencer par décomposer n en base 2.

def puis_rap2(a,n,b=1):
  l=[]
  while n>1:
    l.append(n%2)
    n=n/2
  else:
    x=a
    while l<>[]:
      x=x*x
      if l.pop()==1:
        x=x*a
    else:
      return x*b

Division euclidienne

théorème

Étant donné un couple (a,b) d’entiers relatifs, b> 0, il existe un unique couple (q,r) tel que a = bq + r et 0 \leq r < b.

définition

On dit que q est le quotient et r le reste dans la division de a par b.

démonstration

Appelons amissibles les couples (x,y) vérifiants a = bx + y. Alors (0,a) est admissible. Et, si (x,y) l’est alors f_+(x,y) = (x+1,y-b) et f_-(x,y) = (x-1,y+b) le sont aussi. On cherche un couple (x,y) admissible tel que les deux conditions C_1 : 0 \leq y et C_2 : y < b soient vraies. On a toujours au moins une des deux. On part de (x,y) = (0,a). Si C_1 est fausse, on applique f_-, ce qui augmente y qui était trop petit ; si C_2 est fausse, on applique f_+, ce qui diminue y qui était trop grand. On répète jusqu’à ce que la condition fausse devienne vraie, ce qui est obligé d’arriver car l’ordre usuel sur \mathbb{N} est bien fondé : il n’existe pas de suite infinie strictement décroissante d’entiers naturels. Et alors les deux conditions sont vérifiées.

L’unicité revient à dire, en faisant la différence de deux solutions, que parmi les y vérifiant 0 \leq y < b, le seul multiple de b est 0.

L’algorithme évident contenu dans la démonstration demande q soustraction. En voici une version qui est correcte si l’argument a fourni est un entier naturel et l’argument b un entier naturel non nul.

def quo_res1(a,b):
  q,r=0,a
  while r>=b:
    q,r=q+1,r-b
  else:
    return q,r

De même que pour les puissances, on peut chercher un algorithme plus efficace, qui obtient les bits du quotient q dans l’ordre des poids décroissants. La première boucle détermine en montant le poid du bit de poid maximal. La boucle suivante redescend et met à jour le quotient et le reste.

def quo_res2(a,b):
  k,c=0,b
  while c<=a:
    k,c=k+1,c*2
  q,r,c=0,a,c/2
  for j in range(k):
    if r<c:
      q,c=q*2,c/2
    else:
      q,r,c=q*2+1,r-c,c/2
  return q,r

exercice

Pour les deux fonctions quo_res1 et quo_res2, donner pour chaque boucle les propriétés invariantes et justifier la terminaison.

Dans cet algorithme, les opérations binaires k=k+1, c=c*2, c=c/2, q=q*2 et q=q*2+1 sont considérées comme élémentaires et négligeables. Les opérations importantes sont les soustractions r=r-c. Il y en a au plus k, c’est-à-dire le nombre de bits dans l’écriture de q en base 2 : \lfloor\ln_2 q.

Algorithme d’Euclide

définition

Le plus grand commun diviseur (en abrégé pgcd) de deux entiers relatifs a et b est l’entier naturel d le plus grand parmi ceux qui divisent à la fois a et b. On le note a \wedge b.

Deux entiers relatifs sont premiers entre eux quand leur pgcd est 1.

Le calcul du pgcd de deux entiers est une opération courante en arithmétique, par exemple pour simplifier ou additionner des fractions. On dispose pour cela de l’algorithme d’Euclide qui repose sur les constatations suivantes :

  1. Le pgcd de deux entiers relatifs ne change pas si l’on échange les deux nombres ni si l’on remplace l’un des nombre par son opposé.
  2. Le pgcd de deux entiers naturels distincts ne change pas si l’on remplace le plus grand des deux par la différence des deux.
  3. Le pgcd de deux entiers naturels égaux est égal à ces nombres. Ou encore : si l’un des nombres est nul, le pgcd est égal à l’autre.

exercice

Écrire une fonction qui retourne le pgcd de ses deux arguments supposés entiers naturels, en suivant l’algorithme décrit ci-dessus.

On peut accélerer l’algorithme en remarquant que si l’étape 2 est répétée plusieurs fois sans échanger les deux nombres, c’est que l’on est en train de calculer le reste dans la division euclidienne du plus grand par le plus petit. Cela donne, en utilisant l’opérateur standard a%b pour le reste de a divisé par b (mais on pourrait tout aussi bien utiliser nos fonctions quo_res1 ou quo_res2 :

def pgcd1(x,y):
  while y<>0:
    x,y=y,x%y
  return x

Là encore on pourrait chercher une version binaire. Pour cela on commence éliminer le cas où l’un des nombres est nul. Puis on extrait la plus grande puissance de 2 commune, puis on se ramène au cas de deux nombres impairs :

def pgcd2(x,y):
  if x==y or y==0:
    return x
  elif x==0:
    return y
  c=1
  while x%2==0 and y%2==0:
    x,y,c=x/2,y/2,c*2
  if y%2==0:
    x,y=y/2,x
  while True:
    while x%2==0:
      x=x/2
    if x>y:
      x=(x-y)/2
    elif x<y:
      x,y=(y-x)/2,x
    else:
      return c*x

théorème

Soit a et b deux entiers relatifs.

  1. Quels que soient les entiers relatifs x et y, le nombre ax+by est un mutliple de a \wedge b.
  2. Tout mutliple de a \wedge b peut s’écrire sous la forme ax+by avec x et y entiers relatifs bien choisis.

définition

En particulier, il existe deux entiers relatifs u et v tels que au + bv = a \wedge b. Ceux sont des coefficients de Bezout des nombres a,b.

La partie 2 du théorème et la définition des coefficients de Bezout se généralise en l’étude de l’équation ax+by=c. Les coefficients a,b,c sont des entiers relatifs. Une solution de l’équation est un couple (x,y) d’entiers relatifs qui vérifient l’équation.

Si c n’est pas un multiple de a\wedge b, il n’y a pas de solutions. Si c=a\wedge b\,c', alors l’ensemble des solutions n’est pas vide. Il y a existence mais pas unicité : si (x_0,y_0) est solution, alors pour tout j\in \Z, le couple (x_0-jb',y_0+ja') = (x_0,y_0) + j (-b',a') est encore solution, où a',b' vérifient (a,b)=a\wedge b\,(a',b') (et donc aussi a'\wedge b'=1, ceux sont même les plus grands diviseurs respectivement de a et b premiers entre eux).

C’est le même principe de décomposition en solution particulière de l’équation complète (avec second membre c) plus solution générale de l’équation homogène (ax+by=0) que pour les systèmes linéaires et les équations différentielles linéaires.

Pour trouver une solution particulière de l’équation complète, on se ramène à a'x+b'y=c' avec maintenant a'\wedge b' = 1. Si on connaît des coefficients de Bezout u',v' pour les nombres premiers entre eux a',b', on peut choisir comme solution particulière c'(u',v').