Criba de Eratóstenes

La criba de Eratóstenes (Sieve of Eratosthenes) es un algoritmo para encontrar todos los números primos en un segmento $[1;n]$ usando $O(n \log \log n)$ operaciones.

El algoritmo es muy simple: al comienzo escribimos todos los números entre 2 y $n$ . Marcamos todos los múltiplos propios de 2 (ya que 2 es el menor número primo) como compuestos. Un múltiplo propio de un número $x$ es un número mayor que $x$ y divisible por $x$ . Luego encontramos el siguiente número que no haya sido marcado como compuesto; en este caso es 3. Eso significa que 3 es primo, y marcamos todos los múltiplos propios de 3 como compuestos. El siguiente número no marcado es 5, que es el siguiente número primo, y marcamos todos sus múltiplos propios. Y continuamos este procedimiento hasta haber procesado todos los números de la fila.

En la siguiente imagen se puede ver una visualización del algoritmo para calcular todos los números primos en el rango $[1; 16]$ . Se observa que, con bastante frecuencia, marcamos números como compuestos varias veces.

La idea detrás es esta: Un número es primo si ninguno de los números primos más pequeños lo divide. Como iteramos sobre los números primos en orden, ya marcamos como divisibles todos los números que son divisibles por al menos uno de los números primos. Por lo tanto, si llegamos a una celda y no está marcada, entonces no es divisible por ningún número primo más pequeño y, en consecuencia, tiene que ser primo.

Implementación


int n;
vector<bool> is_prime(n+1, true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; i <= n; i++) {
    if (is_prime[i] && (long long)i * i <= n) {
        for (int j = i * i; j <= n; j += i)
            is_prime[j] = false;
    }
}

Este código primero marca todos los números excepto el cero y el uno como números primos potenciales, y luego comienza el proceso de cribado de los números compuestos. Para ello itera sobre todos los números de $2$ a $n$ . Si el número actual $i$ es un número primo, marca como compuestos todos los números que son múltiplos de $i$ , empezando desde $i^2$ . Esto ya es una optimización respecto de la forma ingenua de implementarlo, y está justificado porque todos los números más pequeños que son múltiplos de $i$ necesariamente también tienen un factor primo menor que $i$ , así que todos ellos ya fueron cribados antes. Como $i^2$ puede desbordar fácilmente el tipo int, se hace la verificación adicional usando el tipo long long antes del segundo bucle anidado.

Con esta implementación el algoritmo consume $O(n)$ de memoria (obviamente) y realiza $O(n \log \log n)$ operaciones (véase la sección siguiente).

Análisis asintótico

Es simple demostrar un tiempo de ejecución de $O(n \log n)$ sin saber nada sobre la distribución de los primos: ignorando la comprobación de is_prime, el bucle interno se ejecuta (a lo sumo) $n/i$ veces para $i = 2, 3, 4, \dots$ , lo que hace que el número total de operaciones en el bucle interno sea una suma armónica de la forma $n(1/2 + 1/3 + 1/4 + \cdots)$ , que está acotada por $O(n \log n)$ .

Demostremos que el tiempo de ejecución del algoritmo es $O(n \log \log n)$ . El algoritmo realizará $\frac{n}{p}$ operaciones por cada primo $p \le n$ en el bucle interno. Por lo tanto, necesitamos evaluar la siguiente expresión:

$\sum_{\substack{p \le n, \\ p \text{ prime}}} \frac n p = n \cdot \sum_{\substack{p \le n, \\ p \text{ prime}}} \frac 1 p.$

Recordemos dos hechos conocidos.

La cantidad de números primos menores o iguales que $n$ es aproximadamente $\frac n {\ln n}$ .
El $k$ -ésimo número primo es aproximadamente igual a $k \ln k$ (esto se sigue del hecho anterior).

Así podemos escribir la suma de la siguiente forma:

$\sum_{\substack{p \le n, \\ p \text{ prime}}} \frac 1 p \approx \frac 1 2 + \sum_{k = 2}^{\frac n {\ln n}} \frac 1 {k \ln k}.$

Acá extraímos el primer número primo 2 de la suma, porque $k = 1$ en la aproximación $k \ln k$ es $0$ y provoca una división por cero.

Ahora evaluemos esta suma usando la integral de la misma función sobre $k$ desde $2$ hasta $\frac n {\ln n}$ (podemos hacer esa aproximación porque, de hecho, la suma se relaciona con la integral como su aproximación por el método de los rectángulos):

$\sum_{k = 2}^{\frac n {\ln n}} \frac 1 {k \ln k} \approx \int_2^{\frac n {\ln n}} \frac 1 {k \ln k} dk.$

La antiderivada del integrando es $\ln \ln k$ . Usando una sustitución y eliminando términos de orden inferior, obtenemos el resultado:

$\int_2^{\frac n {\ln n}} \frac 1 {k \ln k} dk = \ln \ln \frac n {\ln n} - \ln \ln 2 = \ln(\ln n - \ln \ln n) - \ln \ln 2 \approx \ln \ln n.$

