Sieve32Fast - A very fast, memory efficient, multi-threaded Sieve of Eratosthenes

Question

An improved version Sieve32FastV2 is available.

The classical solutions for the Sieve of Eratosthenes fall into 2 camps: one uses a bool[], which is fast but very memory bloated; the other uses a BitArray, which is more sluggish but uses far less memory. Refusing to accept that those are the only 2 possibilities, I’ve created a multi-threaded sieve that is as fast as a bool[] while using only 4% more memory than a BitArray.

I don’t consider the code to be complicated but is quite a bit longer than a simple sieve. If a simple sieve is more of your speed, stick with Sieve31 or Sieve32.

The class is no longer static; however all public methods remain static. The public methods will create a private instance, and that private instance uses 2 embedded private classes as well: Vector and VectorList, which implements List<Vector>.

Related versions

Sieve31 a simple sieve for 31 bit primes, or int. Uses BitArray. Likely candidate for ToList().

Sieve32 a simple sieve for 32 bit primes, or uint. Uses BitArray. Due to memory needs, a bool[] version is not possible.

EBrown’s answer to Sieve31. Uses a bool[]. Over 40% faster than Sieve31 but uses 600% the memory. Due to memory constraints, cannot use ToList().

The version below is as fast if not faster than EBrown’s answer (depends on number of cores) but requires only 4% more memory than Sieve31.

public class Sieve32Fast
{
    private static ArgumentException BadUpperLimitException => new ArgumentException("upperLimit be must greater than or equal to 2.", "upperLimit");

    // NOTE TO 'Sam The Maintainer'
    //
    // As the code deals with primes, composites, indices, and lengths, all of which are integers,
    // it helps to have context over what type of entity a given value represents.
    //
    // A 'Number' will be a uint representing some natural number {0, 1, 2, ..., uint.MaxValue}.
    // The constants Zero, One, Two, Three are such 'Numbers'.
    // If I were to add 1 to a 'Number', I will use something like:
    //
    //      Number + One;
    //
    // An index or length, which are associated with working with arrays or lists, will be a int.
    // If I were to add 1 to an index or subtract 1 from a length, I will use something like:
    //
    //      index + 1;
    //      length - 1;

    private const uint Zero = 0U;
    private const uint One = 1U;
    private const uint Two = 2U;
    private const uint Three = 3U;

    public static IEnumerable<uint> Primes(int upperLimit)
    {
        if (upperLimit < Two) { throw BadUpperLimitException; }
        return Primes((uint)upperLimit);
    }

    public static IEnumerable<uint> Primes(uint upperLimit)
    {
        if (upperLimit < Two) { throw BadUpperLimitException; }

        var instance = new Sieve32Fast(upperLimit);
        return instance.EnumeratePrimes();
    }

    private Sieve32Fast(uint upperLimit)
    {
        _upperLimit = upperLimit;
        _vectors = VectorList.Create(_upperLimit);
    }

    private uint _upperLimit = Zero;
    private VectorList _vectors = null;

    // Note to 'Sam the Maintainer' regarding a Performance Tweak:
    // The frequently called ToIndex and ToNumber Func's require division and multiplying by 2.
    // For billions of calls, this can be slightly expensive.
    // Instead I will bit shift by 1, so that (X / 2) becomes (X >> 1) and (X * 2) becomes (X << 1).
    private static Func<uint, uint, int> ToIndex => (uint number, uint startingNumber) => (int)((number - startingNumber) >> 1);
    private static Func<int, uint, uint> ToNumber => (int bitIndex, uint startingNumber) => (uint)(bitIndex << 1) + startingNumber;

