Nombres premiers: crible des multiplications

Description


Nombres premiers: crible des multiplications


Le principe de l'algorithme consiste à parcourir la table de multiplication pour neutraliser dans le crible les éléments qui correspondent à un produit.

Avantage de l'algorithme: temps d'exécution rapide et programmation simple (à part les dernières optimisations).
Désavantage: utilise un tableau (bool) de tous les nombres candidats. (voir: Optimiasation C)

Cet article correspond à un essai personnel (et peut-être original ?).
Connaitriez-vous d'autres articles similaires ? Ai-je réinventé la roue ?
Proposez-vous un autre nom que "crible des multiplications" ou "crible des produits" ?

C'est par hasard que j'ai remarqué dernièrement que l'étape Mult4.js correspond à peu près à l'algorithme connu sous le nom de Crible de Sundaram (fichier Summary.rtf).
Cet algorithme sera spécifiquement traité sous Mult5.js (classique) et Mult6.js (optimisé).

Une version actualisée de l'introduction aux nombres premiers, complétée de mesures des temps d'exécution, peut être obtenue ici Nombres premiers.
Autres liens sur les nombres premiers:
de division (Contient la définition)
Essais de division avec liste
Crible d'Eratosthène
Crible d'Eratosthène (optimisation C)
Crible d'Atkin
 
 

Mult0.js


Voici pour commencer une interprétation très naïve de l'algorithme basé sur la table de multiplication:
//     i 2  3  4  5  6  7  8  9 ...
//   j
//   2   4  6  8 10 12 14 16 18   
//   3   6  9 12 15 18 21 24 27
//   4   8 12 16 20 24 28 32 36
//   5  10 15 20 25 30 35 40 45
//   6  12 18 24 30 36 42 48 54
//   7  14 21 28 35 42 49 56 63
//   8  16 24 32 40 48 56 64 72
//   9  18 27 36 45 54 63 72 81
//   ...

function Mult(mx) { // Mult0.js: naïf
  var i,j,n,sv=[false,false]; // sv=new Array(mx+1)
  for (n=2; n<=mx; n++) sv[n]=true; // 0,1 non premier
  for (j=2; j<=mx; j++)
    for (i=2; i<=mx; i++) if (i*j<=mx) sv[i*j]=false;
  return sv;
}

Test direct
Ne mettez pas max > 10000, car ce code est très très lent.
Ne vaudrait-il pas mieux d'abandonner ici le développement de cet algorithme ?
 
 

Mult1.js


Malgré tout, nous allons continuer et y apporter plusieurs améliorations "évidentes":

a) La boucle extérieure peut finir avec j*j<=mx.
b) Nous pouvons enlever tous les "doublons" en commençant la boucle interne depuis la diagonale (i=j).
c) La boucle interne peut se terminer lorsque i*j>mx.
//     i 2  3  4  5  6  7  8  9 ...
//   j
//   2   4  6  8 10 12 14 16 18   
//   3      9 12 15 18 21 24 27
//   4        16 20 24 28 32 36
//   5           25 30 35 40 45
//   6              36 42 48 54
//   7                 49 56 63
//   8                    64 72
//   9                       81
//   ...

function Mult(mx) { // Mult1.js:
  var i,j,n,sv=[false,false]; // sv=new Array(mx+1)
  for (n=2; n<=mx; n++) sv[n]=true;
  for (j=2; j*j<=mx; j++)
    for (i=j; i*j<=mx; i++) sv[i*j]=false;
  return sv;
}

Test direct
Comme l'amélioration est fulgurante, nous continuons !
 
 

Mult2.js


La table de multiplication ne doit pas forcément être calculée avec des multiplications; on peut le faire avec des additions:
function Mult(mx) { // Mult2.js:
  var j,n,sv=[false,false]; // sv=new Array(mx+1)
  for (n=2; n<=mx; n++) sv[n]=true;
  for (j=2; j*j<=mx; j++)
    for (n=j*j; n<=mx; n+=j) sv[n]=false;
  return sv;
}

Test direct
 
 

Mult3.js


