Voici une variante du crible d'Eratosthène qui a le mérite d'être itérative :
n := 1
S := dictionnaire de (clé : nombre, valeur : ensemble de nombres) vide
tant que n < n_max faire
n := n + 1
si n est une clé dans S
alors
// on en déduit que n est un nombre composé c
soit c := n
pour chaque p dans S(c), faire
soit cc := c + p
si la clé cc existe dans S
alors
S(cc) := S(cc) union { p }
sinon
S(cc) := { p }
supprimer l'entrée de clé c dans S
sinon
// on en déduit que n est un nombre premier p
soit p := n
S(2p) := { p }
La variable S représente l'état des calculs ; on peut la voir comme le stockage simultané
Digression : cf. A350809 sur OEIS.
Illustration : \( n \) est le temps ; \( S \) stocke chaque cercle (= chaque nombre premier) ainsi que la position courante sur chaque cercle (= le prochain multiple du nombre premier que représente le cercle),
En quelque sorte, on jongle avec les nombres premiers à partir de l'instant \( t = 2 \) :
Je me dis qu'il y a peut-être moyen de mathématiser l'algorithme.
Influencés par les séries génératrices, imaginons que le dictionnaire \( S \) soit un polynôme à deux variables (\( x \) et \( y \)) où les puissances de \( x \) marquent les \( c - n \) tandis que les puissances de \( y \) marquent les \( p \). L'algorithme peut être réalisé comme suit :
n := 1 S := 0 tant que n < n_max faire n := n + 1 S := S / x Soit Q := S(0, y) un polynôme en seulement y si Q != 0 alors // on en déduit que n est un nombre composé S := S - Q + Q(x * y) sinon // on en déduit que n est un nombre premier p S := S + x^n * y^n
Implémentation en PARI :
S(n)=if(n==1,0,my(T=S(n-1)/'x,Q=subst(T,'x,0));T+if(Q==0,'x^n*'y^n,subst(Q,'y,'x*'y)-Q)) for(n=1,15,print1(n," : ",S(n),"\n")) 1 : 0 2 : y^2*x^2 3 : y^3*x^3 + y^2*x 4 : (y^3 + y^2)*x^2 5 : y^5*x^5 + (y^3 + y^2)*x 6 : y^5*x^4 + y^3*x^3 + y^2*x^2 7 : y^7*x^7 + y^5*x^3 + y^3*x^2 + y^2*x 8 : y^7*x^6 + (y^5 + y^2)*x^2 + y^3*x 9 : y^7*x^5 + y^3*x^3 + (y^5 + y^2)*x 10 : y^5*x^5 + y^7*x^4 + (y^3 + y^2)*x^2 11 : y^11*x^11 + y^5*x^4 + y^7*x^3 + (y^3 + y^2)*x 12 : y^11*x^10 + (y^5 + y^3)*x^3 + (y^7 + y^2)*x^2 13 : y^13*x^13 + y^11*x^9 + (y^5 + y^3)*x^2 + (y^7 + y^2)*x 14 : y^13*x^12 + y^11*x^8 + y^7*x^7 + y^2*x^2 + (y^5 + y^3)*x 15 : y^13*x^11 + y^11*x^7 + y^7*x^6 + y^5*x^5 + y^3*x^3 + y^2*x
En particulier, quand ils sont évalués en \( x = 1 \), ces polynômes donnent les nombres premiers :
for(n=1,15,print1(n," : ",subst(S(n),'x,1),"\n")) 1 : 0 2 : y^2 3 : y^3 + y^2 4 : y^3 + y^2 5 : y^5 + y^3 + y^2 6 : y^5 + y^3 + y^2 7 : y^7 + y^5 + y^3 + y^2 8 : y^7 + y^5 + y^3 + y^2 9 : y^7 + y^5 + y^3 + y^2 10 : y^7 + y^5 + y^3 + y^2 11 : y^11 + y^7 + y^5 + y^3 + y^2 12 : y^11 + y^7 + y^5 + y^3 + y^2 13 : y^13 + y^11 + y^7 + y^5 + y^3 + y^2 14 : y^13 + y^11 + y^7 + y^5 + y^3 + y^2 15 : y^13 + y^11 + y^7 + y^5 + y^3 + y^2
Utilité ?
Peut-on mathématiser davantage (notamment le "if") ?
A POURSUIVRE...
LR, 27/02/2023.
Choisir une base \( x \) et essayer de réaliser l'algorithme à l'aide d'écritures en base \( x \) derrière la virgule,
chaque nombre premier \( p \) inférieur ou égal à \( n \) contribuant à fournir des chiffres 1 espacés de \( p \) et déphasés de \( p - (n \mbox{ mod } p) \).
Par exemple, écrire :
Si n = 5
Contribution de 2 : 0.1010101010101010101010101010101010... Contribution de 3 : 0.1001001001001001001001001001001001... Contribution de 5 : 0.0000100001000010000100001000010000... |
Si n = 6
Contribution de 2 : 0.0101010101010101010101010101010101... Contribution de 3 : 0.0010010010010010010010010010010010... Contribution de 5 : 0.0001000010000100001000010000100001... |
et additionner verticalement ? (doute... le fait d'additionner détruit-il de l'information ?). Tentative de formalisation :
\[ S(n) = \sum_{p \leq n} \frac{1}{x^{p - (n \mbox{ mod } p)}} \sum_{k \geq 0} \frac{1}{x^{kp}} \] \[ S(n) = \sum_{p \leq n} \frac{1}{x^{p - (n \mbox{ mod } p)}} \sum_{k \geq 0} \left ( \frac{1}{x^{p}} \right ) ^k \] \[ S(n) = \sum_{p \leq n} \frac{1}{x^{p - (n \mbox{ mod } p)}} \left ( \frac{1}{1 - \left ( \frac{1}{x} \right )^p } \right ) \] \[ S(n) = \sum_{p \leq n} \frac{1}{x^{p - (n \mbox{ mod } p)}} \left ( \frac{x^p}{x^p - 1 } \right ) \] \[ S(n) = \sum_{p \leq n} \frac{x^{(n \mbox{ mod } p)}}{x^p - 1 } \]Intéressant : apparition de polynômes cyclotomiques...
A quoi ressemblent les 10 premières fractions rationnelles S(n), développées ?
S(n)=my(s);forprime(p=2,n,s+=x^(n%p)/(x^p-1));s for(n=1,10,print1(n," : ",S(n),"\n")) 1 : 0 2 : 1/(x^2 - 1) 3 : (x^3 + x^2 + 2*x + 1)/(x^4 + x^3 - x - 1) 4 : (2*x^2 + 2*x + 1)/(x^4 + x^3 - x - 1) 5 : (2*x^7 + 4*x^6 + 5*x^5 + 5*x^4 + 6*x^3 + 5*x^2 + 3*x + 1)/(x^8 + 2*x^7 + 2*x^6 + x^5 - x^3 - 2*x^2 - 2*x - 1) 6 : (x^6 + 3*x^5 + 6*x^4 + 7*x^3 + 7*x^2 + 5*x + 2)/(x^8 + 2*x^7 + 2*x^6 + x^5 - x^3 - 2*x^2 - 2*x - 1) 7 : (x^13 + 4*x^12 + 10*x^11 + 17*x^10 + 24*x^9 + 29*x^8 + 32*x^7 + 33*x^6 + 32*x^5 + 27*x^4 + 20*x^3 + 12*x^2 + 5*x + 1)/(x^14 + 3*x^13 + 5*x^12 + 6*x^11 + 6*x^10 + 5*x^9 + 3*x^8 - 3*x^6 - 5*x^5 - 6*x^4 - 6*x^3 - 5*x^2 - 3*x - 1) 8 : (x^13 + 5*x^12 + 11*x^11 + 18*x^10 + 24*x^9 + 29*x^8 + 33*x^7 + 35*x^6 + 32*x^5 + 26*x^4 + 18*x^3 + 10*x^2 + 4*x + 1)/(x^14 + 3*x^13 + 5*x^12 + 6*x^11 + 6*x^10 + 5*x^9 + 3*x^8 - 3*x^6 - 5*x^5 - 6*x^4 - 6*x^3 - 5*x^2 - 3*x - 1) 9 : (2*x^13 + 6*x^12 + 12*x^11 + 18*x^10 + 24*x^9 + 30*x^8 + 35*x^7 + 35*x^6 + 31*x^5 + 24*x^4 + 16*x^3 + 9*x^2 + 4*x + 1)/(x^14 + 3*x^13 + 5*x^12 + 6*x^11 + 6*x^10 + 5*x^9 + 3*x^8 - 3*x^6 - 5*x^5 - 6*x^4 - 6*x^3 - 5*x^2 - 3*x - 1) 10 : (2*x^12 + 6*x^11 + 12*x^10 + 20*x^9 + 29*x^8 + 35*x^7 + 37*x^6 + 34*x^5 + 28*x^4 + 21*x^3 + 14*x^2 + 7*x + 2)/(x^14 + 3*x^13 + 5*x^12 + 6*x^11 + 6*x^10 + 5*x^9 + 3*x^8 - 3*x^6 - 5*x^5 - 6*x^4 - 6*x^3 - 5*x^2 - 3*x - 1)
Valeurs concrètes des 10 premières \( S(n) \) pour quelques bases \( x \) ?
En base 2 :
for(n=1,10,print1(n, " : ",subst(S(n),x,2), " = ", subst(S(n),x,2)+0.0,"\n")) 1 : 0 = 0.E-38 2 : 1/3 = 0.33333333333333333333333333333333333333 3 : 17/21 = 0.80952380952380952380952380952380952381 4 : 13/21 = 0.61904761904761904761904761904761904762 5 : 827/651 = 1.2703533026113671274961597542242703533 6 : 352/651 = 0.54070660522273425499231950844854070660 7 : 90059/82677 = 1.0892872261935000060476311428813333817 8 : 97441/82677 = 1.1785744523870000120952622857626667634 9 : 112205/82677 = 1.3571489047740000241905245715253335269 10 : 59056/82677 = 0.71429780954800004838104914305066705372
En base 3 :
for(n=1,10,print1(n, " : ",subst(S(n),x,3), " = ", subst(S(n),x,3)+0.0,"\n")) 1 : 0 = 0.E-38 2 : 1/8 = 0.12500000000000000000000000000000000000 3 : 43/104 = 0.41346153846153846153846153846153846154 4 : 25/104 = 0.24038461538461538461538461538461538462 5 : 9127/12584 = 0.72528607755880483153210425937698664971 6 : 2213/12584 = 0.17585823267641449459631277813095994914 7 : 7262719/13754312 = 0.52803215457087202907713595561886337899 8 : 8033845/13754312 = 0.58409646371261608723140786685659013697 9 : 10347223/13754312 = 0.75228939113784826169422360056977041091 10 : 3533045/13754312 = 0.25686817341354478508267080170931123273
Noter la valeur de \( S(n) \) en \( x = 0 \) : \( - \omega(n) \), où \( \omega(n) \) désigne le nombre de nombre premiers distincts dans la décomposition en facteurs premiers de \( n \).
Algorithme pas clair
n := 1 S := 0 tant que n < n_max faire n := n + 1 S := S * x si quoi ? alors n est un nombre composé Comment retrouver dans S la liste des p dont n est composé ? Comment mettre à jour S ? sinon n est un nombre premier S := S + 1 / (x^n - 1)
A POURSUIVRE...
LR, 19/03/2023.