    private IEnumerable<uint> EnumeratePrimes()
    {
        if (_upperLimit < Two) { yield break; }

        yield return Two;
        if (_upperLimit == Two) { yield break; }

        // I call _vectors[0] the rootVector not just because its the very first one, but also because
        // it was intentionally created so that the index of upperLimit's square root is contained within _vectors[0].
        var rootVector = _vectors[0];
        var rootBitIndex = GetSquareRootIndex(_upperLimit);

        // The number of times a bit in all BitArray(s) are accessed:
        //
        //      UpperLimit = int.MaxValue    =>   3,315,151,693 times
        //      UpperLimit = uint.MaxValue   =>   6,701,709,402 times

        for (var bitIndex = 0; bitIndex <= rootBitIndex; bitIndex++)
        {
            if (rootVector[bitIndex])
            {
                var prime32 = ToNumber(bitIndex, Three);
                yield return prime32;
                // If prime, all of its multiples - on all vectors - are composites and should be marked as such.
                MarkCompositesInParallel(bitIndex, prime32);
            }
        }

        // output remaining primes 
        for (var vectorIndex = 0; vectorIndex < _vectors.Count; vectorIndex++)
        {
            var vector = _vectors[vectorIndex];
            var startIndex = (vectorIndex == 0) ? rootBitIndex + 1 : 0;
            // Due to high frequency of access, its ever so slightly faster to have copies created outside the loop
            // rather than called inside the loop directly and repeatedly with vector.BitLength and vector.StartingNumber.
            var length = vector.BitLength;
            var startingNumber = vector.StartingNumber;
            for (var bitIndex = startIndex; bitIndex < length; bitIndex++)
            {
                if (vector[bitIndex]) { yield return ToNumber(bitIndex, startingNumber); }
            }
        }
    }

    private void MarkCompositesInParallel(int bitIndex, uint prime32)
    {
        Parallel.For(0, _vectors.Count, vectorIndex =>
        {
            var vector = _vectors[vectorIndex];
            var startIndex = 0;
            var stopIndex = vector.BitLength - 1;

            if (vectorIndex == 0)
            {
                // startIndex may be calculated way past the BitLength of the vector.
                // That's okay as it will quickly break out of the loop below.
                var square = prime32 * prime32;
                startIndex = ToIndex(square, Three);
            }
            else
            {
                var remainder = vector.StartingNumber % prime32;
                if (remainder != Zero)
                {
                    var targetNumber = vector.StartingNumber + prime32 - remainder;
                    // On the full number scale, every other multiple of prime32 is even and should be skipped 
                    // over for the next multiple, which is an odd number.  We only want odd multiples.
                    if (remainder % Two == Zero)
                    {
                        targetNumber += prime32;
                    }
                    startIndex = ToIndex(targetNumber, vector.StartingNumber);
                }
            }

            // This could be defined once outside the lambda but I want each instance to have their own local copy,
            // primarily due to high frequency of access within the inner loop below.
            var prime31 = (int)prime32;

            // Any multiples of prime31 are composite and their respective flags should be marked as such.
            for (var i = startIndex; i <= stopIndex; i += prime31)
            {
                vector[i] = false;
            }
        });
    }

    private static int GetSquareRootIndex(uint number)
    {
        var squareRoot = (uint)Math.Sqrt(number);
        return ToIndex(squareRoot, Three);
    }

    private class VectorList : List<Vector>
    {
        public uint UpperLimit { get; private set; }

        public static VectorList Create(uint upperLimit)
        {
            var instance = new VectorList(upperLimit);
            instance.CreateVectors();
            return instance;
        }

        private VectorList(uint upperLimit)
        {
            this.UpperLimit = upperLimit;

            // Any upperLimit > 2 should be odd for working with VectorList and Vector(s).
            if ((this.UpperLimit > Two) && (this.UpperLimit % Two == Zero))
            {
                this.UpperLimit--;
            }
        }

        private void CreateVectors()
        {
            var typicalBitLength = CalcTypicalBitLength();
            var typicalNumberRange = (uint)typicalBitLength * Two;
            var count = (UpperLimit / typicalNumberRange) + One;

            for (uint i = Zero, endingNumber = One; i <= count; i++)
            {
                if (endingNumber >= UpperLimit) { break; }
                // The first vector may have to be longer to accomodate the index of UpperLimit's square root.
                var length = (i == 0) ? GetSpecialFirstLength(typicalBitLength) : typicalBitLength;
                var startingNumber = endingNumber + Two;
                var vector = new Vector(startingNumber, length, UpperLimit);
                this.Add(vector);
                endingNumber = vector.EndingNumber;
            }
        }

        private int CalcTypicalBitLength()
        {
            var length = ToIndex(UpperLimit, Three) + 1;

            // Small enough values will result in 1 vector
            const uint smallNumberCutoff = 10000;
            if (UpperLimit < smallNumberCutoff) { return length; }

            // Divide length for later parallelization over many (but not too many) vectors.
            const int tinyFactor = 8;
            var maxVectorCount = tinyFactor * Environment.ProcessorCount;
            length = (length / maxVectorCount) + 1;

            return PaddedLength(length);
        }

