Nos ponemos serios

Tras el experimento anterior, en el que se usaba SQL y Perl para calcular números primos, ahora probaremos un método algo más profesional y bastante más rápido.

Utilizaremos la librería GMP, GNU Multiple Precision, un conjunto de funciones para operar con números enteros y flotantes de cualquier tamaño, donde muchas de las operaciones están implementadas directamente en ensamblador optimizado para varios tipos de procesadores.

Esta librería incorpora la función int mpz_probab_prime_p (mpz_t n, int reps), una función que ejecuta el test de primalidad de Miller-Rabin sobre un número dado. Es decir, en este caso ni aplicamos una criba ni dividimos entre todos los primos anteriores, aunque no tendremos la absoluta certeza de que todos los números obtenidos sean primos.

Este código me ha servido de base para el siguiente:

primes.cpp

[cpp]
#include <stdio.h>
#include <time.h>
#include <gmp.h>

// Using GMP library (http://gmplib.org/manual/)
// similar: http://stackoverflow.com/questions/622/most-efficient-code-for-the-first-10000-prime-numbers
void calculatePrimes ( mpz_t N_max ) {
time_t start,end;
double elapsed;
mpz_t prime, N_primes;
mpz_init(prime);
mpz_init(N_primes);
mpz_set_ui(prime, 1);
mpz_set_ui(N_primes, 0);

time(&start);

while( mpz_cmp(prime, N_max) < 0 ) {
mpz_nextprime(prime, prime);
//printf("%s, ", mpz_get_str(NULL,10,prime));
mpz_add_ui(N_primes, N_primes, 1);
}
// skip the last prime, which is over the range
mpz_sub_ui(N_primes, N_primes, 1);

time(&end);
elapsed = difftime(end, start);

/*printf("%s primes found below %s in %f seconds (last one found over the limit was %s)\n",
mpz_get_str(NULL,10,N_primes),
mpz_get_str(NULL,10,N_max),
elapsed,
mpz_get_str(NULL,10,prime));*/
printf("%s %s %.0f\n",mpz_get_str(NULL,10,N_max),mpz_get_str(NULL,10,N_primes),elapsed);
}

int main() {
int i;
mpz_t N_max;
mpz_init(N_max);
mpz_set_ui(N_max, 100);
for(i=0; i<8; i++) {
calculatePrimes(N_max);
mpz_mul_ui(N_max, N_max, 10);
}
}

[/cpp]

Se compila en Linux con:
g++ -lgmp primes.cpp -o primes

Los resultados obtenidos son algo así:

Máximo número Primos Duración
100 25 0s
1.000 168 0s
10.000 1.229 0s
100.000 9.592 0s
1.000.000 78.498 1s
10.000.000 664.579 12s
100.000.000 5.761.455 126s (~2min)
1.000.000.000 50.847.534 1.275s (~21min)

En el algoritmo en Perl, calcular todos los número primos por debajo de 100 millones llevó 5.929 segundos (1 hora y 39 minutos). En este caso, gracias a GMP, se ha tardado 126 segundos (poco más de 2 minutos) para el mismo rango. Es decir, se consigue lo mismo 47 veces más rápido.

No está mal. Apretaremos aún más las tuercas.