How to generate a couple million prime numbers in C#

Introduction

This is the first in what I am hoping is a series of articles describing some of my coding for Cryptography and Prime Numbers. You can find all of this on the Internet (somewhere) but I am trying to give a newbie a jumping off point. One of the first things you will need for Cryptography and Primality Testing is some prime numbers. A thousand is a good start but a million or more is better. And once the bugs are out, I trust the PC and code more than my own ability to perform perfect data entry for a million Primes. And that is a lot of work. The other option is to buy such a list but I am cheap. And there are free lists available but there are format issues there. (Plus, there might be an error in the file. Unlikely but possible.)

Originally I thought I could just start at 7, take its square root, check for divisors from 2 up to that square root. If no divisors found, the number is prime and put it into a List. If a divisor is found, skip to the next number and repeat. This is fine for a single number that is not too large or for a list of a couple thousand Primes (I suppose), but it gets slow waiting for a lot of numbers to be checked in this way.

Late at night, I read on some Cryptography board a very smart idea someone described. Set an upper bound for the highest number you want to go to, and then starting at the lowest Prime, mark the square of that number as being not Prime and march up to the highest number by increments of twice the Prime. When that lowest Prime has finished, go to the next unmarked Prime and repeat the squaring and setting loop. (If that seems complicated, look at the pictures below. Pictures are how I make sense out of complex problems.)

As for what you get for your speed of running... it will depend on your PC (CPU speed and memory size and speed, how much you have running at the same time, etc.). If the code runs too slow, lower the iMax number.

Background

Avoid testing and setting all the numbers that are cleanly divisible by 2, 3, and 5. This can be done by taking the Modulus 30 of each number. (And it reduces the amount of numbers that need to be examined to 8 out of 30.) A nice graphic of this idea can be seen in Rick Oden's article (here). Imagine a table of 30 columns reduced to remove all columns with even numbers as well as those divisible by by 3 or 5. The first 6 rows of this are:

Every one of those numbers is potentially Prime. For the sake of argument, skip the number "1". I just put it into the table because otherwise there is an empty space there and the top row looks aesthetically pleasing because it holds "N % 30" for every cell beneath it. (Some smart person is going to object: but "1" is not Prime! I say, chill. It's a special case. 1 is the Identity for Multiplication. Just let it slide for now.)

So starting with the number "7", square that number to get "49". We mark that as not Prime. And we also advance by 7 (which is even and not on the table) to "56". Skipping that, we go 7 more and reach "77". And to save time, we will not try to advance by 7. (An odd plus an odd number is always an even (and not on our table). An odd number and twice the Prime in question will be odd and may be on the table.) We add 14 to get to the next number that has a factor of 7, which is "91"... and on we go doing that, gaps of 14, and marking that number as not Prime. If a number is not on our table, we keep going anyway until we reach the maximum number. And our table will look like this:

The next unmarked number in the table is "11". Its square is "121". Mark it not Prime. Add 22 and get to "143", mark it not Prime. And on and on by steps of 22...

If we advance to a marked cell, it's composite and can be skipped. ("49" is our first such cell. And we skip it. The next cell after it is 53, square it and set that, and increment by 106 and it's back to business as usual.)

To code this, allocate a large array of bytes to work in. Bytes_Allocated * 30 = N = Maximum Prime Number to march up to. Each byte in the array corresponds to N / 30. Each bit in a byte corresponds to the 8 possible values of N % 30. (1, 7, 11, 13, 17, 19, 23, 29) for 8 possible values. 1 per bit. So, if we allocate a 1000 bytes, we can cover the Primes up to the number 30000. When the code finishes, any bit left unset will represent a Prime. In this case, I am adding each found Prime to a List.

Using the Code

The following C# code can be cut and paste into Visual Studio. Tie it into your own button click. After that, all the usual - use the code at your own risk and I am not responsible for your possibly getting hooked on chasing Prime Numbers or Primality Testing.

There are no fancy classes here. Just going for speed with some globals, passing by reference, and the optimizations described above for what to even consider Prime.