Ahora, volviendo a la suma original, obtenemos su evaluación aproximada:

$\sum_{\substack{p \le n, \\ p\ is\ prime}} \frac n p \approx n \ln \ln n + o(n).$

Se puede encontrar una demostración más rigurosa (que da una evaluación más precisa, exacta salvo factores constantes) en el libro de Hardy & Wright “An Introduction to the Theory of Numbers” (p. 349).

Distintas optimizaciones de la criba de Eratóstenes

La mayor debilidad del algoritmo es que “recorre” la memoria varias veces, manipulando solo elementos individuales. Esto no es muy amigable con la caché. Y por eso, la constante que queda oculta en $O(n \log \log n)$ es comparativamente grande.

Además, la memoria consumida es un cuello de botella para $n$ grandes.

Los métodos que se presentan abajo permiten reducir la cantidad de operaciones realizadas, y también acortar de forma notable la memoria consumida.

Cribado hasta la raíz

Obviamente, para encontrar todos los números primos hasta $n$ , basta con realizar el cribado solo por los números primos que no superan la raíz de $n$ .


int n;
vector<bool> is_prime(n+1, true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; i * i <= n; i++) {
    if (is_prime[i]) {
        for (int j = i * i; j <= n; j += i)
            is_prime[j] = false;
    }
}

Esa optimización no afecta la complejidad (en efecto, al repetir la demostración de arriba obtenemos la evaluación $n \ln \ln \sqrt n + o(n)$ , que es asintóticamente la misma según las propiedades de los logaritmos), aunque la cantidad de operaciones se reducirá de forma notable.

Cribado solo por los números impares

Como todos los números pares (excepto $2$ ) son compuestos, podemos dejar de revisar los números pares por completo. En su lugar, basta con operar solo con números impares.

Primero, esto nos permite reducir a la mitad la memoria necesaria. Segundo, reduce aproximadamente a la mitad la cantidad de operaciones que realiza el algoritmo.

Consumo de memoria y velocidad de las operaciones

Hay que notar que estas dos implementaciones de la criba de Eratóstenes usan $n$ bits de memoria al utilizar la estructura de datos vector<bool>. vector<bool> no es un contenedor regular que almacena una serie de bool (ya que en la mayoría de las arquitecturas de computadora un bool ocupa un byte de memoria). Es una especialización de vector<T> optimizada en memoria, que solo consume $\frac{N}{8}$ bytes de memoria.

Las arquitecturas de los procesadores modernos trabajan de forma mucho más eficiente con bytes que con bits, ya que normalmente no pueden acceder a los bits de forma directa. Así que por debajo vector<bool> almacena los bits en un bloque grande de memoria continua, accede a la memoria en bloques de unos pocos bytes, y extrae/asigna los bits con operaciones de bits como máscaras (bit masking) y desplazamientos (bit shifting).

Por eso hay cierto overhead al leer o escribir bits con un vector<bool>, y con bastante frecuencia usar un vector<char> (que usa 1 byte por cada entrada, o sea 8 veces la cantidad de memoria) es más rápido.

Sin embargo, para las implementaciones simples de la criba de Eratóstenes, usar un vector<bool> es más rápido. Estamos limitados por la velocidad con la que se pueden cargar los datos en la caché, y por lo tanto usar menos memoria da una gran ventaja. Un benchmark (enlace ) muestra que usar un vector<bool> es entre 1.4x y 1.7x más rápido que usar un vector<char>.

Las mismas consideraciones aplican también a bitset. Es también una forma eficiente de almacenar bits, similar a vector<bool>, así que solo ocupa $\frac{N}{8}$ bytes de memoria, pero es un poco más lenta al acceder a los elementos. En el benchmark de arriba bitset rinde un poco peor que vector<bool>. Otra desventaja de bitset es que hay que conocer el tamaño en tiempo de compilación.

Criba segmentada

De la optimización “cribado hasta la raíz” se sigue que no hace falta mantener todo el arreglo is_prime[1...n] en todo momento. Para cribar basta con conservar los números primos hasta la raíz de $n$ , es decir prime[1... sqrt(n)], partir el rango completo en bloques y cribar cada bloque por separado.

Sea $s$ una constante que determina el tamaño del bloque; entonces tenemos $\lceil {\frac n s} \rceil$ bloques en total, y el bloque $k$ ( $k = 0 … \lfloor {\frac n s} \rfloor$ ) contiene los números del segmento $[ks; ks + s - 1]$ . Podemos trabajar los bloques por turnos, es decir, para cada bloque $k$ recorremos todos los números primos (de $1$ a $\sqrt n$ ) y realizamos el cribado con ellos. Vale la pena notar que hay que modificar un poco la estrategia al manejar los primeros números: primero, todos los números primos de $[1; \sqrt n]$ no deberían eliminarse a sí mismos; y segundo, los números $0$ y $1$ deben marcarse como no primos. Al trabajar sobre el último bloque no hay que olvidar que el último número necesario $n$ no necesariamente está al final del bloque.

