Le crible d'Atkin

Le principe

Le crible d'Atkin est une méthode considérée comme une version optimisée du crible d'Eratosthène pour trouver tous les nombres premiers jusqu'à une borne fixée N. L'algorithme a été conçu récemment par A. O. L. Atkin et Daniel J. Bernstein dont vous trouverez la publication de leurs travaux "Prime sieves using binary quadratic forms" en suivant ce lien.

Afin de comprendre l'approche qu'on eut les auteurs sur la conception du crible, voici la traduction de la section "stratégie" de ce document :
"...L'idée de l'algorithme consiste à énumérer les valeurs de certaines formes quadratiques binaires irréductibles.
Par exemple, un entier positif non divisible par un carrée (square-free) p ∈ 1 + 4Z est premier si et seulement si l'équation 4x² + y² = p a un nombre impair de solutions positives (x, y). Nous couvrons tous les nombres premiers p > 3 comme suit. Pour p ∈ 1 + 4Z nous utilisons 4x² + y² avec x > 0 et y > 0, p ∈ 7 + 12Z nous utilisons 3x² + y² avec x > 0 et y > 0, pour p ∈ 11 + 12Z nous utilisons 3x² - y² avec x > y > 0...."

Nous voyons que les concepts mathématiques qui entrent en jeux sont beaucoup plus complexes qu'avec Eratosthène et il est malheureusement difficile de trouver de l'information pertinente sur le crible car très peu d'article de vulgarisation circulent sur cet algorithme. D'ailleurs ce manque d'information se retrouve également sur Wikipédia, il suffit d'aller sur la version anglaise (Sieve of Atkin) et française (Crible d'Atkin) pour constater que les contenus ne sont pas "synchronisés" et que, par exemple, la date de création du crible est différente entre les deux textes (2004 vs 1999). Peut-être dû a une réédition des travaux ?

Première implémentation

D'un point de vue algorithmique, commençons par implémenter la version trouvée en pseudo code dans le Wikipédia anglais histoire de comprendre comment ça fonctionne. Malgré le problème de contenu cité ci-dessus, le site semble être un bon point de départ pour appréhender Atkin.

private static final int MAX = 100000000;
private static final int SQRT_MAX = (int) Math.sqrt(MAX) + 1;
private static boolean[] array = new boolean[MAX];
//--//
for (int x = 1; x < SQRT_MAX; x++) {
  for (int y = 1; y < SQRT_MAX; y++) {
    int k = 4 * x * x + y * y;
    if ((k < MAX) && ((k % 12 == 1) || (k % 12 == 5))) {
      array[k] = !array[k];
    }
    k = 3 * x * x + y * y;
    if ((k < MAX) && (k % 12 == 7)) {
      array[k] = !array[k];
    }
    if (x > y) {
      k = 3 * x * x - y * y;
      if ((k < MAX) && (k % 12 == 11)) {
        array[k] = !array[k];
      }
    }
  }
}
array[2] = true;
array[3] = true;
for (int n = 5; n <= SQRT_MAX; n++) {
  if (array[n]) {
    int n2 = n * n;
    for (int k = n2; k < MAX; k += n2) {
      array[k] = false;
    }
  }
}

Au niveau des performances nous avons un petit espoir car les temps d'un Atkin non optimisé sont meilleurs que les temps d'un Eratosthène 1er du nom non optimisé puisque nous avons 6 secondes 9 pour Atkin contre 7 secondes 5 pour Eratosthène (dans le cas des premiers 100 millions).
Par contre la première optimisation d'Eratosthène "tape fort" car les temps de la version 2 sont divisés par 2, pourra-t-on faire mieux avec Atkin ? (angoisse)

MAX	Pi(MAX)	Temps (sec)
1 000 000	78 498	0.063
10 000 000	664 579	0.625
100 000 000	5 761 455	6.953

source

Implémentation, analyse des temps

Pour pouvoir optimiser Atkin, essayons d'analyser le programme et voir où il prend le plus de temps, pour cela nous allons d'abord couper la première boucle en 3 afin de pouvoir calculer les temps intermédiaires. Bien entendu l'éclatement de cette boucle n'influence pas le résultat, nous avons bien tous les premiers à la fin de l'exécution du programme.
Nous en profitons également pour redécouper les bornes supérieures des 3 boucles car la borne "racine carrée" de la boucle initiale était un poil trop large et obligeait l'ajout systématique du test "k < MAX". Cette première optimisation ne nous fait pas gagner de temps mais rattrape le temps perdu par le découpage en 3 boucles.

double debut = System.currentTimeMillis();
int xUpper = (int) Math.sqrt(MAX / 4) + 1;
for (int x = 1; x < xUpper; x++) {
  for (int y = 1; y < Math.sqrt(MAX - k1); y++) {
    Si (4x² + y²) % 12 == 1 ou (4x² + y²) % 12 == 5
      --> Alors array[k] = !array[k];
  }
}

double intermediaire1 = System.currentTimeMillis();
xUpper = (int) Math.sqrt(MAX / 3) + 1;
for (int x = 1; x < xUpper; x++) {
  for (int y = 2; y < Math.sqrt(MAX - k1); y++) {
    Si (3x² + y²) % 12 == 7 Alors array[k] = !array[k];
  }
}

double intermediaire2 = System.currentTimeMillis();
xUpper = (int) Math.sqrt(MAX);
for (int x = 1; x < xUpper; x++) {
  for (int y = 1; y < x; y++) {
    Si (3x² - y²) % 12 == 11 Alors array[k] = !array[k];
  }
}

double intermediaire3 = System.currentTimeMillis();
for (int n = 5; n <= SQRT_MAX; n++) {
  Elimination des multiples des carrés n², 2n², 3n², ..., MAX
}
double fin = System.currentTimeMillis();

MAX	Pi(MAX)	Temps (sec)	Boucle 1 (sec)	Boucle 2 (sec)	Boucle 3 (sec)	B. carrées(sec)
1 000 000	78 498	0.062	0.031	0.015	0.0	0.016
10 000 000	664 579	0.656	0.281	0.234	0.094	0.047
10 0000 000	5 761 455	6.984	3.047	2.328	1.062	0.547

Nous voyons que le temps de la boucle 1 est supérieur à la boucle 2, lui-même supérieur à la boucle 3 et lui-même supérieur à la boucle d'élimination des multiples des carrés. De plus la boucle 1 prend 43% du temps total d'exécution de l'algorithme, la boucle 2 33%, la boucle 3 15% et la boucle des carrées 7%.
Nous nous focaliserons donc sur une optimisation des 3 premières boucles car elles sont les plus consommatrices et aussi parce qu'elles se ressemblent (dans le sens une optimisation pour l'une sera sûrement applicable aux autres) et nous laisserons de côté celles des carrées. source

