A non-static Sieve of Eratosthenes class, version 1

Question

Over the years I’ve seen many C# sieves and despite their varying internals all that I have seen invariably are a static method that returns a List<int>, perhaps with the occasional rare exception of an Enumerable<int>. This is fine for small collections. Granted relatively speaking for some people a 1 million item list is small. But if dealing with very large upper limits, say int.MaxValue, you could run into an out-of-memory exception fetching over 105 million items.

So I wrote a sieve class that has these highlights:

• Internally uses BitArray to track odd numbers only.
• Class implements IEnumerable with custom enumerator to walk over BitArray to find true bits.
• Lazy evaluation.
• Methods such as IsPrime(number) and FindNthPrime(sequence). FindNthPrime(-1) returns last prime found, equivalent to FindNthPrime(Count).
• On-demand count is performed in parallel. Counting also caches subtotals in ranges of 250K bits to boost FindNthPrime.
• ToNumber and ToIndex functions add context and clarity for cleaner code.
• And for those who still only want to get a List<int> returned, I included a convenient method GetPrimes.

Usings:

using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Threading;
using System.Threading.Tasks;
using System.Collections.Concurrent;

SieveOfEratosthenes:

    public class SieveOfEratosthenes : IEnumerable
    {
        public static SieveOfEratosthenes CreateInstance(int upperLimit)
        {
            if (upperLimit < 2)
            {
                throw new ArgumentException("Upper Limit be must greater than or equal to 2.");
            }
            return new SieveOfEratosthenes(upperLimit);
        }
    private SieveOfEratosthenes(int upperLimit)
    {
        _upperLimit = upperLimit;
    }

    // In case you were wondering ... the maximum length of _bits will be 1,073,741,823 when _upperLimit = int.MaxValue.
    // That's over 1 billion bits which comes out to be 134,217,728 bytes.
    private BitArray _bits = null;
    private int _knownCount = 0;
    private int _upperLimit = 2;
    private const int _offset = 3;

    // Also in case you were wondering ... the maximum number of entries in _countsCache SortedList:
    //      rangeSize 500K produces 2148 entries
    //      rangeSize 250K produces 4096 entries
    // The SortedList is only populated when needed.
    private SortedList<int, int> _countsCache = new SortedList<int, int>();
    private const int _rangeSize = 250000;

    // Tried methods, delegates, and inline code and all performed equally.
    public  static Func<int, int> ToNumber = delegate(int index) { return (2 * index) + _offset; };
    private static Func<int, int> ToIndex = delegate(int number) { return (number - _offset) / 2; };

    private void CheckBits()
    {
        if (!HasSieveBeenRun)
        {
            RunSieve();
        }
    }

    public void CheckCount()
    {
        if (_knownCount > 0) return;
        _countsCache = new SortedList<int, int>();
        _knownCount = 1;
        if (_upperLimit == 2)
        {
            return;
        }

        CheckBits();

        // Initialize sorted list in sorted order.
        // I create the list before the parallel loop so the parallel writes won't be in conflict later.
        for (int i = 0; i < _bits.Length; i += _rangeSize)
        {
            _countsCache.Add(i, 0);
        }

        var ranges = Partitioner.Create(0, _bits.Length, _rangeSize);
        Parallel.ForEach(ranges, range =>
        {
            _countsCache[range.Item1] = CountPrimesRange(range.Item1, range.Item2);
        });

        foreach (var entry in _countsCache)
        {
            _knownCount += entry.Value;
        }
    }

    public int Count
    {
        get
        {
            // Lazy evaluation. 
            // If _upperLimit = int.MaxValue, we also get SLUGGISH evaluation.  
            CheckCount();
            return _knownCount;
        }
    }

    private int CountPrimesRange(int inclusiveStart, int exclusiveEnd)
    {
        var subtotal = 0;
        for (int i = inclusiveStart; i < exclusiveEnd; i++)
        {
            if (_bits[i])
            {
                subtotal++;
            }
        }
        return subtotal;
    }

    public int Estimate
    {
        get { return (int)(_upperLimit / Math.Log(_upperLimit)); }
    }

    public int FindNthPrime(int sequence)
    {
        // Finding the nth prime may take time to Count first, if it hasn't already been done.
        if (sequence == 1)
        {
            return 2;
        }
        if (sequence == -1)
        {
            sequence = Count;
        }
        if ((sequence < 1) || (sequence > Count))
        {
            throw new ArgumentException(string.Format("Sequence must be between 1 and Known Count of {0} inclusively.", Count.ToString("N0")));
        }

        var tally = FindQuickTally(sequence);
        var counter = tally.Item2; 
        var prime = 0;
        for (int i = tally.Item1; i < _bits.Length; i++)
        {
            if (_bits[i])
            {
                if (++counter == sequence)
                {
                    prime = ToNumber(i);
                    break;
                }
            }
        }

        return prime;
    }

    private Tuple<int, int> FindQuickTally(int sequence)
    {
        var index = 0;
        var totalRangeStart = 1; // remember to include 2 (our first prime)
        foreach (var entry in _countsCache)
        {
            index = entry.Key;
            var totalRangeEnd = totalRangeStart + entry.Value;
            if (sequence <= totalRangeEnd)
            {
                break;
            }
            totalRangeStart = totalRangeEnd;
        }
        return new Tuple<int, int>(index, totalRangeStart);
    }

    IEnumerator IEnumerable.GetEnumerator()
    {
        return (IEnumerator)GetEnumerator();
    }

    public EratosthenesEnumerator GetEnumerator()
    {
        CheckBits();
        return new EratosthenesEnumerator(_bits);
    }

    public static IList<int> GetPrimes(int upperLimit)
    {
        // This is included for people who like getting just a list of primes back.
        // Worst Case: _upperLimit = int.MaxValue
        //    _bits is 134,217,728 bytes or 128 MB.
        //    The output list is 105,097,565 items needing 420,390,260 bytes or 401 MB.
        //    All in all, you need 529 MB RAM for the worst case.
        var sieve = CreateInstance(upperLimit);
        sieve.RunSieve();
        var list = new List<int>(capacity: sieve.Estimate);
        try
        {
            foreach (var prime in sieve)
            {
                list.Add(prime);
            }
        }
        catch (OutOfMemoryException)
        {
            throw new Exception("Too many primes caused an out of memory exception.  Consider iterating over sieve like 'foreach (var prime in sieve)' instead.");
        }
        return list;
    }

    public bool HasSieveBeenRun
    {
        get { return _bits != null; }
    }

    public bool IsComposite(int number)
    {
        return !IsPrime(number);
    }

    public bool IsPrime(int number)
    {
        if ((number < 2) || (number > _upperLimit))
        {
            throw new ArgumentException(string.Format("Number must be between 2 and Upper Limit of {0} inclusively.", _upperLimit.ToString("N0")));
        }
        if (number % 2 == 0) return number == 2;
        CheckBits();
        return _bits[ToIndex(number)];
    }

    public void RunSieve()
    {
        if (HasSieveBeenRun) return;
        _knownCount = 0;
        if (_upperLimit == 2)
        {
            _bits = new BitArray(0);
            return;
        }
        _bits = new BitArray(ToIndex(_upperLimit) + 1, defaultValue: true);
        var upperSqrtIndex = ToIndex((int)Math.Sqrt(_upperLimit));
        for (var i = 0; i <= upperSqrtIndex; i++)
        {
            // If this bit has already been turned off, then its associated number is composite. 
            if (!_bits[i]) continue;
            var number = ToNumber(i);
            // Okay, so number is now known to be prime.  
            // However, any multiples of number are composite and their respective bits should be turned off.
            for (var j = ToIndex(number * number); (j > i) && (j < _bits.Length); j += number)
            {
                _bits[j] = false;
            }
        }
    }

    public int UpperLimit { get { return _upperLimit; } }
}

