CRIBLE D'ERATOSTHÈNE OPTIMISÉ

/*  Crible d'Eratosthène par pgl10  */

#include <stdio.h>   // pour printf
#include <stdlib.h>  // pour system
#include <string.h>  // pour memset
#include <time.h>    // pour clock
#include <math.h>    // pour sqrt

typedef unsigned long ulong;

#define ulmax 4294967295  // (2^32)-1 : valeur maximum des unsigned long 

ulong nmax;    // borne supérieure de l'intervalle [1,nmax] de calcul utilisé
char* crible;  // crible[1+(nmax-1)/16] : une suite de (1+nmax)/2 bits utiles
//   quand le calcul sera fini les nombres pairs ne seront pas marqués
//   et le (1+i)/2-ième bit de crible mis à 1 quand i est premier et 0 si non.
//   crible[0] pour les impairs de 1 à 15, crible[1] pour les impairs de 17 à 31, etc.

void pause() {
    printf("\n");
    system("pause");    
}
void comp(ulong i) { // mettre à 0 le (1+i)/2-ième bit du crible ( i impair et composé )
    ulong q, r;
    char c[9]={'\xFF', '\xFE', '\xFD', '\xFB', '\xF7', '\xEF', '\xDF', '\xBF', '\x7F'};
//  if(i%2==0) return;         // inutile : tous les nombres pairs ne sont pas marqués
    q=i>>4;                    // q=i/16=(x-1)/8 avec x=(1+i)/2 ou x=1+i/2 pour i impair
    r=((1+i)>>1)-(q<<3);       // r=(1+i)/2-8*q=x-8*q avec x=(1+i)/2  ( r de 1 à 8 )
    crible[q]=crible[q]&c[r];  // on met à 0 le bit (1+i)/2 du crible quand i est impair
}
bool prem(ulong i) { // retourner true si le nombre i est premier et false si non
    ulong q, r;
    if(i<2 || i>nmax) return false;
    if(i==2) return true;
    if(i%2==0) return false;
    q=i>>4;                    // q=0 si i de 1 à 15, q=1 si i de 17 à 31, etc.
    r=((1+i)>>1)-(q<<3);
    if(crible[q]>>(r-1)&1) return true;  // si le bit (1+i)/2 de crible est à 1
    return false;
}
void Eratosthene(ulong n) {       // crible d'Eratosthène
    ulong i, j, m, q, r, s, t;
    char c[9]={'\xFF', '\xFE', '\xFD', '\xFB', '\xF7', '\xEF', '\xDF', '\xBF', '\x7F'};
    memset(crible, '\xFF', 1+(n-1)/16);      // initialisation de crible[]
    comp(1);                                 // 1 n'est pas premier
    for(j=9; j<=n; j=j+6) {                  // les multiples de 2 ne sont pas marqués
        q=j>>4;                              // pour marquer les multiples impairs de 3
        r=((1+j)>>1)-(q<<3);                 // on évite d'appeler comp(j) pour gagner
        crible[q]=crible[q]&c[r];            // un peu de temps             
    }
    m=1+(ulong)sqrt((double)n);
    s=4;
    for(i=5; i<m; i=i+s) {                   // on marque les multiples des autres 
        if(prem(i)) {                        // nombres premiers connus en évitant
            if(((i+2)*i)%3) t=4; else t=2;   // les multiples de 2 et 3
            for(j=i*i; j<=n ; j=j+t*i ) { 
                q=j>>4;                      // ici j est toujours impair
                r=((1+j)>>1)-(q<<3);
                if(crible[q]>>(r-1)&1) crible[q]=crible[q]&c[r];
                t=6-t;
//              if(j>ulmax-t*i) {  // pour éviter un dépassement de capacité des ulong
//                  printf("\nLa limite de l'intervalle de calcul est trop grande\n"); 
//                  pause();       // pour gagner du temps on peut éviter ce contrôle, 
//                  exit(1);       // mais il faut veiller à utiliser nmax < 4294709603
//              }
            }
        }
        s=6-s;
    }
}
int main (int argc, char *argv[]) { 
    ulong i, np=0, p1=0, pn;
    char limite[21];
    if(argc!=2) {
        printf("\nPour calculer les nombres premiers de 1 \205 nmax\n");
        printf("\nlancez ce programme en faisant : crible nmax\n");
        pause();
        return 1;
    }
    nmax=atol(argv[1]);
    sprintf(limite,"%lu",nmax);      // on doit y retrouver argv[1]
    printf("\nCalcul des nombres premiers entre 1 et %lu\n", nmax);
    crible=(char *)malloc((1+(nmax-1)/16)*sizeof(char));
    if(strcmp(argv[1],limite)!=0 || crible==NULL) {
        printf("\nLa limite de l'intervalle de calcul est trop grande\n");
        pause();
        return 1;
    }
    clock_t start = clock();
    Eratosthene(nmax);
    printf("\nTemps elapse du calcul : %f s\n", (double)(clock()-start)/CLOCKS_PER_SEC);
    for(i=1; i<=nmax; i=i+1) if(prem(i)){np++; if(p1==0)p1=i; pn=i;} // calcul de np, p1 et pn
    printf("\nNombre de nombres premiers calcul\202s : %lu\n", np);
    printf("\nLe premier nombre premier : %lu\n",p1);
    printf("\nLe %lu-i\212me nombre premier : %lu\n", np, pn);
    free(crible);
    pause();
    return 0;
}