Buscando números primos, con la criba de Eratóstenes

[Pantaláimon]

No me acordaba que había participado aquí. Con lo que me gustan las frikadas numéricas. Aquí una vuelta de tuerca al algoritmo sin perder la esencia (ansi C):

Código C:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
 
int main(void) {
   unsigned top = 1000000000;
   unsigned rtop = sqrt(top);
   unsigned char* nums = calloc( top, sizeof *nums);
   unsigned primsTop = 2.2 * (top / log(top)); /* apaño para ahorrar memoria */
   unsigned* prims = calloc( primsTop, sizeof *prims);
   unsigned n; /* contador primos */
   unsigned i, i2, ii;
   unsigned j, k, l, x;
 
   k = x = nums[4] = nums[6] = nums[8] = 1;
   i = l = n = prims[0] = 2;
   prims[1] = 3;
   l = 5;
   i2 = i + i;
   ii = i * i;
 
   /* calculo de la tabla */
   for (i = 3; i <= rtop; i += 2, i2 += 4) {
      /* El primer primo de la tabla es el primero con valor 0 */
      if(nums[i]) continue;
 
      /* todo numero j tal que l <= j <= i*i-1 con nums[j] == 0 es primo */
      ii = i*i;
      for (j = l; j < ii; j += 2 )
         if(nums[j] == 0)
            prims[n++] = j;
      l = ii+2;
      
      /* Dado k que va de x a n-1, se descartan los compuestos i*prims[k]    *
       * donde se ve que i == prims[x]                                       *
       * (Es una forma de evitar descartar números que ya están descartados) */
      for (k = x; k < n; ++k) {
         j = i*prims[k];
         if (j >= top) break;
         nums[j] = 1;
      }
     
      /* a partir de aquí, se descartan los compuestos de la manera común */
      if (j < top) j += i2;
      for (; j < top; j += i2) {
         if(!nums[j])
            nums[j] = 1;
      }
      ++x;
   }
 
   /* calcular resto de primos */
   for (i = prims[n-1]+2; i < top; i += 2) {
      if (nums[i] == 0) {
         ++n;
         /* prims[n++] = i; */
      }
   }
   
   printf("%d primos\n", n);
   free(nums);
   free(prims);
 
   return 0;
}
 

Un saludo!