Efficient sieving implementation?

[Oswald Kluge]

Dear Reader,

I'm trying to implement the following short algorithm:

Given a number n, remove all the multiples of n
from the list of non-negative numbers from 1 through a limit k.

My solution is to take an array with all entries '1' and iterate
over it in steps of n, setting the current pointer to '0'.

Here's the code:

const LIMIT=k+1;
unsigned char Sieve[LIMIT];

memset(Sieve,1, LIMIT); // initialise all entries to '1'
for(i=n;i<LIMIT ;i+=n) {Sieve[i]=0;}

As k gets large this implementation is really slow for small n,
especially if you have to sieve more than one n. Is there a
more efficient way to implement this? Are there fast memory
operators like memset in c++ for this task?

Thanks for your help!

Best,
Oswald Kluge

[Tom St Denis]

Oswald Kluge wrote:

Dear Reader,

I'm trying to implement the following short algorithm:

Given a number n, remove all the multiples of n
from the list of non-negative numbers from 1 through a limit k.

My solution is to take an array with all entries '1' and iterate
over it in steps of n, setting the current pointer to '0'.

Here's the code:

const LIMIT=k+1;
unsigned char Sieve[LIMIT];

memset(Sieve,1, LIMIT); // initialise all entries to '1'
for(i=n;i<LIMIT ;i+=n) {Sieve[i]=0;}

As k gets large this implementation is really slow for small n,
especially if you have to sieve more than one n. Is there a
more efficient way to implement this? Are there fast memory
operators like memset in c++ for this task?

Sure, instead of

if (Sieve[blah] == 0) { happy_dance(); }

Do

if (Sieve(blah) == 0) { happy_dance(); }

Where Sieve() is a function that returns the value desired.

If your domain is truly large than a table is just inefficient and
wasteful. If the table is really small then your concern is moot.

Tom

[pegguru]

Oswald Kluge wrote:

Dear Reader,

I'm trying to implement the following short algorithm:

Given a number n, remove all the multiples of n
from the list of non-negative numbers from 1 through a limit k.

My solution is to take an array with all entries '1' and iterate
over it in steps of n, setting the current pointer to '0'.

Here's the code:

const LIMIT=k+1;
unsigned char Sieve[LIMIT];

memset(Sieve,1, LIMIT); // initialise all entries to '1'
for(i=n;i<LIMIT ;i+=n) {Sieve[i]=0;}

As k gets large this implementation is really slow for small n,
especially if you have to sieve more than one n. Is there a
more efficient way to implement this? Are there fast memory
operators like memset in c++ for this task?

Thanks for your help!

Best,
Oswald Kluge

It may be possible to improve this for a generalized case (it can be
done for specific cases of multiple "n"s), but it will certainly
involve more memory.

I know of nothing like memset that would be useful in this case. When
you get down to the machine level, any processor I've ever seen (and
I've seen many) will still need to do a loop like your for loop. A
good compiler will turn the for loop into efficient machine code - it's
a very simple bit of code.

Can you provide more information about the range of values for which
you are doing this, the number and types of "n"s, programmatic
constraints on cpu time and memory availability, and the overall
purpose? If you are doing this to find primes, etc, this is not likely
the best route.

Jeff

[Christian Siebert]

I'm trying to implement the following short algorithm:

>
Given a number n, remove all the multiples of n
from the list of non-negative numbers from 1 through a limit k.

Seems to be the "Sieve of Eratosthenes", or at least a very similar
algorithm ...

My solution is to take an array with all entries '1' and iterate
over it in steps of n, setting the current pointer to '0'.

That might be what you want (input and output). On the other hand, a
direct method (like Tom already suggested) can be much faster. And even
if you need this table-driven approach, this usual iteration is not the
only way to achieve this goal.

As k gets large this implementation is really slow for small n,

This is because memory is optimized for streaming access and such
pseudo-random pattern isn't cache-friendly.

especially if you have to sieve more than one n.

Instead of sieving for each n separately, try to combine several passes
for different n's and sieve in blocks of your processors cache size.
Depending on the number of n's, this can give you good performance
boost.

Christian

[bachelor]

it is worth caching if you will call Sieve(x) on the same x several
times. Otherwise don't cache. If you are just not sure try this:
_______________ _______________ _____
const int K = 1000000;
unsigned char sieve[K];

const int N = 3;
int divs[N] = {3, 7, 11};

char div = 0, notdiv = 1, unknown = 2;

void init_Sieve(){ memset(Sieve, unknown, K); }

inline Sieve(int x)
{
if( sieve[x] == unknown)
{
for(int i=0; i<N; i++)
if(x%divs[i] == 0){
sieve[x] = div;
break;
}
}
return sieve[x];
}
_______________ _______________ _________
This way you get sure you don't calculate twice for the same x, with a
trade off -- another couple of CPU cycles for each call Sieve(x) BUT
never calculate Sieve(x) for any never-used x.

[Oswald Kluge]

Can you provide more information about the range of values for which

you are doing this, the number and types of "n"s, programmatic
constraints on cpu time and memory availability, and the overall
purpose? If you are doing this to find primes, etc, this is not likely
the best route.

Jeff

The purpose of this program is to find coprime numbers. If n is prime
then
if you sieve out all multiples of n, then all remaining numbers will
have no
common divisors with n.

If n is composite, I always know what prime factors are occurring, so
after
sieving out all the prime factors you still get all the coprime
numbers.

The range of values starts at k=10^5. I got a pentium 2 GHz with 2GB
RAM.

Oswald

[Tom St Denis]

Oswald Kluge wrote:

Can you provide more information about the range of values for which
you are doing this, the number and types of "n"s, programmatic
constraints on cpu time and memory availability, and the overall
purpose? If you are doing this to find primes, etc, this is not likely
the best route.

Jeff

The purpose of this program is to find coprime numbers. If n is prime
then
if you sieve out all multiples of n, then all remaining numbers will
have no
common divisors with n.

Why not just compute the GCD of the two numbers?

Tom

[pegguru]

Oswald Kluge wrote:
....snip...

The purpose of this program is to find coprime numbers. If n is prime
then
if you sieve out all multiples of n, then all remaining numbers will
have no
common divisors with n.

If n is composite, I always know what prime factors are occurring, so
after
sieving out all the prime factors you still get all the coprime
numbers.

The range of values starts at k=10^5. I got a pentium 2 GHz with 2GB
RAM.

Oswald

I was really looking for a range of values (i.e. both ends). If the
cardinality of the range is not too large, sieving might be worthwhile.
I'm uncertain why you would need to know all of the coprimes of n
within a given range.

If the range of values is not too large, and if you have to do a LOT of
sieving, it might be possible to set up a network in memory, wherein
each number in the range would be connected to its prime factors. To
eliminate numbers with a common factor to n, you would only have to
specify the factors of n itself. This concept would require a fair bit
of precomputation and would use a lot of memory for overhead, so it
would not be useful if your range was, say, hundreds of millions of
numbers.

Are you really certain that you need to know all of the coprimes of n
within a given range? If not, Tom's suggestion of calculating the GCD
of n and a specific other number is excellent.

Jeff

[Logan Shaw]

Christian Siebert wrote:

>I'm trying to implement the following short algorithm:

Given a number n, remove all the multiples of n
from the list of non-negative numbers from 1 through a limit k.

>As k gets large this implementation is really slow for small n,

This is because memory is optimized for streaming access and such
pseudo-random pattern isn't cache-friendly.

>especially if you have to sieve more than one n.

Instead of sieving for each n separately, try to combine several passes
for different n's and sieve in blocks of your processors cache size.

One way to do this, though it may be obvious to some, is to combine n's
into groups by taking the least common multiple. That gives you a pattern
you can repeat. For example, when doing a sieve with 2 and 3, you can do
the straightforward way:

for (i = 0; i < end; i += 2) { a[i] = 0; }
for (i = 0; i < end; i += 3) { a[i] = 0; }

Or instead, you can (mostly) combine these into one pass:

for (i = 0; i < end; i += 6)
{
/* multiples of 2 */
a[i + 0] = 0;
a[i + 2] = 0;
a[i + 4] = 0;

/* multiples of 3 */
a[i + 0] = 0;
a[i + 3] = 0;
}

This can be minimized by combining the redundant assignments:

for (i = 0; i < end; i += 6)
{
a[i + 0] = 0;
a[i + 2] = 0;
a[i + 3] = 0;
a[i + 4] = 0;
a[i + 5] = 0;
}

Now, here's what I think is the fun part. Once you've established
that pattern that repeats after 6 elements, you can, instead of
performing all the assignment statements separately, create a mask:

int lcm = 6;
bool mask[lcm];

for (i = 0; i < lcm; i++) { mask[i] = 1; }

for (i = 0; i < lcm; i += 2) { mask[i] = 0; }
for (i = 0; i < lcm; i += 3) { mask[i] = 0; }

and then apply the mask:

for (i = 0; i < end; i += lcm)
{
for (j = 0; j < lcm; j++)
{
a[i + j] &= mask[j];
}
}

The good thing about this, obviously, is that the inner loop in the
code just above can be optimized to DEATH. Your operating system
or some readily available library may in fact already provide a
routine to AND two regions of memory together efficiently. You
can duplicate the mask until it is some convenient size. You can
make its size a multiple of 8 or 16 or 32 if that helps you align
it in a useful way. And you can combine multiple different n's
together to produce one mask.

Of course, all this is a waste of time if the n's aren't relatively
small compared to k.

- Logan

[Christian Siebert]

Hi Logan,

Instead of sieving for each n separately, try to combine several passes
for different n's and sieve in blocks of your processors cache size.

One way to do this, though it may be obvious to some, is to combine n's
into groups by taking the least common multiple. That gives you a pattern
you can repeat.

nice that you think of this improvement too. I posted this idea some
months ago in conjunction with the QS:
http://groups.google.de/group/sci.cr...90676702d754ba
which gives you an additional speed-up. Though I meant the "more
common" optimization technique in the first place:

Instead of using the whole large sieve at once

memset(Sieve,1, LIMIT); // initialise all entries to '1'
for each n do
for (i=n; i<LIMIT; i+=n) { Sieve[i]=0; }
use_sieve(Sieve , LIMIT);

it is often better to split this large sieve into smaller blocks

start = 0;
while (start < LIMIT) {
memset(SieveBlo ck,1,CACHESIZE) ;
for each n do
for (i = ceil(start/n)*n; i < (start+CACHESIZ E); i+=n)
{ SieveBlock[i-start]=0; }
use_sieve_block (SieveBlock, start, CACHESIZE);
start += CACHESIZE; // proceed with next block
}

Now, here's what I think is the fun part. Once you've established
that pattern that repeats after 6 elements, you can, instead of
performing all the assignment statements separately, create a mask:
... and then apply the mask: ... & ...
Your operating system or some readily available library may in fact
already provide a routine to AND two regions of memory together efficiently.

Why should you waste an additional logical operation, when it is
possible to shift this work into the initialization step. Create the
"mask" like suggested, and then use memcpy() to repeatedly copy it to
your sieve memory...

Christian