From: David Kuestler <kuestler@zeta.org.au>
Subject: Re: Prime Numbers Generator
Date: 28 Apr 1999 00:00:00 GMT
Message-ID: <37270535.78A83EC6@zeta.org.au>
References: <371FC9D4.16853337@aspi.net> <3721B8D4.CC662B6E@zeta.org.au> <37227EFB.ED3B0002@acm.org> <3722CD3A.D22C9A2@zeta.org.au> <37232CB5.F29DAC09@acm.org> <3723FC96.95F8C6C0@zeta.org.au> <3725CAE7.7F2F84E9@NOSPAMhotmail.com>
X-Accept-Language: en
Content-Type: multipart/mixed; boundary="------------703CCF9197824E830D9AD423"
X-Complaints-To: abuse@zeta.org.au
X-Trace: phaedrus.zeta.org.au 925304111 29900 203.26.9.97 (28 Apr 1999 12:55:11 GMT)
Organization: Zeta Internet, http://www.zeta.org.au/
Mime-Version: 1.0
Reply-To: kuestler@NOJUNKMAILzeta.org.au
NNTP-Posting-Date: 28 Apr 1999 12:55:11 GMT
Newsgroups: sci.crypt

This is a multi-part message in MIME format.
--------------703CCF9197824E830D9AD423
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Peter Gunn wrote:

> David Kuestler wrote:
>
> > Jim Gillogly wrote
>
> > [snip]
>
> > Excellent ! This ran just shy of an hour on a 233Mhz K6 and gave the
> > same file size, though the byte order is different ( as expected ).
> > You're certainly a better C programmer than I as you only use integers
> > and pointers where as I used array references with the occasional
> > floating point operation.
> >
> > If you are interested you should submit it to  The Prime Pages . Your
> > code is faster, smaller, more understandable and therefore would be
> > easier to impliment than others.
>
> Have a look at http://www.smdp.freeserve.co.uk/erato.cpp
> for a somewhat faster version... should run in ~10mins or

11.5 minutes ! ( with the file output back in ).

>
> less on a 233Mhz PC.  Basically it uses bitmaps instead
> of character arrays and an 8Mb buffer to speed thing up.
>
> Im sure this could be reduced to less than 2mins runtime
> & 1Mb memory usage if it was hacked a bit further.
>
> I suppose you could even use a 2D graphics accelerator
> to make things *much* faster? That would be interesting...
> might even be worth watching :-)
>
> tata
>
> PG.
>
> PS I took out the code that output the primes to a file :-)

I put them back in to to give me a 'better' timing comparison to other
routines I've run.
I've got a lot to learn about C coding/tuning :-{
I will try Jim's blocking technique and your bitmaps on my 64 bit prime
number field routines.

Thanks to both you and Jim :-)


--------------703CCF9197824E830D9AD423
Content-Type: text/plain; charset=us-ascii;
 name="erato.cpp"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="erato.cpp"

/*
* erato: sieve of eratosthenes for up to 2**32.
* Two steps: first find all primes less than 2**16, then use them to sift
* the remaining numbers in blocks.
*
* Jim Gillogly, 24 Apr 1999
*
* Hacked by Peter Gunn 26 Apr 1999
*
* NOTE: Switch off Gloabal Optimsation & Full Optimisation when using
*       MS VC++ 5, coz its broken!!!
*/

#include <stdio.h>
#include <memory.h>
#include <time.h>

#define MAXSMALL        6542            /* 6542 for primes up to 1<<16. */
#define TOP             (1 << 16)       /* Top of each block to be sifted. */
#define PRFILE          "primes0"       /* Save 'em in a file. */
#define BLOCKS          (1 << 16)       /* How many blocks to do total? */

unsigned long sieve[TOP/32];            /* Turned off if composite. */
#define bit(x) (sieve[(x)>>5]&(1<<(x&31)))
#define cbit(x) sieve[(x)>>5]&=~(1<<(x&31))

/* bitmasks for all primes < 32 in all rotations */
unsigned long sieve2[32][32][TOP/32];
#define csbit(y,r,x) sieve2[y][r][(x)>>5]&=~(1<<(x&31))

//char *  stop = &sieve[TOP];

unsigned long   smalls[MAXSMALL];       /* All primes under 2**16. */
unsigned long   smallnum = 0;

main()
{
   unsigned long   b,b1, i, j, k, nprimes, start, step, *s;
   char            *p, *q, xxx[40];
   FILE            *out;
   unsigned long start_time,stop_time;

   start_time=time((time_t *)NULL);
   memset(sieve,0xff,sizeof(sieve));
   cbit(0); cbit(1);
   b=0;
   do
   {
    b++;
    while (b<TOP && !bit(b)) b++;  /* Scan to next prime. */
    b1=b;
    while ((b1+=b)<TOP)
     cbit(b1);
   } while (b<TOP);

   if ((out = fopen(PRFILE, "w")) == NULL)
   {
    fprintf(stderr, "Phui.\n");
    return 1;
   }

   for (i=0, s=smalls; i<TOP; i++)
   {
    if (bit(i))
    {
     *s=i;
          fwrite(&i, 4, 1, out);

     if (i<32)
     {
      unsigned long rot;
      for (rot=0;rot<32;rot++)
      {
       memset(&sieve2[s-smalls][rot],0xff,sizeof(sieve2[0][0]));
       b1=rot%i;
       csbit(s-smalls,rot,b1);
       while ((b1+=i)<TOP)
        csbit(s-smalls,rot,b1);
      }
     }
     s++;
    }
   }
    smallnum = nprimes = (s - smalls);
    printf("%d small primes.\n", smallnum);
    /* Now sift each 2**16-number block with these small primes. */
    start=0;
    for (k=1; k<BLOCKS; k++)
    {
    start=k*TOP;
       memset(sieve,0xff,sizeof(sieve));
    cbit(0);
    for (i=0, s=smalls; i<smallnum; i++) /* small primes */
    {
     if ((step=*s++)<32)
     {
      unsigned long rot=step-(start%step);
      unsigned long *s1=&sieve[TOP/32-1],
        *s2=&sieve2[s-smalls-1][rot][TOP/32-1];
      while (s1>=sieve)
       *s1--&=*s2--;
     }
     else
     {
      for (b=step-(start%step);b<TOP;b+=step)
       cbit(b);
     }
    }
    for (i=0; i<TOP; i++)
    {
     if (bit(i))
     {
      j=i+start;
//    printf("%d\n",j);
      nprimes++;
//    if (!(nprimes&0xfffff))
//     printf("%ld primes in %ld secs (k=%ld)\n",nprimes,
//            time((time_t *)NULL)-start_time,k);
      fwrite(&j, 4, 1, out);
     }
    }
    }
   stop_time=time((time_t *)NULL);
   printf("%d primes in %ld secs \n", nprimes,stop_time-start_time);
   printf("hit enter\n");
   fclose(out);
    return 0;
}

--------------703CCF9197824E830D9AD423--