5
\$\begingroup\$

I wrote yet another optimized single-threaded Eratosthenes Sieve implementation in C:

erato.c

// #include <stdio.h>
#include <stdlib.h>

#define N 1000000000
#define num_t unsigned long

int main() {
    register char *b = malloc(N * sizeof(char));
    for(num_t i = 0; i ^ N; ++i)
        b[i] = !(i & 1);

    // printf("2\n");
    for(num_t i = 3; i ^ 1 ^ N; i += 2) {
        if(!b[i - 2]) {
            // printf("%llu\n", i);
            const num_t increment = i << 1;
            num_t j = i;
            while(j < N) {
                b[j - 2] = 1;
                j += increment;
            }
        }
    }
    free(b);
    return 0;
}

The output is omitted because it takes most of the time to process.

Optional flags: -O3

Performance

The sieve covers first 1.000.000.000 numbers in

11.14s user 0.32s system 99% cpu 11.465 total

Could you please help me to improve my programming/mathematical skills by reviewing this piece of code?

\$\endgroup\$
10
\$\begingroup\$

Bug

At the end of your program, the b array is all full of 1s, which means that you didn't find any primes. The problem is here:

       num_t j = i;

Because you start your j loop at i, you will mark i (which is prime) as non-prime. You should start your j loop at i*i instead (see below).

Don't get tricky

Why write this:

for(num_t i = 0; i ^ N; ++i)

when you can write this:

for(num_t i = 0; i < N; ++i)

I had to stare at your code for a long time just to figure out whether it was correct, for a simple loop!

Also, you did it here:

for(num_t i = 3; i ^ 1 ^ N; i += 2) {

which should be:

for(num_t i = 3; i < N; i += 2) {

You can afford to waste 2 bytes

Here, you use index-2 to save 2 bytes:

    if(!b[i - 2]) {
        // printf("%llu\n", i);
        const num_t increment = i << 1;
        num_t j = i;
        while(j < N) {
            b[j - 2] = 1;
            j += increment;
        }
    }

There's no need to save 2 bytes. Just use b[index] instead of b[index-2]. You already allocated the full array anyways.

Speed improvements

There are three things I can see that could speed up your program:

  1. In your j loop, you can start j at i*i instead of just i. All non-primes less than i*i will already be marked.
  2. You can stop the i loop at sqrt(N) instead of at N, because once you reach sqrt(N), the second loop won't mark any more entries.
  3. Since you are finding primes up to 10^9, you are using 1GB of memory for the b array. You should get better performance if you use less memory, because the processor cache will be utilized more effectively if your array is smaller.

    a) You can only track odd numbers, which will reduce your memory usage to 1/2 of the original (512MB).

    b) You can use 1 bit per number instead of 1 byte per number. This will reduce your memory to 1/8 of the original. Combined with (a), it will be 1/16 the original, or 64MB. This is a much better than 1GB.

Example using 1 bit per odd number

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define N        1000000000

int main(void)
{
    int       arraySize = ((N + 63)/64 + 1) * sizeof(uint32_t);
    uint32_t *primes    = malloc(arraySize);

    // Each bit in primes is an odd number.
    // Bit 0 = 1, bit 1 = 3, bit 2 = 5, etc.
    memset(primes, 0xff, arraySize);

    // 1 is not a prime.
    primes[0] &= ~0x1;

    int sqrt_N = sqrt(N);
    for(int i = 3; i <= sqrt_N; i += 2) {
        int iIndex = i >> 6;
        int iBit   = (1 << ((i >> 1) & 31));
        if ((primes[iIndex] & iBit) != 0) {
            int increment = i+i;
            for (int j = i * i; j < N; j += increment) {
                int jIndex = j >> 6;
                int jBit   = (1 << ((j >> 1) & 31));
                primes[jIndex] &= ~jBit;
            }
        }
    }

    // Count the number of primes in order to verify that the above worked.
    // Start count at 1 to include 2 as the first prime, since we are only
    // going to count odd primes.
    int count = 1;
    for (int i = 3; i < N; i += 2) {
        int iIndex = i >> 6;
        int iBit   = (1 << ((i >> 1) & 31));
        if (primes[iIndex] & iBit)
            count++;
    }
    printf("%d\n", count);

    free(primes);
    return 0;
}

