Re: Vernon's Prime Sieve

From: vernonner3voltazim (vnemitz_at_pinn.net)
Date: 07/29/04 

Date: 28 Jul 2004 22:20:58 -0700

"Daniel McLaury" <daniel_mcl@hotmail.com> wrote:

> I'll put in a breakpoint later today, after I watch the botball
> nationals.

OK, Thanks!

> I did a quick memory optimization of the sieve program,
> getting this:

Yes, one of the things I was thinking about was reusing
the buffer holding candidates for sieving. I've inserted
a couple comments below
 
> /** Calculate all primes up to M * N
> *** Silly preconditions:
> *** pi(N) < P
> *** M < N
> *** SQRTN >= sqrt(N)
> *** These could be done in code, of course
> *** Total memory usage (N + P) * sizeof(int)
> ***/
>
> #define M /* Set this to whatever you want */
> #define N 65536
> #define SQRTN 256
> #define P 6600
>
> main() {
> unsigned int buffer[N], primes[P], last_prime=0;
> unsigned int i, j, k;

//From your comment at the bottom, I STARTED to suspect
// your primes[P] array is not big enough. Sure, all
// the primes that fit in 16 bits can be used to find
// the primes that fit in 32 bits, but there only
// 203million or so such primes. At four bytes each
// -- Ah, my error! printf() is going to output more
// than 4 bytes per prime, and making a 1.5GB file
// could happen easily. OK, so your program didn't
// wrap around the 32-bit mark.

> for(i = 2; i < N; i++) buffer[i] = i;
>
> for(i= 2; i < N; i++)
> if(buffer[i] && (primes[last_prime++] = i)
> && i < SQRTN)
> for(j = i * i; j < N; j += i) buffer[j] = 0;
>
> for(i = 0; i < last_prime; i++)
> printf("%d\n", primes[i]);

// I'm pretty sure you can make most parts of
// this whole program go significantly faster
// yet by using pointer arithmetic instead of
// array access. See, every time the compiler
// sees something like primes[i], it creates
// an assembly multiplication to figure the
// offset from the start of the array, to the
// particular element. Well, if you are
// walking the array, then, for example, the
// last two lines above could be wriiten as:
// unsigned int *p; // extra declaration needed!
// for(i=0, p=primes; i<last_prime; i++)
// printf("%u\n", *p++); // addition-access to data
// Obviously variants of that example could
// be done for this whole program.
 
> for(k = 1; k <= M; k++) {
> for(i = 1; i < N; i++) buffer[i] = k*N + i;
>
> for(j = 0; (j < last_prime) && primes[j]; j++)
> for(i = (primes[j] - ((k*N) % primes[j]));
> i < N; i += primes[j]) {
> buffer[i] = 0;
> }

//I think there is a way to not need that modulus
//operation, but since I can see it is only used
//for initializing the loop, the time saved might
//not be worth the extra code & memory that would
//needed to replace it.
 
> for(i = 0; i < N; i++)
> if(buffer[i]) printf("%d\n", buffer[i]);
> }
>
> return 0;
> }
>
> This gets the first million primes (from 2 to 15485863)
> in just over six seconds, and should scale linearly
> (i.e. it takes about twice as long to do test twice as
> many numbers). Just set M. I need to find out how to
> get printf to work with unsigned numbers (it started
> printing negative numbers out),

That would be the %u in the example I did above. If you
would like hexadecimal output, then %X is the way to go.

> but it produced about 1.5 GB of primes in about ten
> minutes, at which point the filesystem restrictions
> on the windows machine came into play. I'll run it
> on my (much slower) linux box later, and put in those
> corrections as well.

Thanks for the info! And I am starting to be quite sure
that what happened before was that the malloc(15000) was
insufficient, and the primes you were saving started
overwriting the data in the accumulators. Bad data there
meant that the algorithm couldn't work right.

Relevant Pages

  • Re: greatest multiple algorithm
    ... array stores found primes ... modulo with found primes in array to find new primes ... consumes large amounts of memory ... That doesn't sound like a sieve, ...
    (comp.lang.asm.x86)
  • Re: [Algorithm] Sum of Primes < 1000000
    ... Note also skipping the inner loop for composite values of i. ... The number of primes is O) from what I've seen of primes, which is Osince the log of the square root is simply half the log and the half disappears into the constant factor. ... The inner loop does one imul, three iadds (one's actually the subtract in the loop test and one computes the array index), one icmp, and one boolean assignment. ...
    (comp.lang.java.programmer)
  • Re: More on triangle numbers and primes!
    ... > want to extract or the limit of the array or file you have for storage. ... > spots by identifying the location as a prime or if an ... > for listing out the primes. ... Hmm...Do you have a pattern that runs in accelerating numerical value steps ...
    (sci.math)
  • Re: Request for comments on simple Ada program
    ... >>using a packed array of boolean. ... The program takes input from the command line and calculates 3 ... It simply reports the number of primes ... PROGRAM OUTPUT ...
    (comp.lang.ada)