El algoritmo más conocido para encontrar todos los números primos no mayores que uno dado es el tamiz de Eratóstenes . Funciona muy bien para números de hasta miles de millones, tal vez hasta decenas de miles de millones, si se escriben cuidadosamente. Sin embargo, todos los que les gusta divertirse con los números primos saben que siempre quieres tener tantos como puedas. Una vez, para resolver un problema en un hackerrank, necesitaba una base de datos en memoria de primos de hasta cien mil millones. Con la optimización de memoria máxima, si representa números impares en un tamiz de Eratóstenes con una matriz de bits, su tamaño será de aproximadamente 6 gigabytes, lo que no cabe en la memoria de mi computadora portátil. Hay una modificación del algoritmo que es mucho menos exigente en la memoria (dividiendo el rango original de números en varias piezas y procesando una pieza a la vez):Tamiz segmentado de Eratóstenes , pero es más difícil de implementar, y el resultado completo no cabe en la memoria de todos modos. A continuación, traigo a su atención un algoritmo que es casi tan simple como el tamiz de Eratóstenes, pero que ofrece una doble optimización por memoria (es decir, una base de datos de primos de hasta cien mil millones ocupará aproximadamente 3 gigabytes, que ya deberían caber en la memoria de una computadora portátil estándar).Primero, recordemos cómo funciona el tamiz de Eratóstenes: en una serie de números del 1 al N, los números divisibles por dos se tachan, luego los números divisibles por 3, y así sucesivamente. Los números restantes después de tachar serán simples:2 3 . 5 . 7 . 9 . 11 . 13 . 15 . 17 . 19 . 21 . 23 . 25 . 27 . 29 ...
2 3 . 5 . 7 . . . 11 . 13 . . . 17 . 19 . . . 23 . 25 . . . 29 ...
2 3 . 5 . 7 . . . 11 . 13 . . . 17 . 19 . . . 23 . . . . . 29 ...
Código C #class EratospheneSieve
{
private bool[] Data;
public EratospheneSieve(int length)
{
Data = new bool[length];
for (int i = 2; i < Data.Length; i++) Data[i] = true;
for (int i = 2; i * i < Length; i++)
{
if (!Data[i]) continue;
for (int d = i * i; d < Length; d += i) Data[d] = false;
}
}
public int Length => Data.Length;
public bool IsPrime(int n) => Data[n];
}
En mi computadora portátil, cuenta hasta 15 mil millones de segundos y consume un gigabyte de memoria.El algoritmo que usaremos se llama Wheel Factorization . Para comprender su esencia, escribimos números naturales con una tableta de 30 piezas seguidas: 0 1 2 3 4 5 ... 26 27 28 29
30 31 32 33 34 35 ... 56 57 58 59
60 61 62 63 64 65 ... 86 87 88 89
...
Y tache todos los números divisibles por 2, 3 y 5:. 1 . . . . . 7 . . . 11 . 13 . . . 17 . 19 . . . 23 . . . . . 29
. 31 . . . . . 37 . . . 41 . 43 . . . 47 . 49 . . . 53 . . . . . 59
. 61 . . . . . 67 . . . 71 . 73 . . . 77 . 79 . . . 83 . . . . . 89
...
Se puede ver que en cada fila los números se tachan en las mismas posiciones. Como tachamos solo números compuestos (con la excepción, de hecho, de los números 2, 3, 5), los números primos solo pueden estar en posiciones sin marcar.Como hay ocho números de este tipo en cada fila, esto nos da la idea de representar cada treinta números con un byte, cada uno de los ocho bits indicará la simplicidad del número en la posición correspondiente. ¡Esta coincidencia no solo es estéticamente hermosa, sino que también simplifica enormemente la implementación del algoritmo!
Implementación técnica
En primer lugar, nuestro objetivo es llegar a cien mil millones, por lo que ya no se puede prescindir de los números de cuatro bytes. Por lo tanto, nuestro algoritmo se verá así:class Wheel235
{
private byte[] Data;
public Wheel235(long length)
{
...
}
public bool IsPrime(long n)
{
...
}
...
}
Para facilitar la implementación, redondearemos la longitud al número divisible más cercano por 30. Luegoclass Wheel235
{
private byte[] Data;
public Wheel235(long length)
{
length = (length + 29) / 30 * 30;
Data = new byte[length/30];
...
}
public long Length => Data.Length * 30L;
...
}
A continuación, necesitamos dos matrices: una traduce el índice 0..29 en el número de bit, la segunda, por el contrario, el número de bit en la posición treinta. Usando estas matrices y un poco de aritmética de bits, construimos una correspondencia uno a uno entre números naturales y bits en la matriz de datos. Para simplificar, usamos solo dos métodos: IsPrime (n), que le permite averiguar si el número es primo, y ClearPrime (n), que establece el bit correspondiente a 0, marcando el número n como compuesto:class Wheel235
{
private static int[] BIT_TO_INDEX = new int[] { 1, 7, 11, 13, 17, 19, 23, 29 };
private static int[] INDEX_TO_BIT = new int[] {
-1, 0,
-1, -1, -1, -1, -1, 1,
-1, -1, -1, 2,
-1, 3,
-1, -1, -1, 4,
-1, 5,
-1, -1, -1, 6,
-1, -1, -1, -1, -1, 7,
};
private byte[] Data;
...
public bool IsPrime(long n)
{
if (n <= 5) return n == 2 || n == 3 || n == 5;
int bit = INDEX_TO_BIT[n % 30];
if (bit < 0) return false;
return (Data[n / 30] & (1 << bit)) != 0;
}
private void ClearPrime(long n)
{
int bit = INDEX_TO_BIT[n % 30];
if (bit < 0) return;
Data[n / 30] &= (byte)~(1 << bit);
}
}
Tenga en cuenta que si el número es divisible por 2, 3 o 5, entonces el número de bit será -1, lo que significa que el número es obviamente compuesto. Bueno y, por supuesto, estamos procesando un caso especial n <= 5.Ahora, de hecho, el algoritmo en sí, que casi literalmente repite el tamiz de Eratóstenes: public Wheel235(long length)
{
length = (length + 29) / 30 * 30;
Data = new byte[length / 30];
for (long i = 0; i < Data.Length; i++) Data[i] = byte.MaxValue;
for (long i = 7; i * i < Length; i++)
{
if (!IsPrime(i)) continue;
for (long d = i * i; d < Length; d += i) ClearPrime(d);
}
}
En consecuencia, toda la clase resultante:tiene este aspectoclass Wheel235
{
private static int[] BIT_TO_INDEX = new int[] { 1, 7, 11, 13, 17, 19, 23, 29 };
private static int[] INDEX_TO_BIT = new int[] {
-1, 0,
-1, -1, -1, -1, -1, 1,
-1, -1, -1, 2,
-1, 3,
-1, -1, -1, 4,
-1, 5,
-1, -1, -1, 6,
-1, -1, -1, -1, -1, 7,
};
private byte[] Data;
public Wheel235(long length)
{
length = (length + 29) / 30 * 30;
Data = new byte[length / 30];
for (long i = 0; i < Data.Length; i++) Data[i] = byte.MaxValue;
for (long i = 7; i * i < Length; i++)
{
if (!IsPrime(i)) continue;
for (long d = i * i; d < Length; d += i) ClearPrime(d);
}
}
public long Length => Data.Length * 30L;
public bool IsPrime(long n)
{
if (n >= Length) throw new ArgumentException("Number too big");
if (n <= 5) return n == 2 || n == 3 || n == 5;
int bit = INDEX_TO_BIT[n % 30];
if (bit < 0) return false;
return (Data[n / 30] & (1 << bit)) != 0;
}
private void ClearPrime(long n)
{
int bit = INDEX_TO_BIT[n % 30];
if (bit < 0) return;
Data[n / 30] &= (byte)~(1 << bit);
}
}
Solo hay un pequeño problema: ¡este código no funciona para números mayores de aproximadamente 65 mil millones! Cuando intenta iniciarlo con dichos números, el programa se bloquea con un error:Unhandled Exception: System.OverflowException: Array dimensions exceeded supported range.
El problema es que en las matrices de C # no pueden tener más de 2 ^ 31 elementos (aparentemente porque la implementación interna utiliza un índice de matriz de cuatro bytes). Una opción es crear, por ejemplo, una matriz larga [] en lugar de una matriz byte [], pero esto complicará un poco la aritmética de bits. Por simplicidad, iremos hacia el otro lado, creando una clase especial que simule la matriz deseada, sosteniendo dos matrices cortas en su interior. Al mismo tiempo, le daremos la oportunidad de guardarnos en el disco para que pueda calcular los números primos una vez y luego usar la base de datos nuevamente: class LongArray
{
const long MAX_CHUNK_LENGTH = 2_000_000_000L;
private byte[] firstBuffer;
private byte[] secondBuffer;
public LongArray(long length)
{
firstBuffer = new byte[Math.Min(MAX_CHUNK_LENGTH, length)];
if (length > MAX_CHUNK_LENGTH) secondBuffer = new byte[length - MAX_CHUNK_LENGTH];
}
public LongArray(string file)
{
...
}
public long Length => firstBuffer.LongLength + (secondBuffer == null ? 0 : secondBuffer.LongLength);
public byte this[long index]
{
get => index < MAX_CHUNK_LENGTH ? firstBuffer[index] : secondBuffer[index - MAX_CHUNK_LENGTH];
set
{
if (index < MAX_CHUNK_LENGTH) firstBuffer[index] = value;
else secondBuffer[index - MAX_CHUNK_LENGTH] = value;
}
}
public void Save(string file)
{
...
}
}
Finalmente, combine todos los pasos enversión final del algoritmoclass Wheel235
{
class LongArray
{
const long MAX_CHUNK_LENGTH = 2_000_000_000L;
private byte[] firstBuffer;
private byte[] secondBuffer;
public LongArray(long length)
{
firstBuffer = new byte[Math.Min(MAX_CHUNK_LENGTH, length)];
if (length > MAX_CHUNK_LENGTH) secondBuffer = new byte[length - MAX_CHUNK_LENGTH];
}
public LongArray(string file)
{
using(FileStream stream = File.OpenRead(file))
{
long length = stream.Length;
firstBuffer = new byte[Math.Min(MAX_CHUNK_LENGTH, length)];
if (length > MAX_CHUNK_LENGTH) secondBuffer = new byte[length - MAX_CHUNK_LENGTH];
stream.Read(firstBuffer, 0, firstBuffer.Length);
if(secondBuffer != null) stream.Read(secondBuffer, 0, secondBuffer.Length);
}
}
public long Length => firstBuffer.LongLength + (secondBuffer == null ? 0 : secondBuffer.LongLength);
public byte this[long index]
{
get => index < MAX_CHUNK_LENGTH ? firstBuffer[index] : secondBuffer[index - MAX_CHUNK_LENGTH];
set
{
if (index < MAX_CHUNK_LENGTH) firstBuffer[index] = value;
else secondBuffer[index - MAX_CHUNK_LENGTH] = value;
}
}
public void Save(string file)
{
using(FileStream stream = File.OpenWrite(file))
{
stream.Write(firstBuffer, 0, firstBuffer.Length);
if (secondBuffer != null) stream.Write(secondBuffer, 0, secondBuffer.Length);
}
}
}
private static int[] BIT_TO_INDEX = new int[] { 1, 7, 11, 13, 17, 19, 23, 29 };
private static int[] INDEX_TO_BIT = new int[] {
-1, 0,
-1, -1, -1, -1, -1, 1,
-1, -1, -1, 2,
-1, 3,
-1, -1, -1, 4,
-1, 5,
-1, -1, -1, 6,
-1, -1, -1, -1, -1, 7,
};
private LongArray Data;
public Wheel235(long length)
{
length = (length + 29) / 30 * 30;
Data = new LongArray(length / 30);
for (long i = 0; i < Data.Length; i++) Data[i] = byte.MaxValue;
for (long i = 7; i * i < Length; i++)
{
if (!IsPrime(i)) continue;
for (long d = i * i; d < Length; d += i) ClearPrime(d);
}
}
public Wheel235(string file)
{
Data = new LongArray(file);
}
public void Save(string file) => Data.Save(file);
public long Length => Data.Length * 30L;
public bool IsPrime(long n)
{
if (n >= Length) throw new ArgumentException("Number too big");
if (n <= 5) return n == 2 || n == 3 || n == 5;
int bit = INDEX_TO_BIT[n % 30];
if (bit < 0) return false;
return (Data[n / 30] & (1 << bit)) != 0;
}
private void ClearPrime(long n)
{
int bit = INDEX_TO_BIT[n % 30];
if (bit < 0) return;
Data[n / 30] &= (byte)~(1 << bit);
}
}
Como resultado, obtuvimos una clase que puede calcular números primos que no superan los cien mil millones y guardar el resultado en el disco (tuve este código en mi computadora portátil durante 50 minutos; el tamaño del archivo creado fue de aproximadamente 3 gigabytes, como se esperaba): long N = 100_000_000_000L;
Wheel235 wheel = new Wheel235(N);
wheel.Save("./primes.dat");
Y luego lea desde el disco y use: DateTime start = DateTime.UtcNow;
Wheel235 loaded = new Wheel235("./primes.dat");
DateTime end = DateTime.UtcNow;
Console.WriteLine($"Database of prime numbers up to {loaded.Length:N0} loaded from file to memory in {end - start}");
long number = 98_000_000_093L;
Console.WriteLine($"{number:N0} is {(loaded.IsPrime(number) ? "prime" : "not prime")}!");