EratosthenesEnumerator

public class EratosthenesEnumerator : IEnumerator
{
    private BitArray _bits = null;

    // Enumerators are positioned before the first element until the first MoveNext() call.
    private int _primePosition = -1;
    private int _bitPosition = -1;

    public EratosthenesEnumerator(BitArray bits)
    {
        _bits = bits;
    }

    public bool MoveNext()
    {
        _primePosition++;
        if (_primePosition > 0)
        {
            var found = -1;
            for (var i = _bitPosition + 1; i < _bits.Length; i++)
            {
                if (_bits[i])
                {
                    found = i;
                    break;
                }
            }
            _bitPosition = (found >= 0) ? found : _bits.Length;
        }
        return (_primePosition >= 0) && (_bitPosition < _bits.Length);
    }

    public void Reset()
    {
        _primePosition = -1;
        _bitPosition = -1;
    }

    object IEnumerator.Current
    {
        get
        {
            return Current;
        }
    }

    public int Current
    {
        get
        {
            try
            {
                if (_primePosition == 0)
                {
                    return 2;
                }
                return SieveOfEratosthenes.ToNumber(_bitPosition);
            }
            catch (IndexOutOfRangeException)
            {
                throw new InvalidOperationException();
            }
        }
    }
}

Not Interested In

