Algorithmes arithmétiques

Introduction

Au cours de cet article, nous allons nous implémenter quelques fonctions arithmétiques élémentaires sous Maple : nombre de diviseurs, test de primalité, etc. On s'intéressera dans un premier temps à des solutions naïves, puis on tentera d'aboutir à de meilleures complexités, notamment pour le calcul de la somme des diviseurs où, à l'aide de quelques résultats mathématiques, on optimisera sérieusement les performances de la solution.

Initialisation

restart;
R := 42:

L'entier que nous venons de définir servira de bornes aux tableaux statiques que nous allons allouer. En pratique, on déclarerait des tableaux dynamiquement et non de façon statique, mais une allocation dynamique ne présente aucun intérêt ici. Pour cette raison, on n'utilisera pas R dans les déterminations de complexités car il s'agit d'une constante fictive, mais plutôt n qui est la composante variable à l'appel des procédures.

Approche instinctive

On peut envisager une approche simple du problème qui consiste à tester tous les candidats possibles pour déterminer les diviseurs de l'entier n.

Nombre de diviseurs

Déterminer le nombre de diviseurs d'un entier n revient alors à considérer tous les entiers naturels i entre 1 et n et à tester pour chacun la valeur du reste dans la division euclidienne de n par i : si celle-ci est nulle, i divise n et on peut incrémenter le compteur résultat. Considérons l'implémentation suivante de cette idée :

nbrediv := proc(n)
    local i, resultat;
    resultat := 0;
    if (n = 0) then
        resultat := infinity;
    fi;
    for i from 1 to n do
        if (irem(n, i) = 0) then
            resultat := (resultat + 1);
        fi;
    od;
    resultat;
end:

Sa complexité temporelle est en O(n), c'est à dire que le le nombre d'opérations élémentaires que nécessite la procédure est proportionnel à la valeur de l'entier n passé en paramètre (ici, on a donc n multiplications en tout). Cette solution n'est envisageable jusqu'à un certain seuil, par exemple :

$\large n\,\leq\,10^8$

Valeur pour laquelle il faudrait déjà plus de 10 secondes à un ordinateur courant pour calculer le résultat... Heureusement pour l'humanité, il existe des solutions plus performantes !

Test de primalité

Un entier naturel n est premier s'il admet deux diviseurs : 1 et lui-même. Une procédure triviale se déduit alors de la précédente :

premier := proc(n)
    evalb(nbrediv(n) = 2)
end:

Sa complexité est toujours en O(n) et on en tire les mêmes conclusions que précédemment. Cependant, sans aller chercher trop loin, on sait qu'un diviseur premier, s'il existe, est plus petit que la racine carrée de l'entier considéré, ce qui permet de reprendre l'algorithme précédent en y ajoutant une petite amélioration :

premier := proc(n)
    local i, estPremier;
    estPremier := true;
    if (n = 0 or n = 1) then
        estPremier := false;
    fi;
    for i from 2 to isqrt(n) do
        if (irem(n, i) = 0) then
            estPremier := false;
        fi;
    od;
    estPremier;
end:

Ceci amène à une complexité en $\small O(\sqrt{n})$ plus performante que la précédente, car le calcul peut désormais se faire pour :

$\large n\,\leq\,10^{16}$

Toujours pour un processeur aux alentours du Gigahertz, comme c'est le cas des ordinateurs personnels à notre époque.

Test de parfaite quadrature

Suivant la même approche exhaustive, on peut tester si un entier naturel n est de la forme :

$\large n\,=\,s^2$

En évaluent tous les candidats possibles entre 0 et la partie entière de la racine carrée de n (il est inutile de chercher des candidats après). On peut le faire comme suit :

etrecarre := proc(n)
    local i, estCarre;
    estCarre := false;
    for i from 0 to isqrt(n) do
        if ((i^2) = n) then
            estCarre := true;
        fi;
    od;
    estCarre;
end:

Auquel cas la complexité finale est en $\small O(\sqrt{n})$ .

Vers une solution plus efficace

Il est possible d'aboutir à une solution bien plus performante en faisant appel au crible d'Ératosthène pour obtenir la décomposition en facteurs premiers d'un entier n donné. On utilisera la commande Maple appropriée pour obtenir le résultat avant de s'intéresser à l'implémentation propre du crible.

