Criba de Eratóstenes

[Estoy de vacaciones.... y claro, me aburro. Lo mío es ir haciendo cositas, y ahora me encuentro con demasiadas horas por delante vacías]

Bueno, hoy quiero hablar de la Criba de Eratóstenes. Más información en su wiki. De manera resumida es algo así:
Explicación
La criba de Eratóstenes consiste en calcular los números primos menores que N. Este algoritmo es muy eficiente para calcular números primos pequeños ( menores de 10M ).
Para ello utilizaremos un vector de booleanos. Si una celda vale 1, es primo, si vale cero, no lo es. Bien. Al principio ponemos a 1 todos el vector, de modo que, en todos los números sean "primos". Ahora empezamos con el 2 ( ya que todos los números son divisibles entre 1 ). Por cada numero divisible entre 2 menor que N, lo tachamos (4, 6, 8, 10, .... ). Después hacemos lo mismo con el número 3 (6, 9 12...) . Cuando llegamos al número 4, vemos que está tachado ya que es múltiplo de 2, entonces lo ignoramos. Llegamos al 5, y tachamos todos los múltiplos de 5 menores que N. Y así... si llegamos a un número tachado, pasamos al siguiente. Si llegamos a un número que no está tachado, vamos tachando todos sus múltiplos. Al final, los números NO tachados ( con un 1 ), son los números primos. Vean su wiki para ver la representación visual.

Qué he hecho
Bien, pues he hecho un código en C ( estoy aprendiendo a programar ) para calcular los números primos entre 2 y N, siendo N=10.000.000 y que guarda los resultados en un fichero. En mi PC tarda aproximadamente 1.45 segundos. Me he complicado un poco la vida y lo he que uno pueda . Al final he conseguido poner el valor de N tan grande como quiera ( en realidad no es así, ya que estamos limitados por la memória RAM del sistema )calcular los números primos entre 2 y 450.000.000. El fichero resultante pesa unos 400MB.

Compilación
Para compilar el código fuente, lo grabamos en un fichero "criba.c" y lo compilamos así:
$ gcc criba.c -o criba -lm
Y para ejecutar el código tan solo tenemos que escribir:
$ ./criba
Esto nos creará un ficherito "primos.txt" con los resultados.

Configuración
Pues básicamente jugaremos con las constantes N y M. N es el tamaño del bloque y M en número de bloques. Osea, que se van a calcular los números primos que vayan desde 2 hasta N*M. En mi PC N=1M funciona, pero en otro PC igual podría dar fallo. Disminuyese esta cantidad y auméntese en numero de iteraciones ( M ). Recuerden que el algoritmo se vuelve ineficiente a partir de 10M!


Codigo

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// Each size block
#define N 1000000
// Iterations
#define M 10
//Memory that will be used for store all primes( must be aviable )
#define memory_size M*N*sizeof(long int)*0.08

// Put all values to zero
void initialization(unsigned long int * primes_numbers[]){
int i;
for (i=0;i<memory_size;i++)
(*primes_numbers)[i] = 0;
}

// Save all primes calculated in a file
void save_in_file(unsigned long int * primes_numbers[]){
FILE * fp;
unsigned long long int i;
fp = fopen("primos.txt", "w");
for (i=0;(*primes_numbers)[i];i++)
fprintf(fp, "%ld\n", (*primes_numbers)[i]);
fclose(fp);
}

// Count all the primes calculated
void count(unsigned long int * primes_numbers[]){
unsigned long long int i,j;
j = 0;
for (i=0;(* primes_numbers)[i];i++)
j++;

printf("Total primos encontrados: %lld\n", j);

}

// Check that in the list of prime numbers all are primes
// alert: very slow in very large amount of numbers
void check(unsigned long int * primes_numbers[]){
unsigned long long int i;
for (i=1; (*primes_numbers)[i] != 0; i++){
if ( is_prime((*primes_numbers)[i]) != 1)
printf("FUckkk! %ld\n", (*primes_numbers)[i] );
}
}

// Print all primes calculated
void imprime(unsigned long int * datos[]){
unsigned long long int i;
printf("Llistat:\n");
for (i=0;(*datos)[i]!=0;i++)
printf("%ld\n", (*datos)[i]);
}



int sieve_of_eratosthenes(unsigned long int * primes[], long long int * last_prime_inserted){
char list[N];
unsigned long long int i, j, max;
list[0] = 0;
for (i=2;i<N;i++)
list[i] = 1;

max = sqrt(N)+1;

for (i=2;i<max;i++){
if (list[i])
for (j=2; i*j<N;j++)
list[i*j] = 0;
}
j = 0;
for (i=0;i<N;i++)
if (list[i]){
(*primes)[j] = i;
j++;
}
*last_prime_inserted = j - 1;
return 0;
}

// Calculate the next block of primes
int next(unsigned long int * prime[], unsigned long long int block_number, long long int * last_prime_inserted){
char list[N];
unsigned long long int i, j, min, max;
min = (block_number-1)*N;
max = block_number*N;

//Set all numbers as primes
for (i=0;i<N;i++)
list[i] = 1;
list[0] = 0;

// Mark all numbers alredy calculated
for (i=0;(*prime)[i];i++){
for (j=min/(*prime)[i]; (*prime)[i]*j<max; j++){
if ((*prime)[i]*j>min){
list[((*prime)[i]*j)%N] = 0;
}
}
}

j = *last_prime_inserted;
// Save the new primes
for (i=0;i<N;i++){
if (list[i]){
// printf("%ld j:%ld\n", i+min, j);
(*prime)[j] = i+min;
j++;
}
}
*last_prime_inserted = j - 1;
}

// Calculate all the primes
int calc(int n){
unsigned long long int position_last_prime;
unsigned long int * primes_numbers;
primes_numbers = malloc(memory_size * sizeof(unsigned long int));
initialization(&primes_numbers);
sieve_of_eratosthenes(&primes_numbers, &position_last_prime);
int j;
for (j=2; j<=n; j++){
printf("Calculando %d...\n", j);
next(&primes_numbers, j, &position_last_prime);
}
//imprime(numeros_primos);
printf("Guardando...\n");
save_in_file(&primes_numbers);
//check(&numeros_primos);
count(&primes_numbers);
printf("Done!\n");
free(primes_numbers);
return 0;
}

int main(void){
calc(M);
return 0;
}

int is_prime(int n){
int i, j, prime;
prime = 1;
for (j=2; j<n/2; j++)
if (n % j == 0){
prime = 0;
break;
}
return prime;
}

Fin Codigo

Pues esto es todo, señores. Espero no haberles aburrido demasiado.