This class uses a BitArray and int. PERIOD.

Don’t bother suggesting that I could perform parallel writes with a ConcurrentList<bool>. I am not interested in that here. Yet I would be happy to read a new, original thread created by you where you use ConcurrentList.

Don’t bother suggesting that I could use larger integral types. I am not interested in that here. Yet I would be happy to read a new, original thread created by you where you use uint, long, ulong, or even BigInteger.

The cached counts is a feature that I refuse to remove but am open to improvements. It’s called on-demand, trades-off a small but acceptable amount of memory, and really boosts performance of FindNthPrime.

Concerns

This was my first time to use a function delegate, as well as a custom enumerator. The ToNumber function is declared public static in the SieveOfEratosthenes class, and for the sake of DRY is referenced within the EratosthenesEnum class. Since this was my first time to use either, I don’t know if it’s a proper use or not.

ToIndex is the yin to ToNumber ‘s yang but since the enumerator doesn’t use it, I have its access as private. Should it stay private or should it be public like it twin ToNumber?

Usage

I tried to make the usage feel natural. Consider:

var sieve = SieveOfEratosthenes.CreateInstance(limit);
sieve.RunSieve();
foreach (var prime in sieve)
{
    /* do something */
}
sieve.CheckCount();
Console.WriteLine("Count of primes {0}", sieve.Count);

But since the class implements IEnumerable which calls RunSieve if needed, and Count calls CheckCheck, this performs the same thing:

var sieve = SieveOfEratosthenes.CreateInstance(limit);
foreach (var prime in sieve)
{
    /* do something */
}
Console.WriteLine("Count of primes {0}", sieve.Count);

The only reason I have RunSieve and CheckCount as public methods is for timing them individually with a Stopwatch.

The sad thing is this is a nice class but the things that work best and fastest with it - for a reasonable upper limit - are the things that aren’t needed since you can simply output as a List<int>. And the things that were put in to address really large upper limits are the things that will still be sluggish.

Nonetheless with this class, I can at least iterate over all of the int primes without an out-of-memory exception, which is something I can’t do with the other sieves that I’ve seen.

Answer 1

I did find a spot that could cause problems. The GetPrimes should be changed to keep the declaration for variable list before the try block but move the initialization using an estimated capacity inside the try where the out-of-memory exception is caught.

The out-of-memory exception I get with upper limit of int.MaxValue is REAL. And note the class does NOT use bool[] or int[]. It uses the memory efficient BitArray for odd numbers only – and I still get the exception.

Yes it’s a lot of code for a simple sieve, which I alluded to in the OP (see 2nd to last paragraph). But this sieve has beneficial features that make it more than simple. Let me clearly state: if you can output the primes to List<int>, then you should do so. I whole heartily recommend trying that first. But if you can’t because you get an out-of-memory exception – an ACTUAL one not a theoretical one – then my sieve has some benefits.

Upper Limit int.MaxValue

There will be over 105 million primes found. Below I present 2 problems to solve when using upper limit int.MaxValue. I already know I can’t dump it to a list without getting an exception. So I am forced to iterate over the sieve.

Consider this simple serial counter:

// Approximtely 8.0 seconds
var counter = 0;
foreach (var prime in sieve)
{
    counter++;
}

That takes 8 seconds by itself. Anything else I do will only add to that time.

Problem 1: Randomly List 100 Primes

Offhand I would create a dictionary with 100 entries. The key would hold the index number, and the value would hold the prime that is found at the index. Inside the loop you would have to see if the dictionary contains the current prime. I already know this will take a minimum of 8 seconds.

Unless I do this:

// Approximtely 0.11 seconds
var random = new Random();
var primes = new List<int>();
for (int i = 0; i < 100; i++)
{
    var j = random.Next(sieve.Count);
    primes.Add(sieve.FindNthPrime(j + 1));
}

I don’t need the dictionary. I can just add the desired primes to a list. In less than 0.11 seconds.

Problem 2: List Primes every power of 10

I want to find primes at 1, 10, 100, 1000, etc., up to 100 million. With a simple sieve, this will again take you at least 8 seconds. Yet this takes a blistering 0.015 seconds:

// Approximtely 0.015 seconds
var dict = new Dictionary<int, int>();
var maxPower = Math.Log10(sieve.Count);
for (int i = 0; i <= maxPower; i++)
{
    var j = (int)Math.Pow(10, i);
    dict[j] = sieve.FindNthPrime(j);
}

And sure I get add more memory to my laptop, but if I did that I would probably rewrite this for uint and hit its memory limits!