Tableau des exposants

Pour obtenir la décomposition en facteurs premiers d'un entier, on peut utiliser la fonction ifactor du logiciel (on aurait aussi pu utiliser ifactors) en supposant connue la liste des nombres premiers diviseurs potentiels de cet entier. Déclarons en premier lieu les tableaux que nous allons utiliser :

alpha := array(1..R):
premiers := array(1..R):

Initialisons le tableau des exposants comme suit :

for p from 1 to R do
    alpha[p] := 0;
od:

Nous verrons tout à l'heure comment construire le tableau des premiers nombres premiers. Pour l'instant, on peut se contenter d'utiliser une fonction de substitut :

for p from 1 to R do
    premiers[p] := ithprime(p);
od:

Ceci fait, on récupère la décomposition de l'entier comme suit :

rho := 42:
delta := ifactor(rho);

Et il ne reste plus qu'à traiter son arbre interne pour remplir notre tableau :

p := 1:
for i from 1 to nops(delta) do
    if (nops(op(i, delta)) = 1) then
        valeur := op(1, op(i, delta));
        exposant := 1;
    else
        valeur := op(1, op(1, op(i, delta)));
        exposant := op(2, op(i, delta));
    fi;
    while (premiers[p] < valeur) do
        p := (p + 1);
    od;
    alpha[p] := exposant;
end:

Voilà pour une première approche. Nous verrons ci-après le crible qui évite l'appel à cette fonction ithprime dont on ne connait pas la complexité.

Nombre de diviseurs

Il est désormais aisé de calculer le nombre de diviseurs d'un entier n à partir du tableau construit. On sait que :

$\large d(n)\,=\,\prod_{i=1}^{r} (\alpha_i\,+\,1)$

D'où directement :

nbrediv2 := proc(n)
    local i, resultat;
    resultat := 1;
    for i from 1 to isqrt(n) do
        resultat := resultat * (alpha[i] + 1);
    od;
end:

On vérifie bien que :

evalb(nbrediv(rho) = nbrediv2(rho));

renvoie true ; mais la complexité de cette nouvelle mouture est en $\small O(\sqrt{n})$ ce qui est bien plus performant que du O(n) comme l'illustrent ces quelques temps d'exécution, en secondes sur une machine à 2.0 Ghz :

time(nbrediv(1000000));

3.295

time(nbrediv2(1000000));

0.009

Test de primalité

On utilisant à nouveau cette représentation, le test de primalité se fait aussi en $\small O(\sqrt{n})$ , mais on peut donner dans la concision :

premier2 := proc(n)
    evalb(nbrediv2(n) = 2)
end:

Il est néanmoins possible de faire d'une pierre n coups en implémentant un crible qui nous permettra de tester la primalité en O(1) (à terme) et de remplir le tableau des premiers nombres premiers sans passer par la fonction ithprime du logiciel.

Crible d'Eratosthène

Le but de ce crible est de déterminer les nombres premiers plus petits qu'un certain entier. Mais ici, nous allons détourner légèrement l'idée initiale afin de stocker non-pas un booléen sur la primalité d'un entier, mais plutôt le nombre de diviseurs de l'entier en question. Le principe de l'algorithme demeure le même : pour chaque entier, on incrémente de un le compteur du nombre de diviseurs pour chacun de ses multiples.

Pour commencer, nous allons remplir le tableau des premiers nombres premiers en l'initialisant à l'infini pour que le résultat soit toujours compatible avec la procédure précédente :

for p from 1 to R do
    premiers[p] := infinity;
od:

Puis on déclare un tableau au nom plutôt explicite :

nbDiv := array(1..R):

Que l'on initialise avec la valeur false :

for p from 1 to R do
    nbDiv[p] := 1;
od:

Et que l'on renseigne comme suit : pour chaque entier dans l'ordre croissant qui n'a pas été marqué au préalable comme non-premier (et qui est donc premier, puisque ses éventuels diviseurs sont plus petits que lui), on « crible » tous ses multiples en les marquant comme "non-premiers". Seuls les survivants trouveront refuge dans le tableau résultat. Ainsi :

