RNR's blog

Help decoding Yarin sieve
By RNR, history, 5 years ago, In English

I know the standard sieve algorithm for finding primes in a range [0..upperbound]

Here is a standard implementation

for(int i=2;i<=N;i++) 
    if(prime[i])
        for(int j=i*i;j<=N;j+=i) prime[j]=0;

But when memory is tight, here is the famous implementation by yarin

#define MAXSIEVE 100000000 // All prime numbers up to this
#define MAXSIEVEHALF (MAXSIEVE/2)
#define MAXSQRT 5000 // sqrt(MAXSIEVE)/2
char a[MAXSIEVE/16+2];
#define isprime(n) (a[(n)>>4]&(1<<(((n)>>1)&7))) // Works when n is odd

int i,j;
memset(a,255,sizeof(a));
a[0]=0xFE; // 254
for(i=1;i<MAXSQRT;i++)
  if (a[i>>3]&(1<<(i&7)))
    for(j=i+i+i+1;j<MAXSIEVEHALF;j+=i+i+1)
      a[j>>3]&=~(1<<(j&7));

I understood that we are only storing odd numbers and char can store 8 bits and we are storing it using j>>3 and unsetting using ~(1<<(j&7)). When considering only odd numbers, you let index 1 be 3, index 2 be 5 etc.

But I didn't understand the strange loop incrementations and initialization in for(j=i+i+i+1;j<MAXSIEVEHALF;j+=i+i+1)

Any help would be greatly appreciated. Thank you in advance.

  • Vote: I like it
  • 0
  • Vote: I do not like it

»
5 years ago, # |
← Rev. 2   Vote: I like it +6 Vote: I do not like it

First thing to notice here is that the check a[i>>3]&(1<<(i&7)) isn't about the primality of $$$i$$$, but the primality of $$$2i + 1$$$. So we can rewrite our code like this:

memset(a, 255, sizeof(a));
a[0] = 0xFE; // 254

for (int i = 1; i < MAXSQRT; i++)
  if (isprime(2 * i + 1))
    for (int j = 3 * i + 1; j < MAXSIEVEHALF; j += 2 * i + 1)
      mark 2 * j + 1 as not prime;

Now let's do the substitutions $$$u = 2i + 1$$$ and $$$v = 2j + 1$$$.

for (int u = 3; u < sqrt(MAXSIEVE); u += 2)
  if (isprime(u))
    for (int v = 3 * u; v < MAXSIEVE; v += 2 * u)
      mark v as not prime;

Now this is quite standard. The only thing to notice here is that $$$u$$$ and $$$v$$$ both skip the even numbers, that's why you see u += 2 and v += 2 * u, contrary to the standard implementation.

In case you're wondering about a[0], it's set to 254 to exclude 1 from the primes as a special case.

In fact I don't get why the author didn't write the code like I did. The code you give is basically unreadable, littered with crap like <<<>>0&1<<4>>(n). You might be thinking that it is somehow slower to call isprime(u) and triggering a division by 2, but modern compilers are really smart. They can do some really non-trivial optimizations. By the way the same goes for the "i++ is slower than ++i" crowd. At least if you use them as the increment in a simple for-loop, the compiler is smart enough to realize that you don't care about it's return value.

→ Reply
  • »
    »
    5 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    Thank you -is-this-fft-, I mean it.

    We can make it faster by initializing v=u*uthat becomes j=2*i*i+2*i

    → Reply