Base de données des nombres premiers jusqu'à cent milliards sur le genou

image

L'algorithme le plus connu pour trouver tous les nombres premiers non supérieurs à un nombre donné est le tamis d'Ératosthène . Cela fonctionne très bien pour des nombres allant jusqu'à des milliards, peut-être jusqu'à des dizaines de milliards, s'ils sont soigneusement écrits. Cependant, tous ceux qui aiment s'amuser avec les nombres premiers savent que vous voulez toujours en avoir autant que possible. Une fois, pour résoudre un problème sur un hackerrank, j'avais besoin d'une base de données en mémoire de nombres premiers jusqu'à cent milliards. Avec une optimisation maximale par la mémoire, si le tamis d'Eratosthène présente des nombres impairs avec un tableau de bits, sa taille sera d'environ 6 gigaoctets, ce qui ne correspondait pas à la mémoire de mon ordinateur portable. Il y a une modification de l'algorithme qui est beaucoup moins exigeante en mémoire (division de la plage de nombres d'origine en plusieurs morceaux et traitement d'une pièce à la fois) -tamis segmenté d'Eratosthène , mais il est plus difficile à mettre en œuvre, et le résultat global ne rentrera de toute façon pas dans la mémoire. Ci-dessous, j'attire votre attention sur un algorithme qui est presque aussi simple que le tamis Eratosthène, mais qui donne une double optimisation par la mémoire (c'est-à-dire qu'une base de données de nombres premiers jusqu'à cent milliards occupera environ 3 gigaoctets, qui devraient déjà entrer dans la mémoire d'un ordinateur portable standard).

Pour commencer, nous rappelons comment fonctionne le tamis d'Ératosthène: dans un tableau de nombres de 1 à N, les nombres divisibles par deux sont barrés, puis les nombres divisibles par 3, et ainsi de suite. Les chiffres restants après la traversée seront 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  ...

Code 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];
}


Sur mon ordinateur portable, il compte jusqu'à 15 milliards de secondes et consomme un gigaoctet de mémoire.

L'algorithme que nous utiliserons s'appelle Wheel Factorization . Pour comprendre son essence, nous écrivons des nombres naturels avec une tablette de 30 pièces consécutives:

 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
...

Et biffez tous les nombres divisibles par 2, 3 et 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
...

On peut voir que dans chaque ligne, les numéros sont barrés aux mêmes positions. Comme nous n'avons barré que des nombres composites (à l'exception, en fait, des nombres 2, 3, 5), les nombres premiers ne peuvent être que dans des positions non marquées.

Puisqu'il y a huit de ces nombres dans chaque ligne, cela nous donne l'idée de représenter chaque trente nombres avec un octet, dont chacun des huit bits indiquera la simplicité du nombre dans la position correspondante. Cette coïncidence est non seulement esthétiquement belle, mais simplifie également considérablement la mise en œuvre de l'algorithme!

image

Implémentation technique


Premièrement, notre objectif est d'atteindre cent milliards, de sorte que les nombres à quatre octets ne peuvent plus être supprimés. Par conséquent, notre algorithme ressemblera à ceci:

class Wheel235
{
    private byte[] Data;

    public Wheel235(long length)
    {
        ...
    }
    public bool IsPrime(long n)
    {
        ...
    }
    ...
}

Pour faciliter la mise en œuvre, nous arrondirons la longueur au nombre le plus proche divisible par 30. Ensuite,

class Wheel235
{
    private byte[] Data;

    public Wheel235(long length)
    {
        // ensure length divides by 30
        length = (length + 29) / 30 * 30; 
        Data = new byte[length/30];
        ...
    }

    public long Length => Data.Length * 30L;

    ...
}

Ensuite, nous avons besoin de deux tableaux: l'un traduit l'index 0..29 en nombre de bits, le second, au contraire, le nombre de bits en position trente. En utilisant ces tableaux et un peu d'arithmétique, nous construisons une correspondance biunivoque entre les nombres naturels et les bits dans le tableau de données. Pour simplifier, nous utilisons seulement deux méthodes: IsPrime (n), qui vous permet de savoir si le nombre est premier, et ClearPrime (n), qui définit le bit correspondant à 0, marquant le nombre n comme composé:

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);
    }
}

Notez que si le nombre est divisible par 2, 3 ou 5, alors le nombre de bits sera -1, ce qui signifie que le nombre est évidemment composite. Eh bien et, bien sûr, nous traitons un cas spécial n <= 5.

Maintenant, en fait, l'algorithme lui-même, qui répète presque littéralement le tamis d'Eratosthène:

    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 conséquence, l'ensemble de la classe résultante:

Ressemble à ça
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 Wheel235(long length)
    {
        // ensure length divides by 30
        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);
    }
}


Il n'y a qu'un petit problème: ce code ne fonctionne pas pour des nombres supérieurs à environ 65 milliards! Lorsque vous essayez de le démarrer avec de tels nombres, le programme se bloque avec une erreur:

Unhandled Exception: System.OverflowException: Array dimensions exceeded supported range.

Le problème est que dans les tableaux c # ne peut pas avoir plus de 2 ^ 31 éléments (apparemment parce que l'implémentation interne utilise un index de tableau à quatre octets). Une option consiste à créer, par exemple, un tableau long [] au lieu d'un tableau d'octets [], mais cela compliquera un peu l'arithmétique des bits. Pour plus de simplicité, nous irons dans l'autre sens, en créant une classe spéciale qui simule le tableau souhaité, contenant deux tableaux courts à l'intérieur. Dans le même temps, nous lui donnerons la possibilité de nous enregistrer sur le disque afin que vous puissiez calculer une fois les nombres premiers, puis réutiliser la base de données:

    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)
        {
            ...
        }

    }


Enfin, combinez toutes les étapes de

version finale de l'algorithme
class 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)
    {
        // ensure length divides by 30
        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);
    }

}


En conséquence, nous avons obtenu une classe qui peut calculer des nombres premiers ne dépassant pas cent milliards et enregistrer le résultat sur le disque (j'ai eu ce code sur mon ordinateur portable pendant 50 minutes; la taille du fichier créé était d'environ 3 gigaoctets, comme prévu):

        long N = 100_000_000_000L;
        Wheel235 wheel = new Wheel235(N);
        wheel.Save("./primes.dat");

Et puis lisez à partir du disque et utilisez:

        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")}!");

All Articles