i := 1:
for p from 2 to R do
    nbDiv[p] := nbDiv[p] + 1;
    if (nbDiv[p] = 2) then
        premiers[i] := p;
    i := (i + 1);
    fi;
    j := p + p;
    while (j <= R) do
        nbDiv[j] := nbDiv[j] + 1;
        j := (j + p);
    od;
od:

Déterminer si un nombre est premier se fait alors en temps constant, de même que pour le nombre de diviseurs d'un entier :

nbrediv3 := proc(n)
    eval(nbDiv[n]);
end:

premier3 := proc(n)
    evalb(nbDiv[n] = 2)
end:

Complexité

On peut évaluer la complexité de ce crible on profitant de quelques probabilités : il existe environ $\small n \over p$ multiples de p inférieurs à n. Comme pour chaque mutliple on effectue de l'ordre d'une opération, on en déduit un ordre de grandeur $\small \eta$ du nombre d'opérations :

$\large \eta\,=\,\sum_{p=1}^{n} {\frac n p}\,=\,n \sum_{p=1}^n {\frac 1 p}\,\simeq\,n \ln(n)$

D'où une complexité finale en $\small O(n \ln(n))$ .

Somme des diviseurs

On s'intéresse désormais à la somme des diviseurs d'un entier donné. On va commencer par donner une variante directe du crible précédent qui permet une première résolution efficace avant de profiter d'un résultat établi pour proposer une solution un peu plus performante...

Variante du crible

Et si nous reprenions le principe du crible mais en ajoutant les diviseurs au lieu de les compter ? Cela pourrait donner une déclaration comme celle-ci :

sommeDiv := array(1..R):

Suivie d'une initialisation comme celle-là :

for p from 1 to R do
    sommeDiv[p] := 0;
od:

Et il ne reste plus qu'à "alléger" un peu le crible précédent tout en conservant l'essentiel :

i := 1:
for p from 1 to R do
    sommeDiv[p] := sommeDiv[p] + p;
    j := p + p;
    while (j <= R) do
        sommeDiv[j] := sommeDiv[j] + p;
        j := (j + p);
    od;
od:

L'écriture de la procédure associée est alors aisée :

robert := proc(n)
    sommeDiv[n];
end:

sommediviseurs := robert:

Et on vérifie bien par exemple :

sommediviseurs(5);

Voici donc une première solution en $\small O(n \ln(n))$ à ce problème.

Les mathématiques au secours

Réflexion

Cependant, quelques mathématiques permettent d'établir le résultat suivant (on ne le démontrera pas ici) : si on a

$\large n\,=\,\sum_{k=1}^{n} {{p_k}^{\alpha_k}}$

alors la somme des diviseurs est donnée par

$\large s(n)\,=\,\prod_{k=1}^{n} {{p_k^{\alpha_k+1}\,-\,1} \over {{p_k}\,-\,1}}$

où s(n) est la somme des diviseurs de n. Cette solution est moins simple mais plus efficace que les précédentes. L'idée est de retenir pour chaque n le plus grand diviseur premier de n ainsi que son exposant, c'est à dire un couple : $\small (p_n, \gamma_n)$ . L'intérêt de la chose est de permettre de reformuler la relation précédente sous forme d'une récurrence (car en informatique on aime bien les relations récursives) en posant si n > 1 :