List<UInt32> listPrimes = new List<UInt32>();
private void button8_Click(object sender, EventArgs e)
{
    UInt32 iBump = 2;
    UInt32 iPrime = 0;
    UInt32 iPrime2 = 0;
    UInt32 iPrimeSqr = 0;
    UInt32 iMax   = 20000;      // 0 sec
    iMax         = 200000;      // 0.2 sec
    iMax         = 2000000;     // 0.8 sec
    iMax         = 20000000;    // 8.3 sec; gens List of first 8.8M primes
    //iMax         = 200000000; // 103 sec; gens List of first 70M primes
    UInt32 iMax8  = iMax * 8;
    UInt64 iiMax8  = iMax * 8;
    if ( iiMax8 > UInt32.MaxValue )
    {
        System.Diagnostics.Debug.WriteLine ( 
          "Too many numbers requested for this data type!");
        return;
    }

    // Start out all = false
    byte[] aPrimesBool = new byte [ iMax ];

    UInt32 iFactor = 0;

    DateTime dt0 = DateTime.Now;
    DateTime dt1 = DateTime.Now;

    // Clear the old primes list
    listPrimes.Clear();
    listPrimes.Add ( 2 );
    listPrimes.Add ( 3 );
    listPrimes.Add ( 5 );

    // Every odd number is potentially prime; set every byte (8 bits at a time)
    for ( iFactor = 0; iFactor < iMax; iFactor += 1 )
    {
        aPrimesBool [ iFactor ] = 0xff;
    }
    dt1 = DateTime.Now;
    System.Diagnostics.Debug.WriteLine ( "start after : " + 
            ( dt1 - dt0).TotalSeconds + " sec");

    // Init _iGR and i_GM to 7
    _iGW = 0; // == 7 / 30
    _iGR = 7; // == 7 % 30
    for ( iPrime = 7; iPrime < iMax8; iPrime += iBump )
    {
        if ( iPrime % 1000000 == 1 )
        {
            dt1 = DateTime.Now;
            System.Diagnostics.Debug.WriteLine ( "reached: " + 
              iPrime + " after : " + ( dt1 - dt0).TotalSeconds + " sec");
        }

        if ( GetCodePrime ( ref aPrimesBool ) ) // iPrime passed via _iGR and _iGW
        {
            listPrimes.Add ( iPrime );
                    
            iPrime2   = iPrime * 2;
            iPrimeSqr = iPrime * iPrime;

            // Every multiple of this prime is not prime, loop and mark them as not prime
            for ( iFactor = iPrimeSqr; iFactor < iMax8; iFactor += iPrime2 )
            {
                ClrCodePrime ( ref iFactor, ref aPrimesBool );
            }
        }

        // Alternate upping to the next possible prime by 2 and 4
        // ( 7 + 4, 11 + 2, 13 + 4, 17 + 2, 19 + 4,
        //        23 + 2, 25 + 4, 29 + 2, 31 + 4, 35 + 2, ... )
        if ( iBump == 2 )
            iBump = 4;
        else
            iBump = 2;

        // Counters for GetCodePrime()
        if ( _iGR < 29 )
        {
            _iGR += iBump; // up the remainder
        }
        else
        {
            _iGR = 1;
            _iGW++;
        }
    }

    dt1 = DateTime.Now;
    System.Diagnostics.Debug.WriteLine ( "done: " + iPrime );
    System.Diagnostics.Debug.WriteLine ( "count = " + listPrimes.Count );

    // wipe out the byte array
    aPrimesBool = null;
    System.GC.WaitForPendingFinalizers();
}


UInt32 _iCR = 0; // Clear Remainder
UInt32 _iCW = 0; // Clear Whole
private void ClrCodePrime ( ref UInt32 iNum, ref byte[] bA )
{
    _iCR = iNum % 30;
    switch ( _iCR )
    {
        case  1: 
            _iCW = iNum / 30; // Only calculate
              // the whole number if the remainder is worth it
            bA [ _iCW ] &= 0xfe; // Clear bit 2^0 == 1
            return;

        case  7: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0xfd; // Clear bit 2^1 == 2
            return;

        case 11: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0xfb; // Clear bit 2^2 == 4
            return;

        case 13: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0xf7; // Clear bit 2^3 == 8
            return;

        case 17: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0xef; // Clear bit 2^4 == 16
            return;

                
        case 19: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0xdf; // Clear bit 2^5 == 32
            return;

                
        case 23: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0xbf; // Clear bit 2^6 == 64
            return;

                
        case 29: 
            _iCW = iNum / 30;
            bA [ _iCW ] &= 0x7f; // Clear bit 2^7 == 128
            return;
    }

    // if remainder is anything else, ignore the number because
    // it is divisible by 2, 3, or 5.
    return;
}

UInt32 _iGR = 0; // Get remainder
UInt32 _iGW = 0; // Get whole
private bool GetCodePrime ( ref byte[] bA )
{
    //
    // Returns true if ( _iGR  + 30 * _iGW ) is prime
    //
    switch ( _iGR )
    {
        case  1: // 1 + 30 * _iGW 
            return ( bA[_iGW] & 0x01 ) == 0x01; 

        case  7: // 7 + 30 * _iGW 
            return ( bA[_iGW] & 0x02 ) == 0x02;

        case 11: // 11 + 30 * _iGW 
            return ( bA[_iGW] & 0x04 ) == 0x04;

        case 13: // 13 + 30 * _iGW 
            return ( bA[_iGW] & 0x08 ) == 0x08;

        case 17: // 17 + 30 * _iGW 
            return ( bA[_iGW] & 0x10 ) == 0x10;

                
        case 19: // 19 + 30 * _iGW
            return ( bA[_iGW] & 0x20 ) == 0x20;

                
        case 23: // 23 + 30 * _iGW
            return ( bA[_iGW] & 0x40 ) == 0x40;

                
        case 29: // 29 + 30 * _iGW 
            return ( bA[_iGW] & 0x80 ) == 0x80;
    }

    return false;
}

I have shied away from multi-threading to make it a little easier to understand. If you find a way to make this code run faster, by all means please share what you changed.

Points of Interest

Some cute bit clearing is going on in the ClrCodePrime routine. And because it is not good classes, be careful when changing stuff. It is very possible to break something accidentally. (Part of the whole thing with why classes are better than this dinosaur style of code. I maintain that the old lizard runs faster. But if you have a faster class, share and we can all decide if it is better.)

With generating a lot of Primes, it is less efficient to compute each Prime individually (and look beneath it to see that it has no Factors) and more efficient for each Prime to set all of the factorable Composites above it to some requested Maximum.

History

• Initial article = 9/16/2011.