4
\$\begingroup\$

I have written many variations of the Sieve of Eratosthenses, which is the fastest way to generate a large collection of primes. (@slepic asked for clarification, which I provide in an answer below. My intended statement is that a sieve in general is much faster than naive methods at generating lots of primes; not that the Sieve of Eratosthenses is the fastest ever.)

If you later want to query the gathered primes by count or at a specific index, the performance of the sieve is lacking compared to a list. So I thought, why not make a prime table that uses a high performance sieve to generate the primes, but later move those primes into a list (memory permitting).

I originally wrote this as an answer to someone else's post, but much of my goals, objectives, code, and features differed so much that I am posting for my own review.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Diagnostics;
using System.Collections;

namespace Prime_Table_Core
{
    // What's in a name?  Variable/parameter names for any Int32 were chosen to denote context. 
    //
    //    number: any Int32 on the "number line" to be evaluated as prime, composite, or neither.
    //    prime : a subset of number where the Int32 is a prime.
    //    index : an Int32 used as the positional index into _knownPrimes list.
    //    value : no specific context or restriction on this Int32.

    public static class PrimeTable
    {
        private static readonly List<int> _knownPrimes = new List<int>() { 2 };

        public static bool IsInitialized { get; private set; } = false;
        public static TimeSpan LastDuration { get; private set; } = TimeSpan.Zero;

        // If you want to work directly with just the known primes, no need for streaming
        // since the table is already in memory.
        public static IReadOnlyList<int> KnownPrimes => _knownPrimes;

        public static int KnownPrimeCount => _knownPrimes.Count;
        public static int LastKnownPrime => _knownPrimes.Last();
        public static int LastKnownIndex => _knownPrimes.Count - 1;

        // Track the very last number checked using GetNextUnknownPrime() or Initialize().
        // This number could be greater than LastKnownPrime.
        private static int _lastNumberChecked = 2;

        private static Func<int, bool> HasMoreNumbers = number => (int.MaxValue - number) > 2;
        private static Func<int, int> DoubleIt = value => value << 1;
        private static Func<int, int> HalveIt = value => value >> 1;
        private static Func<int, bool> IsEven = value => value % 2 == 0;

        public static int GetIndexAtOrBefore(int number)
        {
            if (number < 2)
            {
                return -1;
            }

            InitializeIfNeeded();

            if (number >= LastKnownPrime)
            {
                return LastKnownIndex;
            }

            var upperIndex = LastKnownIndex;
            var lowerIndex = 0;
            var midIndex = HalveIt(upperIndex + lowerIndex);

            // Instead of a while(true), let's completely avoid an infinite loop.
            // The for loop won't use it's index variable other than to prevent
            // the loop from being infinite.  But as a debugging bonus, you can use
            // "iteration" to see how many iterations were needed for a lookup.
            for (var iteration = 1; iteration < _knownPrimes.Count; iteration++)
            {
                if (number == _knownPrimes[midIndex])
                {
                    return midIndex;
                }

                if ((upperIndex - lowerIndex) <= 1)
                {
                    return (number > _knownPrimes[upperIndex]) ? upperIndex : lowerIndex;
                }

                if (number > _knownPrimes[midIndex])
                {
                    lowerIndex = midIndex;
                }
                else
                {
                    upperIndex = midIndex;
                }

                midIndex = HalveIt(upperIndex + lowerIndex);
            }

            return -1;  // for safety's sake, but really is unreachable.
        }

        public static int GetIndexBefore(int number) => (number <= 2) ? -1 : GetIndexAtOrBefore(number - 1);
        public static int GetIndexAfter(int number) => (number == int.MaxValue) ? -1 : GetIndexAtOrAfter(number + 1);
        public static int GetIndexAtOrAfter(int number)
        {
            var index = GetIndexAtOrBefore(number);
            if (index == -1)
            {
                return 0;
            }
            if (_knownPrimes[index] == number)
            {
                return index;
            }
            return ++index < KnownPrimeCount ? index : -1;
        }