Ne traitons la ligne j que lorsque j est premier et remplaçons le test j*j<=mx par j<=sqrt(mx):
function Mult(mx) { // Mult3.js:
  var j,n,r=Math.floor(Math.sqrt(mx)),sv=[false,false]; // sv=new Array(mx+1)
  for (n=2; n<=mx; n++) sv[n]=true;
  for (j=2; j<=r; j++)
    if (sv[j]) for (n=j*j; n<=mx; n+=j) sv[n]=false;
  return sv;
}

Test direct
 
 

Mult4.js


Le prochain but est de ne consirérer que les nombres impairs.

Comme pour "Sieve2.js", faisons correspondre à un élément du tableau sv non pas
l'index p lui-même, mais, pour p>0, la valeur impaire pp=1+2*p (donc p=(pp-1)/2):
//  p:  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 ...
// pp:  2  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 ...
//
// Adaptons la table de multiplication pp,     et les indexes p correspondants:
//   ii  3  5  7  9  11  13  15  17  19 ...     i 1  2  3  4  5  6   7   8   9 ...
//  jj                                         j
//   3   9 15 21 27  33  39  45  51  57 ...    1  4  7 10 13 16 19  22  25  28 ...
//   5     25 35 45  55  65  75  85  95        2    12 17 22 27 32  37  42  47
//   7        49 63  77  91 105 119 133        3       24 31 38 45  52  59  66
//   9           81  99 117 135 153 171        4          40 49 58  67  76  85
//  11              121 143 165 187 209        5             60 71  82  93 104
//  13                  169 195 221 247        6                84  97 110 123
//  15                      225 255 285        7                   112 127 142
//  17                          289 223        8                       144 161
//  19                              361        9                           180
//  ...                                        ...

On constate qu'on peut directement "travailler" sur les indexes, ils augmentent de 1+2*j dans les lignes, et de 4*(1+j)=n dans la diagonale d:
function Mult(mx) { // Mult4.js:
  var h=1+Math.floor((mx-1)/2),d,j,jj,p,sv=[]; // sv=new Array(h)
  for (k=0; k<h; k++) sv[k]=true;
  for (j=1,d=4; d<h; ++j,d+=4*j) // traitement par lignes
    if (sv[j]) {
      jj=1+2*j;
      for (p=d; p<h; p+=jj) sv[p]=false;
    }
  return sv;
}

Test direct
Cette fonction semble être plus rapide que celles traitées sous Nombres premiers: Crible d'Eratosthène.
 
 

Mult5.js


Comme annoncé au début, étudions l'algorithme du Crible de Sundaram.

Soit In = 2*n + 1 un entier impair quelconque >= 3, avec n >= 1.
Soit Ip = In + 2*k un entier impair quelconque >= In, avec k >= 0.

On peut obtenir tous les produits de nombres impairs
In*Ip = In*(In+2*k) = (2*n+1)*(2*n+1+2*k)
en incrémentant n et k.

Utilisons comme précédemment:
//  p: 0 1 2 3 4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 ...
// pp: 2 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 ...
//
// Voici le début du tableau obtenu pour pp=In*Ip, et les indexes p correspondants:
//   nn  3  5   7   9  11  13  15  17  19        n   1  2  3  4   5   6   7   8   9
// kk                                           k 
//  3    9                                      1    4
//  5   15 25                                   2    7 12
//  7   21 35  49                               3   10 17 24
//  9   27 45  63  81                           4   13 22 31 40
// 11   33 55  77  99 121                       5   16 27 38 49  60
// 13   39 65  91 117 143 169                   6   19 32 45 28  71  84
// 15   45 75 105 135 165 195 225               7   22 37 52 67  82  97 112
// 17   51 85 119 153 187 221 255 289           8   25 42 59 76  93 110 127 144
// 19   57 95 133 171 209 247 285 323 361       9   28 47 66 85 104 123 142 161 180

On obtient les tableaux transposés décrits sous Mult4.