Edit: Even less memory usage

Pete Kirkham suggested using even less memory by using 2 out of every 6 numbers. In other words, using 1 bit per every 3 numbers instead of 1 bit per every 2 numbers. At first I was skeptical because this required using a division in the inner loop. However, after coding it up, it turned out to be faster. The code is quite a bit trickier however, because the inner loops need to avoid any multiples of 3, because all multiples of 3 are no longer in the primes array:

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define N        1000000000

int main(void)
{
    int       arraySize = (N/24 + 1);
    uint32_t *primes    = malloc(arraySize);

    // The bits in primes follow this pattern:
    //
    // Bit 0 = 5, bit 1 = 7, bit 2 = 11, bit 3 = 13, bit 4 = 17, etc.
    //
    // For even bits, bit n represents 5 + 6*n
    // For odd  bits, bit n represents 1 + 6*n
    memset(primes , 0xff, arraySize);

    int sqrt_N = sqrt(N);
    for(int i = 5; i <= sqrt_N; i += 4) {
        int iBitNumber = i / 3 - 1;
        int iIndex = iBitNumber >> 5;
        int iBit   = 1 << (iBitNumber & 31);
        if ((primes[iIndex] & iBit) != 0) {
            int increment = i+i;
            for (int j = i * i; j < N; j += increment) {
                int jBitNumber = j / 3 - 1;
                int jIndex = jBitNumber >> 5;
                int jBit   = 1 << (jBitNumber & 31);

                primes[jIndex] &= ~jBit;

                j += increment;
                if (j >= N)
                    break;

                jBitNumber = j / 3 - 1;
                jIndex = jBitNumber >> 5;
                jBit   = 1 << (jBitNumber & 31);

                primes[jIndex] &= ~jBit;

                // Skip multiple of 3.
                j += increment;
            }
        }
        i += 2;
        iBit <<= 1;
        if ((primes[iIndex] & iBit) != 0) {
            int increment = i+i;
            for (int j = i * i; j < N; j += increment) {
                int jBitNumber = j / 3 - 1;
                int jIndex = jBitNumber >> 5;
                int jBit   = 1 << (jBitNumber & 31);

                primes[jIndex] &= ~jBit;

                // Skip multiple of 3.
                j += increment;

                j += increment;
                if (j >= N)
                    break;

                jBitNumber = j / 3 - 1;
                jIndex = jBitNumber >> 5;
                jBit   = 1 << (jBitNumber & 31);

                primes[jIndex] &= ~jBit;
            }
        }
    }

    // Initial count includes 2, 3.
    int count=2;
    for (int i=5;i<N;i+=6) {
        int iBitNumber = i / 3 - 1;
        int iIndex = iBitNumber >> 5;
        int iBit   = 1 << (iBitNumber & 31);
        if (primes[iIndex] & iBit) {
            count++;
        }
        iBit <<= 1;
        if (primes[iIndex] & iBit) {
            count++;
        }
    }
    printf("%d\n", count);

    free(primes);
    return 0;
}

I'm sure the variant that works on multiples of 30 would be even faster than this, although the code would be even more complicated.

Timings

Here are the speeds of the various programs I ran on my computer:

Original program        : 13.33 seconds
OP fixed (i*i, sqrt(n)) :  7.75 seconds
1 bit per odd number    :  3.90 seconds
1 bit per 3 numbers     :  2.77 seconds