        public static bool IsPrime(this int number)
        {
            // First, dispense with easy cases.
            if (number < 2) { return false; }
            if (IsEven(number)) { return number == 2; }

            InitializeIfNeeded();

            var index = 0;

            // Second, quickly check against _knownPrimes and _lastNumberChecked.
            if (number <= LastKnownPrime)
            {
                index = GetIndexAtOrBefore(number);
                return _knownPrimes[index] == number;
            }
            if (number <= _lastNumberChecked)
            {
                return false;
            }

            // Third, perform naive primality test using known primes.
            var sqrt = (int)Math.Sqrt(number);

            for (index = 0; index < _knownPrimes.Count; index++)
            {
                if (number % _knownPrimes[index] == 0)
                {
                    return false;
                }
                if (_knownPrimes[index] > sqrt)
                {
                    return true;
                }
            }

            //  Fourth, perform naive primality test on Odds beyond LargestKnownPrime
            for (var possibleDivisor = _lastNumberChecked + 2; possibleDivisor <= sqrt; possibleDivisor += 2)
            {
                if (number % possibleDivisor == 0)
                {
                    return false;
                }
            }

            // Finally, it must be prime.
            return true;
        }

        // This method will stream the known primes first, followed by the unknown ones.
        public static IEnumerable<int> GetPrimes()
        {
            InitializeIfNeeded();

            foreach (var prime in _knownPrimes)
            {
                yield return prime;
            }

            for (; ; )
            {
                var next = GetNextUnknownPrime();
                if (next.HasValue)
                {
                    yield return next.Value;
                }
                else
                {
                    yield break;
                }
            }
        }

        // This method bypasses the known primes and starts streaming the unknown ones, if any.
        public static IEnumerable<int> GetUnknownPrimes()
        {
            InitializeIfNeeded();

            for (; ; )
            {
                var next = GetNextUnknownPrime();
                if (next.HasValue)
                {
                    yield return next.Value;
                }
                else
                {
                    yield break;
                }
            }
        }

        public static int? GetNextUnknownPrime()
        {
            if (!HasMoreNumbers(_lastNumberChecked))
            {
                LastDuration = TimeSpan.Zero;
                return null;
            }

            int result = -1;

            InitializeIfNeeded();

            var sw = Stopwatch.StartNew();

            for (var candidate = _lastNumberChecked + 2; ; candidate += 2)
            {
                if (IsPrime(candidate))
                {
                    _lastNumberChecked = candidate;
                    result = candidate;
                    break;
                }
                _lastNumberChecked = candidate;
                if (!HasMoreNumbers(candidate))
                {
                    // Do this here instead of inside for condition so that
                    // we do not overflow past Int.MaxValue, or worse,
                    // wrap around to Int.MinValue.
                    break;
                }
            }

            if (result > 1)
            {
                _knownPrimes.Add(result);
            }

            sw.Stop();
            LastDuration = sw.Elapsed;
            return result;
        }

        // This will only initialize _knownPrimes once.
        public static void InitializeIfNeeded()
        {
            const int DefaultUpperLimit = 1_500_001;    // produces   114_155 primes in 0.01 seconds
            if (!IsInitialized)
            {
                Initialize(DefaultUpperLimit);
            }
        }

        // You may Initialize and re-Initialize to your heart's content.
        // Depending upon upperLimit, this may take a split second or half a minute or longer based
        // upon your CPU and RAM.
        public static void Initialize(int upperLimit)
        {
            const int MinimumUpperLimit = 1000;

            if (upperLimit < MinimumUpperLimit)
            {
                throw new ArgumentException($"{nameof(upperLimit)} must be {MinimumUpperLimit} or greater.");
            }

            var sw = Stopwatch.StartNew();

            GenerateSieve(upperLimit);

            sw.Stop();
            LastDuration = sw.Elapsed;
            IsInitialized = true;
        }

        // The intent is to start off with a small, very fast sieve to build the _knownPrimes up to a point.
        // While a BitArray uses less memory, it is also slower than bool[].
        // Once this method completes, the array is set to null and memory can be GC'd.
        // If responsiveness is your goal, then a "reasonable" upperLimit is one that executes 
        // in less than 0.25 seconds on your hardware.
        private static void GenerateSieve(int upperLimit)
        {
            lock (_knownPrimes)
            {
                _knownPrimes.Clear();
                _knownPrimes.Add(2);

                // Evens all done.  Now check only odd numbers for primality

                if (IsEven(upperLimit))
                {
                    upperLimit++;
                }

                const int offset = 1;
                Func<int, int> ToNumber = index => DoubleIt(index) + offset;
                Func<int, int> ToIndex = number => HalveIt(number - offset);

                // initial flags are false
                var flags = new BitArray(ToIndex(upperLimit) + 1, true);
                flags[0] = false;

                var upperSqrtIndex = ToIndex((int)Math.Sqrt(upperLimit));

                for (var i = 1; i <= upperSqrtIndex; i++)
                {
                    // If this bit has already been turned off, then its associated number is composite. 
                    if (!flags[i]) { continue; }
                    var number = ToNumber(i);
                    _knownPrimes.Add(number);
                    // Any multiples of number are composite and their respective flags should be turned off.
                    for (var j = ToIndex(number * number); j < flags.Length; j += number)
                    {
                        flags[j] = false;
                    }
                }

                // Output remaining primes once flags array is fully resolved:
                for (var i = upperSqrtIndex + 1; i < flags.Length; i++)
                {
                    if (flags[i])
                    {
                        _knownPrimes.Add(ToNumber(i));
                    }
                }

                _lastNumberChecked = upperLimit;
            }
        }
    }
}