Como se discutió antes, la implementación típica de la criba de Eratóstenes está limitada por la velocidad con la que se pueden cargar los datos en las cachés de la CPU. Al partir el rango de números primos potenciales $[1; n]$ en bloques más chicos, nunca hace falta mantener varios bloques en memoria al mismo tiempo, y todas las operaciones son mucho más amigables con la caché. Como ya no estamos limitados por las velocidades de la caché, podemos reemplazar el vector<bool> por un vector<char> y ganar un poco más de rendimiento, ya que los procesadores pueden manejar lecturas y escrituras de bytes de forma directa y no necesitan apoyarse en operaciones de bits para extraer bits individuales. El benchmark (enlace ) muestra que, en esta situación, usar un vector<char> es unas 3x más rápido que usar un vector<bool>. Una advertencia: esos números pueden variar según la arquitectura, el compilador y los niveles de optimización.

Acá tenemos una implementación que cuenta la cantidad de primos menores o iguales que $n$ usando cribado por bloques.


int count_primes(int n) {
    const int S = 10000;
 
    vector<int> primes;
    int nsqrt = sqrt(n);
    vector<char> is_prime(nsqrt + 2, true);
    for (int i = 2; i <= nsqrt; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            for (int j = i * i; j <= nsqrt; j += i)
                is_prime[j] = false;
        }
    }
 
    int result = 0;
    vector<char> block(S);
    for (int k = 0; k * S <= n; k++) {
        fill(block.begin(), block.end(), true);
        int start = k * S;
        for (int p : primes) {
            int start_idx = (start + p - 1) / p;
            int j = max(start_idx, p) * p - start;
            for (; j < S; j += p)
                block[j] = false;
        }
        if (k == 0)
            block[0] = block[1] = false;
        for (int i = 0; i < S && start + i <= n; i++) {
            if (block[i])
                result++;
        }
    }
    return result;
}

El tiempo de ejecución del cribado por bloques es el mismo que el de la criba de Eratóstenes regular (salvo que el tamaño de los bloques sea muy chico), pero la memoria necesaria se reduce a $O(\sqrt{n} + S)$ y obtenemos mejores resultados de caché. Por otro lado, habrá una división por cada par de un bloque y un número primo de $[1; \sqrt{n}]$ , y eso será bastante peor para tamaños de bloque más chicos. Por lo tanto, es necesario mantener un equilibrio al elegir la constante $S$ . Obtuvimos los mejores resultados para tamaños de bloque entre $10^4$ y $10^5$ .

Encontrar primos en un rango

A veces necesitamos encontrar todos los números primos en un rango $[L,R]$ de tamaño chico (p. ej. $R - L + 1 \approx 1e7$ ), donde $R$ puede ser muy grande (p. ej. $1e12$ ).

Para resolver un problema así, podemos usar la idea de la criba segmentada. Pregeneramos todos los números primos hasta $\sqrt R$ , y usamos esos primos para marcar todos los números compuestos en el segmento $[L, R]$ .


vector<char> segmentedSieve(long long L, long long R) {
    // generar todos los primos hasta sqrt(R)
    long long lim = sqrt(R);
    vector<char> mark(lim + 1, false);
    vector<long long> primes;
    for (long long i = 2; i <= lim; ++i) {
        if (!mark[i]) {
            primes.emplace_back(i);
            for (long long j = i * i; j <= lim; j += i)
                mark[j] = true;
        }
    }
 
    vector<char> isPrime(R - L + 1, true);
    for (long long i : primes)
        for (long long j = max(i * i, (L + i - 1) / i * i); j <= R; j += i)
            isPrime[j - L] = false;
    if (L == 1)
        isPrime[0] = false;
    return isPrime;
}

La complejidad temporal de este enfoque es $O((R - L + 1) \log \log (R) + \sqrt R \log \log \sqrt R)$ .

También es posible no pregenerar todos los números primos:


vector<char> segmentedSieveNoPreGen(long long L, long long R) {
    vector<char> isPrime(R - L + 1, true);
    long long lim = sqrt(R);
    for (long long i = 2; i <= lim; ++i)
        for (long long j = max(i * i, (L + i - 1) / i * i); j <= R; j += i)
            isPrime[j - L] = false;
    if (L == 1)
        isPrime[0] = false;
    return isPrime;
}

Obviamente, la complejidad es peor: $O((R - L + 1) \log (R) + \sqrt R)$ . Sin embargo, en la práctica sigue ejecutándose muy rápido.

Modificación de tiempo lineal

Podemos modificar el algoritmo de forma tal que solo tenga complejidad temporal lineal. Este enfoque se describe en el artículo Criba lineal. Sin embargo, este algoritmo también tiene sus propias debilidades.