Note that I did not include the time it took to count the primes to verify the result. I modified each program to count the number of primes just to make sure that they all worked, which is how I discovered that the original program had a bug in it.

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ It's worth trying i * i <= N rather than i <= sqrt_N to avoid the sqrt() call \$\endgroup\$ – Tavian Barnes Dec 4 '15 at 21:18
  • \$\begingroup\$ Wow! This optimization is impressive! I was sure my code was perfect, but something was wrong with it. That something is huge this time. But isn't i < N slower than i ^ 1 ^ N? This is a repeated operation. \$\endgroup\$ – theoden8 Dec 4 '15 at 21:43
  • \$\begingroup\$ @TavianBarnes, I guess i * i will be carried out 1000000000 times, whereas sqrt will be performed only once. \$\endgroup\$ – theoden8 Dec 4 '15 at 21:45
  • 1
    \$\begingroup\$ @theoden going from 2 to 6 means a factor of three in the outer loop, and also makes the array more compressed, so you are using less memory bandwidth. 30 would would make the calculations more complicated, but sometimes you can trade off by only using 30 bits of your words, and your initial memset is replaced by a for loop marking off every fifth bit. \$\endgroup\$ – Pete Kirkham Dec 4 '15 at 22:55
  • 1
    \$\begingroup\$ @PeteKirkham I tried the 2 of 6 way and it was indeed faster than the odd number way. I'll post the code and results above. \$\endgroup\$ – JS1 Dec 4 '15 at 23:47
6
\$\begingroup\$
  1. Don't use register: Most likely it won't have any effect (aside from forbidding address-of on the variable).
    And if it does have an effect, did you actually profile it?

  2. sizeof(char) is always 1, which makes using it pointless. Anyway, even if it was a different type you should use sizeof expression.

  3. You didn't check whether malloc failed. Why are you using it anyway, instead of a statically allocated array?

  4. i ^ N for i != N? That isn't a submission to the IOCCC, so I'm a bit puzzled why you would do that...

  5. You could save half the space for your big array if you omitted all even numbers. And reduce that to an eightth by bit-packing, though that need not be an advantage, profile that.

  6. The body of your if-statement should be a single for-loop (incorporating bug-fix from JS1):

    for(num_t increment = i << 1, j = i * i; j < N; j += increment)
        b[j - 2] = 1;
    
  7. Since C99, return 0; is implicit for main.

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ Sometimes it is very interesting to understand how bad your code is :) (3) Allocation in stack memory leads to size restrictions, actually. (4) Isn't i ^ N faster than i != N? This operation is carried out 1000000000 times. (5) Bit-packing is a good idea, wouldn't that force to use own-written operators which would be slower than standard ones? \$\endgroup\$ – theoden8 Dec 4 '15 at 21:37
  • 2
    \$\begingroup\$ I said static, not automatic. And while i ^ N is one character less than i != N, it should lead to the exact same code. Yes, bit-packing means you have to mask and shift on each access (read as well as write), which means more code. So the question is whether the extra instructions or the reduced memory-demand weigh heavier. \$\endgroup\$ – Deduplicator Dec 4 '15 at 21:43
  • \$\begingroup\$ Is alloca is that you mean? It works for me until N is more than 10 ^ 8. Pretty much, though not sufficient for 10^9. \$\endgroup\$ – theoden8 Dec 4 '15 at 23:09
  • \$\begingroup\$ No, static allocation is things like globals, and locals which are marked static. alloca is stack-allocation too, just like automatic variables. \$\endgroup\$ – Deduplicator Dec 5 '15 at 0:42
  • 1
    \$\begingroup\$ @cad: I'm saying sizeof expression, instead of sizeof(type). Whether you decide to add additional parentheses to the expression or not is your decision. \$\endgroup\$ – Deduplicator Dec 5 '15 at 15:04
4
\$\begingroup\$

Only minor things to add as previous reviewers did a fine job.

The following would be an infinite loop in N was odd. Since the definition of N did not comment of this limitation, code should not assume N is even.

for(num_t i = 3; i ^ 1 ^ N; i += 2) {

Mis-statement: This code does not generate "1.000.000.000 prime numbers in" so many seconds. Instead it generate primes with values less than 1.000.000.000.

#define num_t unsigned long is interesting. Yet since the type is used as an array index, any type wider that size_t is a problem and unsigned long may be wider than size_t Consider a comment near the definition to that effect or some compile time test.

| improve this answer | |
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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