Sur le web, les "meilleurs" codes que j'ai trouvés correspondent à peu près à:
function Mult(mx) { // Mult5.js: Sundaram du web
  var h=1+Math.floor((mx-1)/2),k,kk,n,nn,p,pp,sv=[]; // sv=new Array(h)
  for (p=0; p<h; p++) sv[p]=true;
  for (n=1; n<h; n++) { // traitement par colonnes
    nn=2*n+1;
    for (k=n; k<h; k++) {
      kk=2*k+1;
      pp=nn*kk;
      p=(pp-1)/2;
      if (p>=h) break;
      sv[p]=false;
    }
  }
  return sv;
}

Test direct Sundaram
 
 

Mult6.js


Limitons la boucle extérieure à nn*nn<=mx et traitons la ligne n que lorsque sv[n]==true (nn est premier).
function Mult(mx) {
  var h=1+Math.floor((mx-1)/2),k,kk,n,nn,p,pp,sv=[]; // sv=new Array(h)
  for (p=0; p<h; p++) sv[p]=true;
  for (n=1,nn=3; nn*nn<=mx; n++,nn+=2) { // nn=2*n+1
    if (sv[n])
      for (k=n,kk=nn; k<h; k++,kk+=2) { // kk=2*k+1
        p=(nn*kk-1)/2;
        if (p>=h) break;
        sv[p]=false;
      }
  }
  return sv;
}

Puis simplifions les calculs de la boucle intérieure, ce qui revient à faire les calculs directement sur les indexes.
function Mult(mx) { // Mult6: Sundaram optimisé
  var h=1+Math.floor((mx-1)/2),n,nn,p,sv=[]; // sv=new Array(h)
  for (p=0; p<h; p++) sv[p]=true;
  for (n=1,nn=3; nn*nn<=mx; n++,nn+=2) // nn=2*n+1
    if (sv[n]) for (p=(nn*nn-1)/2; p<h; p+=nn) sv[p]=false;
  return sv;
}

Test direct Sundaram optimisé
 
 

Mult7.js


Revenons à Mult4 et éliminons maintenant tous les nombres divisibles par 2 et 3, (voir aussi Cycle0.js sous Nombres premiers: Crible avec cycle additif)
A partir de p=5, on a p+=c[i%2], avec le "cycle" c=[2,4].
Adaptons également les tables de multiplication:
// n:  0 1 2 3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 ...
// p:  2 3 5 7 11 13 17 19 23 25 29 31 35 37 41 43 47 49 53 55 59 61 65 67 71 ...
//
//  ii 5  7  11  13  17  19  23  25  29  31   35   37   41   43   47  ...
// jj
//  5 25 35  55  65  85  95 115 125 145 155  175  185  205  215  235
//  7    49  77  91 119 133 161 175 203 217  245  259  287  301  329
// 11       121 143 187 209 253 275 319 341  385  407  451  473  517
// 13           169 221 247 299 325 377 403  455  481  533  559  611
// 17               289 323 391 425 493 527  595  629  697  731  799
// 19                   361 437 475 551 589  665  703  779  817  893
// 23                       529 575 667 713  805  851  943  989 1081
// 25                           625 725 775  875  925 1025 1075 1175
// 29                               841 899 1015 1073 1189 1247 1363
// 31                                   961 1085 1147 1271 1333 1457
// 35                                       1225 1295 1435 1505 1645
// 37                                            1369 1517 1591 1739
// 41                                                 1681 1763 1927
// 43                                                      1849 2021
// 47                                                           2209
//     ...
//
//   i 2  3   4   5   6   7   8   9  10  11   12   13   14   15   16    u
//  j                                                                   v   [x,y]
//  2  9 12  19  22  29  32  39  42  49  52   59   62   69   72   79    8  [ 3, 7]
//  3    17  26  31  40  45  54  59  68  73   82   87   96  101  110   24  [ 5, 9]
//  4        41  48  63  70  85  92 107 114  129  136  151  158  173   16  [ 7,15]
//  5            57  74  83 100 109 126 135  152  161  178  187  204   40  [ 9,17]
//  6                97 108 131 142 165 176  199  210  233  244  267   24  [11,23]
//  7                   121 146 159 184 197  222  235  260  273  298   56  [13,25]
//  8                       177 192 223 238  269  284  315  330  361   32  [15,31]
//  9                           209 242 259  292  309  342  359  392   72  [17,33]
// 10                               281 300  339  358  397  416  455   40  [19,39]
// 11                                   321  362  383  424  445  486   88  [21,41]
// 12                                        409  432  479  502  549   48  [23,47]
// 13                                             457  506  531  580  104  [25,49]
// 14                                                  561  588  643   56  [27,55]
// 15                                                       617  674  120
// 16                                                            737
//     ...