Deuxième implémentation

Pour optimiser l'algorithme, nous allons nous attaquer à l'opérateur modulo qui représente en informatique une opération couteuse en temps.
L'idée est de pouvoir se passer du test "modulo 12" présent dans les 3 boucles et pour cela nous allons essayer de résoudre les équations modulaires ou du moins trouver les solutions sous une forme plus simple "informatiquement parlant".

Que pouvons-nous dire des équations 4x² + y² = 1 / 5 [12], 3x² + y² = 7[12], 3x² - y² = 11[12] ?

Elles ont 2 inconnues x et y.
Elles sont modulaires.
Ce sont des équations diophantiennes, c'est a dire les coefficients sont des nombres entiers et les solutions recherchées sont également entières.
Elles sont non linéaires (à cause du degré 2).

Ne connaissant pas les outils mathématiques permettant de les résoudre, nous allons "tricher" en déroulant les équations à la main afin de voir sous quelle forme sont les solutions puis essayer d'en dégager une règle. C'est une sorte de reverse engineering des équations....

L'équation 3x² + y² = 7[12] de la 2eme boucle semble être un bon candidat pour commencer. source

    for (int x = 1; x < xUpper; x++) {
      int k1 = 3 * x * x;
      for (int y = 2; y < Math.sqrt(MAX - k1); y++) {
        int k = k1 + y * y;
        if (k % 12 == 7) {
          System.out.println("x = " + x + " y = " + y);
        }
      }
    }    x = 55 y = 9988
x = 55 y = 9992
x = 55 y = 9994
x = 55 y = 9998
x = 57 y = 2
x = 57 y = 4
x = 57 y = 8
x = 57 y = 10
x = 57 y = 14

Nous voyons que les solutions répondent aux règles suivantes :

x doit être impair.
y prend les valeurs suivante 2, 4, 8, 10, 14... c'est à dire que l'écart entre les solutions y semble suivre la séquence +2, +4, +2, +4....

