Algorithmes arithmétiques en Maple

Dans cet article, nous allons implémenter quelques fonctions arithmétiques sous Maple : nombre de diviseurs, test de primalité, somme des diviseurs, etc. On s'intéressera dans un premier temps à des solutions naïves, puis on tentera d'aboutir à de meilleures complexités en temps, 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.

Algorithmes naïfs

Tout d'abord, initialisons notre feuille de calcul :

restart;
R := 42:

L'entier \(R\) que nous venons de définir servira de bornes aux tableaux statiques que nous allons allouer. Notez qu'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.

Nombre de diviseurs

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\). Déterminer le nombre de diviseurs de \(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é en temps est \(O(n)\), c'est à dire que le le nombre d'opérations élémentaires que nécessite la procédure est directement proportionnel à la valeur de l'entier \(n\) passé en paramètre (ici, nous avons \(n\) multiplications en tout). Cette solution n'est envisageable que pour \(n \leq 10^8\), valeur pour laquelle il faudrait déjà plus de 10 secondes à un ordinateur actuel 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 \(O(\sqrt{n})\) plus performante que la précédente, car le calcul peut désormais se faire pour \(n \leq 10^{16}\), toujours pour un processeur actuel aux alentours du Ghz.

Test de parfaite quadrature

Suivant la même approche exhaustive, on veut tester si un entier naturel \(n\) est de la forme \(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 \(O(\sqrt{n})\).

Algorithmes plus performants

Il est possible d'aboutir à une solution bien plus performante en faisant appel au crible d'Eratosthè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é. (Celle-ci n'est pas précisée dans la documentation au moment où nous rédigeons ces lignes en mars 2007. Edit: ni en octobre 2021.)

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 :

\begin{equation*} d(n) = \prod_{i=1}^{r}(\alpha_i + 1) \end{equation*}

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));
true

Mais la complexité de cette nouvelle mouture est en \(O(\sqrt{n})\) ce qui est bien plus performant que du \(O(n)\) comme l'illustrent les temps d'exécution suivants (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 \(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 du 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.

Première version

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 sous un nom de variable explicite :

nbDiv := array(1..R):

Que l'on initialise avec la valeur booléenne 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:

On peut évaluer la complexité de ce crible on profitant de quelques probabilités : il existe environ \(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 \(\eta\) du nombre d'opérations :

\begin{equation*} \eta = \sum_{p=1}^{n}{n \over p} = n \left(\sum_{p=1}^n{1 \over p}\right) \sim n \ln n \end{equation*}

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

Somme des diviseurs

On peut maintenant s'intéresser à 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.

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);
6

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

Crible encore plus performant

Quelques mathématiques permettent de dériver un crible encore plus performant. Celui-ci repose sur le théorème suivant (qu'on ne démontrera pas ici) :

Si \(n = \sum_{k=1}^{n}{{p_k}^{\alpha_k}}\), alors \(s(n) = \prod_{k=1}^{n} {{p_k^{\alpha_k+1} -1} \over {{p_k} - 1}}\).

\(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 \((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 \(1 < n\), alors \(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 :

\begin{align} \phi_n & = {{{p_n}^{\gamma_n + 1}-1} \over {p_n - 1}} & \delta_n & = {n \over {p_n^{\gamma_n}}} \end{align}

De sorte que \(s(n) = \phi_n s(\delta_n)\). Une fois ces simplifications posées, on va tâcher de déterminer les différents \(\phi_n\) et \(\delta_n\) en criblant à nouveau, mais d'une façon un peu plus élaborée. 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 &lt;= R) do
        phiActuel := ((p * pk - 1) / (p - 1));
        if (pk = p) then
            m := 2;
        else
            m := 1;
        fi;
        while (m * pk &lt;= R) do
            i := m * pk;
            estPremier[i] := false;
            phi[i] := phiActuel;
            delta[i] := m;
            m := (m + 1);
        od;
        pk := pk * p;
    od;
end:

Une 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 les puissances de \(p\) que l'on a nommées 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é du nouveau crible

Le calcul de la complexité de cette version n'est pas ce que l'on pourrait qualifier de trivial et nous allons faire appel à certains résultats que je ne sais pas établir, mais qui sont beaux en soi et donnent une complexité convenable. Je regrette cette approche un peu physicienne, mais les règles de manipulation de la notation \(O\) m'y invitent ;-)

Pour un \(p\) fixé, observons l'évolution du nombre de calculs effectués pour chaque exposant. En notant \(\eta\) le nombre d'opérations : pour \(\gamma = 1\), nous avons \(R/p\) opérations ; pour \(\gamma = 2\), nous avons \(R/{p^2}\) opérations, etc. 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 approximer par excès ce nombre d'opérations par la somme d'une série géométrique :

\begin{equation*} \eta = \frac{R}{p} \sum_{k=0}^{+\infty} \frac{1}{p^k} = {R \over p} {1 \over {1 - {\frac 1 p}}} \sim {R \over p} \end{equation*}

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 :

\begin{equation*} \mathbb{P}[p \ \mathrm{premier}] = {1 \over {\ln p}} \end{equation*}

On peut ainsi approximer la somme des inverses des nombres premiers entre \(1\) et \(R\) :

\begin{align} \sum_{p \ \mathrm{premier} \leq R} \frac{1}{p} & = \sum_{p = 2}^R \frac{1}{p} \mathbb{1}[p \ \mathrm{premier}] \\ & \approx \mathbb{E} \left( \sum_{p = 2}^R \frac{1}{p} \mathbb{1}[p \ \mathrm{premier}] \right) \\ & = \sum_{p = 2}^R \frac{1}{p} \mathbb{E}\left(\mathbb{1}[p \ \mathrm{premier}]\right) \\ & = \sum_{p=2}^R \frac{1}{p} \mathbb{P}[p \ \mathrm{premier}] \\ & \approx \sum_{p=2}^{R}{\frac 1 {p \ln p}} \end{align}

Que l'on peut brutalement approximer par l'intégrale associée en primitivant \(\ln \ln p\). Au final, on a une complexité de \(O(n \ln \ln n)\), ce qui achève une démonstration plus fatiguante que la rédaction de l'algorithme lui-même ! In fine, on retiendra qu'une approche mathématique du problème à permis d'aboutir à une solution plus performante, même si la différence ne sera sensible que pour des valeurs de \(n\) de l'ordre du milliard.

© Stéphane Caron — All content on this website is licensed under the CC BY 4.0 license.
π