Formule de Lehmer : calcul de pi(x)

Télécharger le projet
Commenter

Description

Pour explorer les nombres premiers il y a deux fonctions basiques, le calcul du n-ième nombre premier et le calcul de pi(x) : le nombre de nombres premiers compris entre 1 et x. On s'intéresse ici à la fonction pi(x). Une méthode très simple pour calculer pi(x) consiste à effectuer un crible d'Ératosthène puis à compter les nombres premiers entre 1 et x.

Après la méthode de Legendre, d'autres méthodes sont apparues pour en faire le calcul plus efficacement. La méthode de D. H. Lehmer, 1905-1991, est l'une des plus anciennes parmi celles-ci. Le programme présenté ici est une seconde version, assez simple, de la méthode de Lehmer pour mieux comprendre sa formulation qui n'est pas immédiatement évidente. En prenant diverses précautions pour avoir de bonnes performances, cela permet de calculer pi(10^12) en 9 secondes, pi(10^13) en 80 secondes, pi(10^14) en 12 minutes et pi(10^15) en 2 heures et 8 minutes. A noter que 10^15 c'est un million de gigas et qu'à ce niveau il n'est plus question de calculer, d'archiver ou de compter les nombres premiers un par un. Les méthodes de Legendre et de Meissel y sont disponibles pour les comparer à celle de Lehmer. Attention, la récurrence de phi(x,a) utilise ici la récursivité, ce qui peut provoquer un problème d'exécution pour un calcul de pi(x) avec x très grand. Sous Windows le programme fonctionne en mode 64 bits.

Voir : http://fr.wikipedia.org/wiki/Fonction_de_compte_des_nombres_premiers pour plus d'information.

Dans une variante on a mesuré le temps d'exécution de chaque partie de la méthode de Lehmer. Pour un exemple de calcul de pi(10^14) en 739 secondes, l'instruction phi(x,a) utilise 581 secondes et l'instruction erato(w) utilise 146 secondes. Le reste de la fonction lehmer(x) est relativement beaucoup plus rapide.

Pour obtenir une version encore plus efficace on peut essayer d'utiliser des entiers de grande taille quand c'est nécessaire, faire du multithreading et limiter la récurrence de la fonction phi(x,a) avant a=4. On pourrait aussi utiliser la récurrence de phi(x,a) sans faire de récursivité. On trouve sur Internet les sources de Kim Walisch ou de Dana Jacobsen qui non seulement ont des performances excellentes mais qui développent aussi plusieurs autres méthodes dont celle de Lagarias, Miller et Odlyzko et celles de Deléglise et Rivat ou de Deléglise seul qui sont plus récentes.

L'URL http://mathworld.wolfram.com/PrimeCountingFunction.html permet vérifier les calculs obtenus.
Voir aussi les fichiers ci-joints : aide.txt, Legendre.txt, Meissel.txt, Lehmer.txt, Lehmer.pdf, lehmer.jpg et pi(x).txt. 

/*************************************************
*                                                *
*  Calcul de la fonction pi(x) -  @author pgl10  *
*                                                *
*     On utilise la méthode de D. K. Lehmer      *
*                                                *
*    pi(x) est le nombre de nombres premiers     *
*                                                *
*           compris entre 1 et x                 *
*                                                *
*************************************************/

#include <iostream>
#include <stdint.h>
#include <vector>
#include <ctime>

namespace {
// p[] : les premiers nombres premiers
   std::vector<int64_t> p;
// crible[] : le crible d'Ératosthène
   char* crible;
// nprems[i] : nombre de premiers impairs 
// compris entre 1 et 128*(i+1)-1
   int64_t* nprems;
   long double eps = 0.0000000000001;
   long double one = 1.0 + eps;
}

int64_t isqrt(int64_t n) {
// Partie entière de la racine carrée de l'entier n
// if(n < 0) return 0;
   long double x = sqrt((long double)n);
   return int64_t(x*one);
}

int64_t iroot(int64_t n, int64_t k) {
// Partie entière de la racine k-ième de l'entier n
// if(n < 0) return 0;
   long double t = 1.0 / (long double)k;
   long double x = pow((long double)n, t);
   return int64_t(x*one);
}