        private int GetSpecialFirstLength(int length)
        {
            // For the very first vector, aka the root vector, make sure the index of the upper limit's square root is in the first vector.
            var rootIndex = GetSquareRootIndex(UpperLimit);
            if (rootIndex < length) { return length; }
            return PaddedLength(rootIndex + 1);
        }

        private static int PaddedLength(int length)
        {
            // BitArray internally uses 32 bit int[] so align upwards to a 32 bit boundary, 
            // i.e. pad the end of length (in bits) to consume a full 32 bit int.
            var remainder = length % 32;
            return (remainder == 0) ? length : length + 32 - remainder;
        }
    }

    // A Vector is aware of its bits and length, as well as its starting and ending number.
    // A Vector is unaware of other vectors or that it is a member of a collection of vectors or that it has an index into such collections.
    private class Vector
    {
        private BitArray _bits = null;

        public Vector(uint startNumber, int length, uint upperLimit)
        {
            StartingNumber = startNumber;
            var endNumber = startNumber + (Two * (length - 1));
            if (endNumber > upperLimit)
            {
                length = ToIndex(upperLimit, StartingNumber) + 1;
            }
            _bits = new BitArray(length, defaultValue: true);
        }

        public bool this[int index] { get { return _bits[index]; } set { _bits[index] = value; } }

        public int BitLength => _bits.Length; 

        public uint StartingNumber { get; private set; }

        public uint EndingNumber => ToNumber(_bits.Length - 1, StartingNumber);
    }
}

About BitArray

Curious readers may want to review Microsoft’s Reference Source on the BitArray Class. At the very top is this:

A vector of bits. Use this to store bits efficiently, without having to do bit shifting yourself.

Peeking around you’ll see that the actual backing data for BitArray is an int[]:

private int[] m_array;

Likewise the get indexer boils down to this call:

return (m_array[index / 32] & (1 << (index % 32))) != 0;

And the relevant setter code is:

if (value) {
    m_array[index / 32] |= (1 << (index % 32));
} else {
    m_array[index / 32] &= ~(1 << (index % 32));
}

So inside the class there is a whole lot of shaking going, which is why a BitArray is slower than a bool[]. It’s also why an individual BitArray is not thread safe, i.e. don’t try to operate over partitioned ranges within the BitArray. Besides being unsafe, it may actually be slower than single-threaded (I tried).

However, a List, or in my code List, can be safer and faster because each thread works on an individual list item.

Concerns of High Volume of Operations (HVO)

Trust me: I don’t take pleasure in tweaking code to shave ½ a second here or there. I’d rather write my code to where the logic is clear to read and doesn’t cause the reader to pause with “Why is he doing it that way?” When I work with lists of 1M, 10M, or even 100 million items, I usually have my code written straightforward, and I may include certain safety checks. The problem with working with the sieve is that it requires a high volume of operations.

• An upperLimit of int.MaxValue requires access to a BitArray 3,315,151,693 times.
• An upperLimit of uint.MaxValue requires access to a BitArray 6,701,709,402 times.

Each little operation and each conditional check adds up. Group several together and you’ve added 1-3 seconds to the execution time. Likewise, if you can safely omit some, you can trim off some execution time.

Performance would slower if the Vector indexer checked if the index was in range. Vector is a private class and calls to the indexer are well controlled by me, so I omitted a range check.

Dividing by 2 for billions of times can be a drag. Instead I bit shift with >> 1.

Rather than rely upon implicit casting, which may obscure how many casts are being performed, I will use explicit casts in any spot I deem relevant. This is something that I would do extremely rarely in my code, if ever, other than working with sieves.

That said, there are spots of code where I am not so nitpicky. My concerns with performance kick in any time I must process over a BitArray. If I am setting up VectorList, and on my 8-core laptop there would be no more than 64 Vector’s, I don’t sweat over a performance tweak. I don’t want to shave nanoseconds for 64 operations. But I will do it for 3.3 billion.

Performance and Memory Usage