Remarque: j'ai bien sûr fait un petit programme pour afficher ces deux tableaux.

Le nombre de candidats<=mx qui restent est est un peu difficile à calculer:
h = ((mx-5-a)/3+((a>1)? 4: 3)) , où a = (mx-5)%6.

Dans le tableau des indexes, après avoir noté sous u,v l'augmentation des éléments de la diagonale, observons que u,v augmentent alternativement de 8,16.
Puis, notons l'augmentation cyclique [x,y] des éléments dans les lignes: nous constatons que x=2*j-1 et que y vaut alternativement 4*j-1, 4*j-3.

function Mult(mx) { // Mult7.js: mx>9
  var a=(mx-5)%6,h=((mx-5-a)/3+((a>1)?4:3)),d,j,k,sv=[],u,v,xy=[]; // sv=new Array(h)
  for (k=0; k<h; k++) sv[k]=true;
  for (j=2,d=9,u=0,v=8; d<h; d+=(++j&1)? u+=8: v+=16) {
    if (sv[j]) {
      if (j&1) {xy[0]=2*j-1; xy[1]=4*j-3;} else {xy[0]=4*j-1; xy[1]=2*j-1;}
      for (k=d; k<h; k+=xy[k&1]) sv[k]=false; // k&1 équivalent à k%2
    }
  }
  return sv;
}

Test direct

Une fois le crible calculé, on peut en extraire la liste des nombres premiers p (à partir de 5) avec la boucle
for (i=2,p=5; i<h; p+=cyc[i&1],i++) if (sv[i]) { /* p is prime */ }
 
 

Conclusions


Malgré un début très décevant, nous arrivons finalement à un code tout à fait "honnête", qui supporte la comparaison avec les "meilleurs".
 
 

Mesures


Le tableau ci-dessous n'a qu'une valeur comparative.
Les mesures dépendent de l'ordinateur, du navigateur et de leur charge.
Sur le même ordinateur, j'ai chaque fois effectué une série de tests et indiqué le meilleur résultat.
Voir aussi Nombres premiers (fichier Summary.rtf)

Max: 1000 10000 100000 1000000 10000000 100000000
n:    168  1229   9592   78498   664579   5761455    memory:

Div0    2    63   3365       -        -         -    0
Div1    1    42   1770  145128        -         -    0
Div2    0     2     35     309     7579         -    0
Div3    0     2     31     305     6799         -    0
Div4    0     1     30     296     6564         -    0

List0   0     1     26     323     7301         -    int[n]
List1   0     1     19     172     2557         -    int[n]

Sieve0  0     1     11      79      659      7104    bool[Max]
Sieve1  0     1      9      63      466      5003    bool[Max]
Sieve2  0     1      2      20      142      1495    bool[Max/2]
Sieve3  0     0      1      23      138      1455    bool[Max/2]

Atkin0  0     0      5      49      407      4847    bool[Max]
Atkin1  0     0      5      47      384      4777    bool[Max]
Atkin2  0     0      5      37      360      4154    bool[Max]

Mult0  12   456 140532       -        -         -    bool[Max]
Mult1   0     0      4      48      787     10139    bool[Max]
Mult2   0     0      5      36      776     10091    bool[Max]
Mult3   0     0      3      25      226      2616    bool[Max]
Mult4   0     0      2      17      111      1307    bool[Max/2]
Mult5 Su      0      3      51      491      5356    bool[Max/2]
Mult6 Su-Opt  0      1      19      119      1363    bool[Max/2]
Mult7   0     0      1       8       77       850    bool[Max/3]

 
 
Testé avec:
Firefox            23.0     OK 
Google Chrome      28.0     OK 
Internet Explorer  10.0.7   OK 
Opera              12.16    OK 
Safari              5.1.7   OK 

Codes Sources

A voir également

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.