bool is_prime(int64_t n) {
// L'entier impair n est-il premier ?
// if(n<2) return false;
// if(n%2==0) return n==2;
   if(n%3==0) return n==3;
   int64_t d=5;
   while(d*d<=n) {
      if(n%d) d+=2; else return false;
      if(n%d) d+=4; else return false;
   }
   return true;   
}

void primes(int64_t n) {
// Calcul des n premiers nombres premiers
// p[i] pour i = 1 à n ( p[1] = 2 )
// if(n<=0) return;
   p.clear();
   p.push_back(0);
   p.push_back(2);
   if(n==1) return;
   p.push_back(3);
   if(n==2) return;
   int64_t i=1, k=2, s=4;
   while(k<n) {
      i=i+s;
      if(is_prime(i)) {p.push_back(i); k++;}
      s=6-s;
   }
}

int64_t pi(int64_t x) {
// Calcul de pi(x)
// utilisé seulement quand x est petit
   int pp[24] = {0, 0, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 
                 5, 6, 6, 6, 6, 7, 7, 8, 8, 8, 8, 9};
   if(x<24) return pp[x];
   int64_t i, k, m;
   m=x-1;
   k=9;
   for(i=29; i<=m; i=i+2) if(is_prime(i)) k++;
   if((x%2>0) && is_prime(x)) k++;
   return k;
}

void erato(int64_t n) {
// Crible d'Ératosthène ( de 1 à n )
// crible[] contient les bits utiles du crible.
// le crible ne concerne que les nombres impairs :
// crible[0] pour 1 à 15, crible[1] pour 17 à 31, etc.
   char c0[8]={'xFE', 'xFD', 'xFB', 'xF7', 'xEF', 'xDF', 'xBF', 'x7F'};
   char c1[8]={'x01', 'x02', 'x04', 'x08', 'x10', 'x20', 'x40', 'x80'};
   int64_t i, i2, j, s, t;
   s=1+(n-1)/16;
   crible = new char[s];
   if(crible==NULL) {std::cerr << "Erreur" << std::endl; return;}
   // Initialisation de crible[] : tous les marqueurs à 1
   memset(crible, 'xFF', s);
   // 1 n'est pas premier
   crible[0]='xFE';
   // On marque les multiples impairs de 3
   for(j=9; j<=n; j=j+6) crible[j>>4]=crible[j>>4]&c0[(j/2)&7];
   i=5;
   i2=25;
   s=2;
   // On marque les multiples de i non multiples de 2 et 3
   while(i2<=n) {
      if((crible[i>>4]&c1[(i/2)&7]) != 'x00')
         for(j=i2, t=s; j<=n; j=j+t*i, t=6-t)
            crible[j>>4]=crible[j>>4]&c0[(j/2)&7];
      i=i+s;
      i2=i*i;
      s=6-s;
   }
   // pour terminer on va compter le nombre de nombres premiers impairs
   // compris entre 1 et 128*(j+1)-1 pour j = 0, 1, ... , n/128
   // ce qui correspond à grouper les octets de crible[] par 8
   int64_t m = 1+n/128;
   nprems = new int64_t[m];
   if(nprems==NULL) {std::cerr << "Erreur" << std::endl; return;}
   int64_t a = 0;
   for(int64_t j=0; j<m; j++) {
      nprems[j] = a;
      int64_t im = 128*(j+1);
      for(int64_t i=128*j+1; i<im; i+=2)
         if((crible[i>>4]&c1[(i/2)&7]) != 'x00') nprems[j]++;
      a = nprems[j];
   }
}

int64_t px(int64_t x) {
// Calcul de pi(x)
// utilisé seulement quand x est petit ou moyen
// après avoir effectué un crible d'Ératosthène de taille suffisante
   char c1[8]={'x01', 'x02', 'x04', 'x08', 'x10', 'x20', 'x40', 'x80'};
   if(x == 2) return 1;
   int64_t k = 1;
   int64_t x1 = x/128;
   if(x1 > 0) k = k + nprems[x1-1];
   // dans le groupe de 8 octets qui a le marqueur de x
   for(int64_t i=128*x1+1; i<=x; i+=2)
      if((crible[i>>4]&c1[(i/2)&7]) != 'x00') k++;
   return k;
}