Les solutions étant généralisées, nous pouvons donc optimiser la boucle 2 :

    int[] sequence2 = { 4, 2 };

    for (int x = 1; x < xUpper; x += 2) {
      int index = 0;
      int k1 = 3 * x * x;
      for (int y = 2; y < Math.sqrt(MAX - k1); y += sequence2[(++index & 1)]) {
        int k = k1 + y * y;
        array[k] = !array[k];
      }
    }

Et les améliorations sont au rendez-vous. source

MAX	Boucle 2 (sec)
1 000 000	0.0
10 000 000	0.047
100 000 000	0.5

L'optimisation donne de bon résultat puisque la boucle est 5 fois plus rapide que la précédente implémentation mais ne vous emballez pas car cette boucle se prêtait bien à cette modification...

Troisième implémentation

Maintenant que nous avons vu comment améliorer les temps de la boucle 2, nous pouvons faire de même sur les 2 autres boucles.
Nous obtenons les règles suivantes :

Pour la boucle 1

Il n'y a pas de règles spécifiques pour x.
Concernant y, l'équation admet "2 jeux de solution" selon que x est divisible par 3 ou non.
Lorsque x est divisible par 3, y prend les valeurs 1, 5, 7, 11, 13... c'est dire que l'écart semble suivre la séquence +4, +2, +4, +2....
Lorsque x est non divisible par 3, les solutions y sont impairs.

Pour la boucle 3

Il n'y a pas de règles spécifiques pour x.
Concernant y, l'équation admet "2 jeux de solution" selon que x est pair ou non.
Lorsque x est pair, l'écart des solutions y semble suivre la séquence +2, +4, +2, +4...
Lorsque x est impair, l'écart des solutions y semble suivre la séquence +4, +2, +4, +2....

On remarquera la séquence [+2,+4,...] (et son équivalent [+4,+2,...]) présent dans toutes les solutions y, nous en profitons pour mutualiser cet élément dans l'implémentation optimisée.
Et voici les résultats tant attendus...

MAX	Pi(MAX)	Temps (sec)	Boucle 1 (sec)	Boucle 2 (sec)	Boucle 3 (sec)	B. carrées(sec)
1 000 000	78 498	0.016	0.0	0.016	0.0	0.0
10 000 000	664 579	0.187	0.078	0.031	0.016	0.062
100 000 000	5 761 455	2.187	0.859	0.391	0.344	0.593

Enjoy ! Nous avons réussi à obtenir des meilleurs temps qu'un Eratosthène optimisé. A titre de comparaison, l'Eratosthène "deuxième implémentation" est à 3.063 secondes pour une borne à 100 millions alors qu'Atkin est à 2.187, nous avons donc un gain d'un peu moins de 50%. source

Quatrième implémentation

Continuons sur la lancé et essayons d'atteindre le milliard en implémentant Atkin avec un tableau de bit, c'est exactement la même approche que nous avions eu avec la 3ème implémentation d'Erastotene. Nous aurons besoin de 2 opérations de bit supplémentaires par rapport à Erastotene à savoir : unSetBit permettant de mettre à 0 un bit du tableau et toggleBit permettant de switcher un bit. source

MAX	Pi(MAX)	Temps (sec)	Boucle 1 (sec)	Boucle 2 (sec)	Boucle 3 (sec)	B. carrées(sec)
1 000 000	78 497	0.016	0.0	0.0	0.0	0.016
10 000 000	664 578	0.172	0.047	0.015	0.031	0.079
100 000 000	5 761 454	3.172	1.422	0.531	0.39	0.829
1 000 000 000	50 847 533	37.453	16.109	7.157	6.015	8.172

Et la c'est le drame, Atkin devient plus long ! Mais pourquoi ??
Pour l'instant c'est encore le flou mais la complexité de l'algorithme nous amène à 2 pistes :

Int vs long
Pour atteindre le milliard certains types java sont passés de int à long. Peut-être que ce changement est plus impactant que prévu et ressort dans les calculs plus complexes d'Atkin ? (je pense aux x² ou y²). C'est sans compter que les temps ont été mesurés sur une vielle machine 32 bits.
Trop d'appels aux opérations de manipulation des bits.
Les opérations setBit, unSetBit, toggleBit sont vraisemblablement couteuses en temps, peut-être qu'Atkin "abuse" et les appelle trop souvent ?

En conclusion

Nous voyons qu'Atkin peut aller plus vite qu'Eratosthène lorsqu'il est "savamment" optimisé mais que d'effort pour arriver à ce résultat.
Le code est beaucoup moins élégant qu'Eratosthène, tout ça pour gagner du temps.