Performance is dependent upon the number of processor cores. My tests were done on an 8-core laptop. Memory results are taken from VS 2015 debugger to provide a general indication of memory usage. Execution time is more precise as it is measured with a Stopwatch object in Release mode.

For upperLimit of int.MaxValue:

• Sieve31 took 37.97 seconds and used 143 MB memory.

• EBrown’s answer took 22.18 seconds and used 1 GB memory.

• Sieve32Fast above took 20.38 seconds and used 149 MB.

For upperLimit of uint.MaxValue:

• Sieve32 took 77.10 seconds and used 283 MB memory.
• Sieve32Fast above took 41.23 seconds and used 293 MB memory.

Questions/Concerns

Magic numbers to const: Is it Overkill?

There are instances where I have chosen a value out of the blue. I name these as constants immediately before the one time I use them. See smallNumberCutoff or tinyFactor. The name helps add context to the value, and the fact that it’s a const is an indicator that it’s a magic number substitute. I think this is properly done.

But then there’s the constants Zero, One, Two, and Three. Are they magic numbers, or even if they aren't, is it overkill to use this?

I wonder if its overkill. Maybe my eyes are so use to reading 2 and 3 when dealing with primes. But it does seem helpful with the context that the value is a uint representing a Number versus an int representing an Index. So this tells me I am dealing with a context of uint Number:

var count = (UpperLimit / typicalNumberRange) + One;

whereas this tells me I am dealing int Index:

var stopIndex = _bits.Length – 1;

Visually one could say that I am inconsistent with +/- 1. But they are different, albeit in a very subtle way. The top of the class has very prominent comments to 'Sam the Maintainer' explaining this.

Should there be a Vector.Id property?

Even though it would have the exact same as its index in a VectorList? It works without it, but what it would provide is a means to replace this:

for (var vectorIndex = 0; vectorIndex < _vectors.Count; vectorIndex++)
{
    var startIndex = (vectorIndex == 0) ? rootBitIndex + 1 : 0;
    var vector = _vectors[vectorIndex];
    // more stuff
}

With this:

foreach (var vector in _vectors)
{
    var startIndex = (vector.Id == 0) ? rootBitIndex + 1 : 0;
    // more stuff
}

Should VectorList have a RootVector property?

It would simply return the item at index 0. Together with a Vector.Id, this would eliminate parent methods from having to reference a VectorList by index. It wouldn’t totally prevent it, nor do I think it should. But it would mean accessing VectorList becomes a foreach while stepping over a BitArray is a for (var bitIndex = 0; thing.

Using Three instead of rootVector.StartingNumber

The root vector is special. For the root vector, be it called rootVector or _vectors[0], I use code like:

var prime32 = ToNumber(bitIndex, Three);

From a purely technical standpoint, it would be ‘more proper’ to use:

var prime32 = ToNumber(bitIndex, rootVector.StartingNumber);

But for the root vector the StartingNumber will always be Three, and will never be anything but Three. I just feel better passing in the constant Three rather than fetching StartingNumber. While this is less proper, would it be totally improper?

Answer 1

Micro-Optimizations

One very micro-optimization you can make is, instead of subtracting numbers, add the negative variants.

For example:

var stopIndex = vector.BitLength - 1;

Is faster as:

var stopIndex = vector.BitLength + -1;

It's just the nature of the beast. Subtraction takes more cycles than addition, so if you add the negative complements of your values, you gain speed. (When calculating the first 100,000,000 primes I saw times drop from 3808ms to 3653ms, averaged over 5 runs.)

The same can be said about division and multiplication. Multiplying by the reciprocal is faster than dividing. (These are only beneficial if the calculation for the reciprocal happens fewer times than the multiplication.)

---

Another extremely micro-optimization you can make is to replace certain modulo operations with bitwise and operations.

Consider:

(this.UpperLimit & Two) == Zero)

This is slightly faster as:

((this.UpperLimit & 0x01U) == Zero)

This doesn't have much effect on the results, but I saw times drop further from 3653ms to 3622ms (5 run average).

---

Readability

Apart from the few micro-optimizations above, one of the issues I see with your code is just how packed it is.

For example, this is hard to read:

// output remaining primes 
for (var vectorIndex = 0; vectorIndex < _vectors.Count; vectorIndex++)
{
    var vector = _vectors[vectorIndex];
    var startIndex = (vectorIndex == 0) ? rootBitIndex + 1 : 0;
    // Due to high frequency of access, its ever so slightly faster to have copies created outside the loop
    // rather than called inside the loop directly and repeatedly with vector.BitLength and vector.StartingNumber.
    var length = vector.BitLength;
    var startingNumber = vector.StartingNumber;
    for (var bitIndex = startIndex; bitIndex < length; bitIndex++)
    {
        if (vector[bitIndex]) { yield return ToNumber(bitIndex, startingNumber); }
    }
}

It's ever so slightly easier to follow with a little more whitespace:

// output remaining primes 
for (var vectorIndex = 0; vectorIndex < _vectors.Count; vectorIndex++)
{
    var vector = _vectors[vectorIndex];
    var startIndex = (vectorIndex == 0) ? rootBitIndex + 1 : 0;
    // Due to high frequency of access, its ever so slightly faster to have copies created outside the loop
    // rather than called inside the loop directly and repeatedly with vector.BitLength and vector.StartingNumber.
    var length = vector.BitLength;
    var startingNumber = vector.StartingNumber;

    for (var bitIndex = startIndex; bitIndex < length; bitIndex++)
    {
        if (vector[bitIndex])
        {
            yield return ToNumber(bitIndex, startingNumber);
        }
    }
}

---

General Concerns

As far as your general questions:

But then there’s the constants Zero, One, Two, and Three. Are they magic numbers, or even if they aren't, is it overkill to use this?

The problem with having these constants is that the name of the constant doesn't imply what it means, but instead what it is. This means that any usage of them within the code is meaningless, as the only thing we know about them is that Zero is a uint valued at 0. If this is all the magic numbers are for, then having a constant adds absolutely no value to the code. It's only noise.

When you use the Find All References feature, it should be clear that each instance of that constant, property, field, method, class, struct, etc. is used in a specific scenario. However, these constants are not. So what is the point? I don't actually care about what references the constant Zero, so why is it a constant? Would I ever need to change it? No, so there's no value added.

Should there be a Vector.Id property?

Should VectorList have a RootVector property?

Using Three instead of rootVector.StartingNumber

All three of these have the same response:

Properties exist to allow is to idiomatically access information within them. By not using/creating properties, you fall into the trap of having to remember just exactly what the operation means. So, for:

var prime32 = ToNumber(bitIndex, Three);

Why did I choose Three for this? Ah, right, because that's what rootVector.StartingNumber is. However, now you have to either: add a comment explaining such, or remind yourself every time you get there what that value means. It adds no value to the code. You should use the property rootVector.StartingNumber, simply because it's self-explanatory. (It's also no slower than using the constant Three, which we should remove anyway as it adds no value to the code either.)

Generally speaking, you should use whatever practice adds the most value to your code, not whatever practice is the simplest. For example, you have the functions:

private static Func<uint, uint, int> ToIndex => (uint number, uint startingNumber) => (int)((number - startingNumber) >> 1);
private static Func<int, uint, uint> ToNumber => (int bitIndex, uint startingNumber) => (uint)(bitIndex << 1) + startingNumber;

Because of the value they add. Let's follow the same idioms throughout. :)

---

Best-Practice Concerns

Concerns about your BadUpperLimitException are pretty obvious as well. You have, unfortunately, fallen into a bad-practice here. Using a field to store all the Exception information, and merely throw _field; is a poor idea. Why? Because there are certain features you can no longer get.

private static ArgumentException BadUpperLimitException => new ArgumentException("upperLimit be must greater than or equal to 2.", "upperLimit");

public static IEnumerable<uint> Primes(int upperLimit)
{
    if (upperLimit < Two) { throw BadUpperLimitException; }
    return Primes((uint)upperLimit);
}

What happens of int upperLimit is refactored to a different name? Now the parameter in the ArgumentException is wrong.

However, if we rewrite it as so:

public static IEnumerable<uint> Primes(int upperLimit)
{
    if (upperLimit < Two)
    {
        throw new ArgumentException($"{nameof(upperLimit)} be must greater than or equal to 2.", $"{nameof(upperLimit)}");
    }

    return Primes((uint)upperLimit);
}

We can refactor int upperLimit with no issues. (This may not seem like a benefit, but it adds significant value to the code in the form of idiomatic code. It's much more clear what is going on with this example, than with your original code.)