int64_t phi(int64_t x, int64_t a) {
// Calcul de phi(x,a)
   if(x == 0) return 0;
   if(a == 4) return (x+1)/2-(x/3+1)/2-(x/5+1)/2+(x/15+1)/2-
                     (x/7+1)/2+(x/21+1)/2+(x/35+1)/2-(x/105+1)/2;
   if(a == 3) return (x+1)/2-(x/3+1)/2-(x/5+1)/2+(x/15+1)/2;
   if(a == 2) return (x+1)/2-(x/3+1)/2;
   if(a == 1) return (x+1)/2;
   if(a > 2+(x-4)/3) return 1;
   return phi(x, a-1) - phi(x/p[a], a-1);
}

int64_t legendre(int64_t x) {
// Calcul de pi(x) : méthode de Legendre
   if(x < 80) return pi(x);
   int64_t a = pi(isqrt(x));
   int64_t ps = (int64_t)p.size();
   bool done = false;
   if(a > ps-1) {
      primes(a); 
      done = true;
   }
   int64_t t = phi(x, a) + a - 1;
   if(done) if(ps) primes(ps);
   return t;
}

int64_t meissel(int64_t x) {
// Calcul de pi(x) : méthode de Meissel
   if(x < 80) return pi(x);
   int64_t b = pi(isqrt(x));
   int64_t c = pi(iroot(x, 3));
   int64_t ps = (int64_t)p.size();
   bool done = false;
   if(b > ps-1) {
      primes(b); 
      done = true;
   }
   int64_t s = phi(x,c) + (b+c-2) * (b-c+1) / 2;
   for(int64_t i=c+1; i<=b; i++) {
      int64_t w = x / p[i];
      if(i == c+1) erato(w);
      s = s - px(w);
   }
   delete [] crible;
   delete [] nprems;
   if(done) if(ps) primes(ps);
   return s;
}

int64_t lehmer(int64_t x) {
// Calcul de pi(x) : méthode de Lehmer
   if(x < 80) return pi(x);
   int64_t a = pi(iroot(x, 4));
   int64_t b = pi(iroot(x, 2));
   int64_t c = pi(iroot(x, 3));
   int64_t ps = (int64_t)p.size();
   bool done = false;
   if(b > ps-1) {
      primes(b); 
      done = true;
   }
   int64_t s = phi(x,a) + (b+a-2) * (b-a+1) / 2;
   for(int64_t i=a+1; i<=b; i++) {
      int64_t w = x / p[i];
      if(i == a+1) erato(w);
      s = s - px(w);
      int64_t m = px(isqrt(w));
      if(i <= c)
         for(int64_t j=i; j<=m; j++) 
            s = s - px(w/p[j]) + j - 1;
   }
   delete [] crible;
   delete [] nprems;
   if(done) if(ps) primes(ps);
   return s;
}

int main(int argc, char *argv[]) {
   if(argc != 2) {
      std::cout << "Il faut faire : lehmer n" << std::endl;
      return 0;
   }
   int64_t n = 0;
   for(int i=0; argv[1][i]!=''; i++) {
      if((int(argv[1][i])<48)||(int(argv[1][i])>57)) {
         std::cout << "Commande invalide" << std::endl;
         return 0;
      }
      n = n * 10 + int(argv[1][i]) - 48;
   }
   clock_t start = clock();
   std::cout << "pi(" << n << ") = " << lehmer(n);
   std::cout << "   en " << double(clock()-start)/double(CLOCKS_PER_SEC) << " s" << std::endl;
   return 0;
}

Conclusion :

Pour calculer pi(10^15) la fonction lehmer(x) utilise déjà le fonctionnement de la mémoire virtuelle de Windows, ce qui rend très difficile d'aller plus loin. Y a-t-il des améliorations simples qui permettraient d'aller plus vite ou plus loin assez facilement ? Vos commentaires sont les bienvenus. Merci.
Télécharger le projet

Codes Sources

A voir également

Du même auteur (pgl10)