This was written in .NET Core 3.0, but also ported to full Framework 4.8. The full Framework is about 50% slower on the same hardware.

Once the prime table is generated, you may query against the list of what I call known primes. But you may also continue to discover unknown primes, if any, that once discovered are then added to the known primes.

You may quickly initialize a larger number of known primes using the Initialize(upperLimit) method. If speedy responsiveness is your main objective, then a good upperlimit should be something that returns in 0.25 seconds or less on your particular hardware. If you want to max out all of Int32, you can do that too but it may take quite a while to generate all 105+ million primes.

An example of it in use:

PrimeTable.Initialize using assorted upper limits:
   Upper Limit = 1000001, PrimeCount = 78498, LastPrime = 999983, Duration: 00:00:00.0064373  (includes JIT time)
   Upper Limit = 1500001, PrimeCount = 114155, LastPrime = 1499977, Duration: 00:00:00.0043673
   Upper Limit = 2000001, PrimeCount = 148933, LastPrime = 1999993, Duration: 00:00:00.0072214
   Upper Limit = 5000001, PrimeCount = 348513, LastPrime = 4999999, Duration: 00:00:00.0180426
   Upper Limit = 10000001, PrimeCount = 664579, LastPrime = 9999991, Duration: 00:00:00.0330480
   Upper Limit = 17000001, PrimeCount = 1091314, LastPrime = 16999999, Duration: 00:00:00.0573246
   Upper Limit = 20000001, PrimeCount = 1270607, LastPrime = 19999999, Duration: 00:00:00.0648279
   Upper Limit = 50000001, PrimeCount = 3001134, LastPrime = 49999991, Duration: 00:00:00.1564291

Demo of index usage to KnownPrimes:
   GetIndexAtOrBefore(55551) = 5636, KnownPrimes[5636] = 55547
   GetIndexAtOrAfter (55551) = 5637, KnownPrimes[5637] = 55579

Demo fetching next 10 unknown primes:
   PrimeCount = 3001135, LastPrime = 50000017, Duration: 00:00:00.0004588  (includes JIT time)
   PrimeCount = 3001136, LastPrime = 50000021, Duration: 00:00:00.0000044
   PrimeCount = 3001137, LastPrime = 50000047, Duration: 00:00:00.0000188
   PrimeCount = 3001138, LastPrime = 50000059, Duration: 00:00:00.0000065
   PrimeCount = 3001139, LastPrime = 50000063, Duration: 00:00:00.0000180
   PrimeCount = 3001140, LastPrime = 50000101, Duration: 00:00:00.0000048
   PrimeCount = 3001141, LastPrime = 50000131, Duration: 00:00:00.0000071
   PrimeCount = 3001142, LastPrime = 50000141, Duration: 00:00:00.0000193
   PrimeCount = 3001143, LastPrime = 50000161, Duration: 00:00:00.0000097
   PrimeCount = 3001144, LastPrime = 50000201, Duration: 00:00:00.0000148

PrimeTable.Initialize(int.MaxValue):
   Upper Limit = 2147483647, PrimeCount = 105097565, LastPrime = 2147483647, Duration: 00:00:12.8353907
   GetIndexAtOrBefore(55551) = 5636, KnownPrimes[5636] = 55547
   GetIndexAtOrAfter (55551) = 5637, KnownPrimes[5637] = 55579
   GetIndexAtOrAfter (2147483647) = 105097564, KnownPrimes[105097564] = 2147483647
   GetIndexAfter (2147483647) = -1
   GetNextUnknownPrime() = <null>

Press ENTER key to close

There are 3 ways to enumerate over a large collection of primes:

  1. Use the KnownPrimes table, a read-only list.
  2. GetUnknownPrimes() skips over the known primes and streams the unknown to you.
  3. GetPrimes() will first stream the known primes to you, followed by the unknown.

Other features:

Since performance is a curiosity, there is a LastDuration property to inform you how long the sieve took to generate, or how long the last GetNextUnknownPrime took.

Anything using the known prime's index does not discover any unknown primes. This includes the IsPrime method, which is a tad long as it tries to first check against the known primes before resorting to a naive implementation.

\$\endgroup\$
1
  • 2
    \$\begingroup\$ "which is the fastest way to generate a large collection of primes". That deserves some reference... \$\endgroup\$
    – slepic
    Dec 1, 2019 at 8:10

1 Answer 1

1
\$\begingroup\$

I am providing an answer to my post in order to address a comment by @slepic regarding the first sentence in the OP. @slepic asked for clarification to this statement:

I have written many variations of the Sieve of Eratosthenses, which is the fastest way to generate a large collection of primes.

First of all, what I meant was that in order to generate a lot of primes that a sieve is faster than using naive methods. There may be sieves faster than Eratosthenses, but a sieve will be much faster than not using a sieve. That was my intended statement and hopefully addresses the clarification that was requested.

My PrimeTable can be easily modified to demonstrate this. First, I changed this line in PrimeTable.cs:

public static bool IsInitialized { get; private set; } = true;

But hit a quirk because the only prime I have to start with is 2, and my later logic assumes the last known prime is odd. I could change that logic but I chose instead to change this line:

private static readonly List<int> _knownPrimes = new List<int>() { 2, 3 };

Which also required me to change a field, which was upgraded with softer coding:

private static int _lastNumberChecked = LastKnownPrime;

With those few changes, I then wrote a method to generate primes:

private static void SlowerGrowth()
{
    Console.WriteLine("Display 'slower' growth without sieve.");
    // Account for JIT
    var prime = PrimeTable.GetNextUnknownPrime();  
    var preCount = PrimeTable.KnownPrimeCount;  

    var step = TimeSpan.FromMilliseconds(10);
    var limit = TimeSpan.FromSeconds(1);
    var progressMark = step;

    var total = TimeSpan.Zero;
    var count = 0;

    while (total < limit)
    {
        prime = PrimeTable.GetNextUnknownPrime();
        var elapsed = PrimeTable.LastDuration;
        total += elapsed;

        if (total >= progressMark || total >= limit)
        {
            count++;
            Console.WriteLine($"   Count = {(PrimeTable.KnownPrimeCount - preCount)}, Largest = {PrimeTable.LastKnownPrime}, Elapsed = {total}"); //, Step = {step}, Mark = {progressMark}");
            if (count == 5 || total >= limit)
            {
                step = 10 * step;
                progressMark = step;
                count = 0;
            }
            else
            {
                progressMark += step;
            }
        }
    }
}

Which produced this output:

WITHOUT A SIEVE (NAIVE CHECKS)

Display 'slower' growth without sieve.
   Count = 16427, Largest = 181211, Elapsed = 00:00:00.0100004
   Count = 29658, Largest = 346079, Elapsed = 00:00:00.0200006
   Count = 41234, Largest = 496007, Elapsed = 00:00:00.0300001
   Count = 52233, Largest = 642197, Elapsed = 00:00:00.0400015
   Count = 62740, Largest = 783707, Elapsed = 00:00:00.0500005
   Count = 104720, Largest = 1366609, Elapsed = 00:00:00.1000005
   Count = 178155, Largest = 2427463, Elapsed = 00:00:00.2000005
   Count = 243973, Largest = 3406421, Elapsed = 00:00:00.3000012
   Count = 306982, Largest = 4363897, Elapsed = 00:00:00.4000024
   Count = 365978, Largest = 5270231, Elapsed = 00:00:00.5000013
   Count = 619977, Largest = 9280757, Elapsed = 00:00:01.0000003

I followed up by running a few different size sieves, to get these results:

WITH A SIEVE

PrimeTable.Initialize using assorted upper limits:
   Upper Limit = 10000001, PrimeCount = 664579, LastPrime = 9999991, Duration: 00:00:00.0340529  (includes JIT time)
   Upper Limit = 20000001, PrimeCount = 1270607, LastPrime = 19999999, Duration: 00:00:00.0618941
   Upper Limit = 200000001, PrimeCount = 11078937, LastPrime = 199999991, Duration: 00:00:00.9063038

Using ballpark numbers, it took the naive methods almost 1 second to generate around 620K primes with the largest near 9.3 million. Using a sieve, it took only 0.035 seconds to find the same (plus 40K more). For 1 seconds using a sieve, I could find over 11 million primes which is over 17X more than using naive methods.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.