$\large s(n)\,=\,{s({n \over {p_n^{\gamma_n}}}).{{{p_n}^{\gamma_n\,+\,1}\,-\,1} \over {p_n\,-\,1}}$

Par la suite, histoire de simplifier un peu ces notations que Maple a du mal à mettre en forme, on notera :

$\large \phi_n\,=\,{{{p_n}^{\gamma_n+1}\,-\,1} \over {p_n\,-\,1}}$

$\large \delta_n\,=\,{n \over {p_n^{\gamma_n}}}$

D'où :

$\large s(n)\,=\,\phi_n\,\cdot\,s(\delta_n)$

Une fois ces petites simplifications posées, on va tâcher de déterminer les différents $\small \phi_n$ et $\small \delta_n$ en criblant à nouveau, mais d'une façon un peu plus élaborée.

Implémentation

On commence par les classiques déclarations :

estPremier := array(1..R):
phi := array(1..R):
delta := array(1..R):

Puis on initialise les tableaux :

for p from 1 to R do
    estPremier[p] := true;
od:

Et on attaque le problème :

for p from 2 to R do
    if (estPremier[p] = false) then
        next;
    fi;
    phi[p] := p + 1;
    delta[p] := 1;
    pk := p;
    while (pk <= R) do
        phiActuel := ((p * pk - 1) / (p - 1));
        if (pk = p) then
            m := 2;
        else
            m := 1;
        fi;
        while (m * pk <= R) do
            i := m * pk;
            estPremier[i] := false;
            phi[i] := phiActuel;
            delta[i] := m;
            m := (m + 1);
        od;
        pk := pk * p;
    od;
end:

Une petite explication s'impose. Comme dans les cribles précédents, on a besoin de savoir si un entier est premier, et en réalité la boucle ne traite que les nombres premiers, d'où le test en tête de boucle. Une fois que l'on est sûr que l'entier p considéré est bien premier, on renseigne ses champs directement (on connait bien la valeur de s(n) pour un nombre premier) puis on traite chacune des puissances de p que l'on a noté pk.

Pour une puissance donnée, on calcule l'expression géométrique associée, puis on évalue tous les multiples d'un facteur m de cette puissance. Attention : à la première itération, on doit commencer par le double de p puisqu'on a déjà renseigné les champs de p. Ceci fait, pour chaque multiple, on renseigne les deux champs de sorte que le nombre premier actuel soit le prochain sur la liste de remontée. Comme on a parcouru les nombres premiers dans l'ordre croissant, on est sûr que cette liste chaînée sera bien ordonnée dans l'ordre inverse, et le calcul des s(n) peut alors se faire itérativement dans un tableau :

S := array(1..R):
S[1] := 1:
for n from 2 to R do
    S[n] := phi[n] * S[delta[n]];
od:

Et il ne reste plus qu'à définir jovialement la procédure de lecture :

sommediviseurs2 := proc(n)
    S[n];
end:

Complexité

Le calcul de la complexité de cette version n'est pas ce que l'on pourrait qualifier de trivial et je vais devoir faire appel à certains résultats que je ne sais pas établir, mais qui sont beaux en soit et me donnent une complexité convenable. Je regrette cette approche parfois un peu physicienne, mais les règles de manipulation de la notation O y invitent parfois un peu...

Pour un p fixé, observons l'évolution du nombre de calculs effectués pour chaque exposant. En notant $\small \eta$ le nombre d'opérations :

$\large \gamma\,=\,1 \qquad \Rightarrow \qquad \eta\,=\,{\frac R p}$

$\large \gamma\,=\,2 \qquad \Rightarrow \qquad \eta\,=\,{\frac R {p^2}}$

Et ainsi de suite. On voit que, plus l'exposant est grand, moins il y a d'opérations, sachant qu'une condition d'arrêt survient tôt ou tard. On peut donc approximer par excès ce nombre d'opérations par la somme d'une série géométrique de raison $\small \frac 1 p$ :

$\large \eta\,=\,{\frac R p} {\frac 1 {1\,-\,{\frac 1 p}}} \ \simeq \ {R \over p}$

Ce qui nous avance un peu. Reste à calculer la somme des inverses des nombres premiers, ce qui est moins évident. On introduit alors la probabilité qu'un entier p a d'être premier, et qui vaut :

$\large \rho\,=\,{\frac 1 {\ln p}}$

Et on peut alors approximer la somme $\small \sigma$ des inverses des nombres premiers entre 1 et R :

$\large \sigma\,=\,\sum_{p=2}^{R} {\frac 1 {p \ln p}}$

Que l'on peut brutalement approximer par l'intégrale associée en primitivant $\small {\frac 1 {p \ln p}}$ par $\small \ln (\ln (p))$ . Au final, on a une complexité de $\small {\cal O}(n \ln(\ln(n)))$ , ce qui achève cette démonstration.

Conclusion

En conclusion, on remarquera quand même qu'une approche mathématique intéressante du problème à permis d'aboutir à une solution plus performante, même si la différence ne sera sensible que pour de très grandes valeurs de n (de l'ordre du milliard par exemple).