Prev: B+ tree
Next: Arrays and like objects
From: user923005 on 13 Jun 2008 23:43
On Jun 13, 12:03 pm, user923005 <dcor...(a)connx.com> wrote:
> On Jun 13, 10:27 am, jzakiya <jzak...(a)gmail.com> wrote:
>
>
>
>
>
> > This is to announce the release of my paper "Ultimate Prime Sieve --
> > Sieve of Zakiiya (SoZ)" in which I show and explain the development of
> > a class of Number Theory Sieves to generate prime numbers.   I used
> > Ruby 1.9.0-1 as my development environment on a P4 2.8 Ghz laptop.
>
> > You can get the pdf of my paper and Ruby and Python source from here:
>
> >http://www.4shared.com/dir/7467736/97bd7b71/sharing.html
>
> > Below is a sample of one of the simple prime generators. I did a
> > Python version of this in my paper (see Python source too).  The Ruby
> > version below is the minimum array size version, while the Python has
> > array of size N (I made no attempt to optimize its implementation,
> > it's to show the method).  See my paper for what/why is going on here.
>
> > class Integer
> >    def primesP3a
> >       # all prime candidates > 3 are of form  6*k+1 and 6*k+5
> >       # initialize sieve array with only these candidate values
> >       # where sieve contains the odd integers representatives
> >       # convert integers to array indices/vals by  i = (n-3)>>1
> > =(n>>1)-1
> >       n1, n2 = -1, 1;  lndx= (self-1) >>1;  sieve = []
> >       while n2 < lndx
> >          n1 +=3;   n2 += 3;   sieve[n1] = n1;  sieve[n2] = n2
> >       end
> >       #now initialize sieve array with (odd) primes < 6, resize array
> >       sieve[0] =0;  sieve[1]=1;  sieve=sieve[0..lndx-1]
>
> >       5.step(Math.sqrt(self).to_i, 2) do |i|
> >          next unless sieve[(i>>1) - 1]
> >          # p5= 5*i,  k = 6*i,  p7 = 7*i
> >          # p1 = (5*i-3)>>1;  p2 = (7*i-3)>>1;  k = (6*i)>>1
> >          i6 = 6*i;  p1 = (i6-i-3)>>1;  p2 = (i6+i-3)>>1;  k = i6>>1
> >          while p1 < lndx
> >              sieve[p1] = nil;  sieve[p2] = nil;  p1 += k;  p2 += k
> >          end
> >       end
> >       return [2] if self < 3
> >       [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
> >    end
> > end
>
> > def primesP3(val):
> >     # all prime candidates > 3 are of form  6*k+(1,5)
> >     # initialize sieve array with only these candidate values
> >     n1, n2 = 1, 5
> >     sieve = [False]*(val+6)
> >     while  n2 < val:
> >         n1 += 6;   n2 += 6;  sieve[n1] = n1;   sieve[n2] = n2
> >     # now load sieve with seed primes 3 < pi < 6, in this case just 5
> >     sieve[5] = 5
>
> >     for i in range( 5, int(ceil(sqrt(val))), 2) :
> >        if not sieve[i]:  continue
> >        #  p1= 5*i,  k = 6*i,  p2 = 7*i,
> >        p1 = 5*i;  k = p1+i;  p2 = k+i
> >        while p2 <= val:
> >           sieve[p1] = False;  sieve[p2] = False;  p1 += k;  p2 += k
> >        if p1 <= val:  sieve[p1] = False
>
> >     primes = [2,3]
> >     if val < 3 : return [2]
> >     primes.extend( i for i in range(5, val+(val&1), 2)  if sieve[i] )
>
> >     return primes
>
> > Now to generate an array of the primes up to some N just do:
>
> > Ruby:      10000001.primesP3a
> > Python:   primesP3a(10000001)
>
> > The paper presents benchmarks with Ruby 1.9.0-1 (YARV).  I would love
> > to see my various prime generators benchmarked with optimized
> > implementations in other languages.  I'm hoping C/C++ gurus will do
> > good implementations.  The methodology is very simple, since all I do
> > is additions, multiplications, and array reads/writes.
>
> > I would also like to the C implementations benchmarked against the
> > versions create by Daniel J Bernstein of the Sieve of Atkin (SoA). The
> > C code is here:
>
> >http://cr.yp.to/primesgen.html
>
> I might give it a go.  Bernstein's primegen is not your best
> competition.  This is:http://www.primzahlen.de/files/referent/kw/sieb.htm
>
> Here is output using a single 3 GHz CPU:
>
> Microsoft Windows [Version 6.0.6001]
> Copyright (c) 2006 Microsoft Corporation.  All rights reserved.
>
> C:\math\sieve\ecprime\x64\Release 64>ecprime 100000000000
>
> 100%
> primes: 4118054813
> time: 60.153 sec
> C:\math\sieve\ecprime\x64\Release 64>
>
> So your mark is to beat sieving through 100 billion in a minute.
>
> I might give your algorithm a crack (no promises though).  I set
> follow-ups to news:comp.programming

Now that I have read your paper, it seems to me that this is quite an
old idea.

If we use a product of some small primes, we can reduce our storage
and operations by a small, constant factor.

Consider this thread, from 1998:
http://groups.google.com/group/comp.lang.asm.x86/browse_thread/thread/9672ddba66665cb6/2219d96e15c891fa?lnk=st&q=sieve+author%3ATerje#2219d96e15c891fa

Here is a 20 year old implementation by Terje M. in C that uses the
idea with 2*3*5:

//
// Terje Mathisen wrote this beautiful sieve algorithm in Feb, 1998
//
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

unsigned char *buf;

unsigned char mask_tab[30] = {
0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0,
0, 0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0, 128
};

unsigned max_prime_number;

int test_dm(unsigned d, unsigned m)
{
return mask_tab[m] && ((buf[d] & mask_tab[m]) == 0);
}

void set_dm(unsigned d, unsigned m)
{
buf[d] |= mask_tab[m];
}

unsigned divmod30(unsigned index, unsigned *pm)
{
// d = index / 30;
unsigned d;
__int64 t = (__int64) index * 2290649225;
d = *(((unsigned *) &t) + 1);
d >>= 4;
*pm = index - d * 30;
return d;
}

unsigned next_prime(unsigned p)
{
unsigned d,
m;

if (p & 1)
p++;
p++;
d = divmod30(p, &m);
while (m < 30) {
if (p > max_prime_number)
return 0;
if (test_dm(d, m))
return p;
m += 2;
p += 2;
}
m = 1;
d++;
while (buf[d] == 0xff) {
p += 30;
d++;
if (p > max_prime_number)
return 0;
}
do {
if (test_dm(d, m))
return p;
m += 2;
p += 2;
} while (p <= max_prime_number);
return 0;
}

unsigned sieve(unsigned alloc, unsigned max)
{
unsigned maxtest = (unsigned) sqrt((double) max);
unsigned curr,
step,
dstep[30],
next_m[30],
dcurr,
mcurr;
unsigned prime,
i;
unsigned count;
unsigned zero_bits[256];

max_prime_number = max;

prime = 5;
for (;;) {
prime = next_prime(prime);
if (prime > maxtest)
break;

step = prime * 2;
curr = prime * prime;
dcurr = divmod30(curr, &mcurr);

for (i = 1; i < 30; i += 2) {
unsigned d,
m,
s = i;
do {
s += step;
d = divmod30(s, &m);
} while (mask_tab[m] == 0);
dstep[i] = d;
next_m[i] = m;
}

do {
if ((buf[dcurr] & mask_tab[mcurr]) == 0)
set_dm(dcurr, mcurr);

dcurr += dstep[mcurr];
mcurr = next_m[mcurr];
} while (dcurr < alloc);
}
// Count the primes!
printf("...Counting!\n");


for (curr = 0; curr < 256; curr++) {
unsigned c = curr;
for (count = 8; c; count--) {
c &= c - 1;
}
zero_bits[curr] = count;
}

for (count = 2, curr = 0; curr < (max / 30); curr++) {
count += zero_bits[buf[curr]];
}
for (prime = curr * 30; prime < max;) {
prime = next_prime(prime);
if (!prime)
break;
count++;
}

return count;
}

int main(int argc, char *argv[])
{
unsigned count,
alloc,
s,
max,
show;
time_t begin;
if (argc > 1)
max = (unsigned) atof(argv[1]);
else
max = 1000L;

if (argc > 2)
show = (unsigned) atol(argv[2]);
else
show = 100;

alloc = (max + 29) / 30 + 4;
buf = (unsigned char *) calloc(alloc, 1);
if (buf == NULL) {
printf("Cannot handle this many numbers!\n");
return EXIT_FAILURE;
}
printf("Searching prime numbers to : %ld", max);
begin = time(NULL);

count = sieve(alloc, max);

printf(" %ld prime numbers found in %ld secs\n",
count, time(NULL) - begin);

s = 1;
for (;;) {
if (s >= 5) {
s = next_prime(s);
if (!s || (s >= show))
break;
} else if (s >= 3)
s = 5;
else if (s >= 2)
s = 3;
else
s = 2;
printf("%8d", s);
}
return 0;
}


P.S.
Miller-Rabin shows this to be composite immediately:
921133374485363450332059416899657
How long will your algorithm take?

It is not very important to have a fast prime/note-prime test for
small numbers (e.g. 32 bit numbers). Any method will do, including
brute force. It when the numbers are hundreds of digits that it
becomes very important how you try to find the answers.
From: user923005 on 13 Jun 2008 23:46
On Jun 13, 8:43 pm, user923005 <dcor...(a)connx.com> wrote:
> On Jun 13, 12:03 pm, user923005 <dcor...(a)connx.com> wrote:
>
>
>
>
>
> > On Jun 13, 10:27 am, jzakiya <jzak...(a)gmail.com> wrote:
>
> > > This is to announce the release of my paper "Ultimate Prime Sieve --
> > > Sieve of Zakiiya (SoZ)" in which I show and explain the development of
> > > a class of Number Theory Sieves to generate prime numbers.   I used
> > > Ruby 1.9.0-1 as my development environment on a P4 2.8 Ghz laptop.
>
> > > You can get the pdf of my paper and Ruby and Python source from here:
>
> > >http://www.4shared.com/dir/7467736/97bd7b71/sharing.html
>
> > > Below is a sample of one of the simple prime generators. I did a
> > > Python version of this in my paper (see Python source too).  The Ruby
> > > version below is the minimum array size version, while the Python has
> > > array of size N (I made no attempt to optimize its implementation,
> > > it's to show the method).  See my paper for what/why is going on here.
>
> > > class Integer
> > >    def primesP3a
> > >       # all prime candidates > 3 are of form  6*k+1 and 6*k+5
> > >       # initialize sieve array with only these candidate values
> > >       # where sieve contains the odd integers representatives
> > >       # convert integers to array indices/vals by  i = (n-3)>>1
> > > =(n>>1)-1
> > >       n1, n2 = -1, 1;  lndx= (self-1) >>1;  sieve = []
> > >       while n2 < lndx
> > >          n1 +=3;   n2 += 3;   sieve[n1] = n1;  sieve[n2] = n2
> > >       end
> > >       #now initialize sieve array with (odd) primes < 6, resize array
> > >       sieve[0] =0;  sieve[1]=1;  sieve=sieve[0..lndx-1]
>
> > >       5.step(Math.sqrt(self).to_i, 2) do |i|
> > >          next unless sieve[(i>>1) - 1]
> > >          # p5= 5*i,  k = 6*i,  p7 = 7*i
> > >          # p1 = (5*i-3)>>1;  p2 = (7*i-3)>>1;  k = (6*i)>>1
> > >          i6 = 6*i;  p1 = (i6-i-3)>>1;  p2 = (i6+i-3)>>1;  k = i6>>1
> > >          while p1 < lndx
> > >              sieve[p1] = nil;  sieve[p2] = nil;  p1 += k;  p2 += k
> > >          end
> > >       end
> > >       return [2] if self < 3
> > >       [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
> > >    end
> > > end
>
> > > def primesP3(val):
> > >     # all prime candidates > 3 are of form  6*k+(1,5)
> > >     # initialize sieve array with only these candidate values
> > >     n1, n2 = 1, 5
> > >     sieve = [False]*(val+6)
> > >     while  n2 < val:
> > >         n1 += 6;   n2 += 6;  sieve[n1] = n1;   sieve[n2] = n2
> > >     # now load sieve with seed primes 3 < pi < 6, in this case just 5
> > >     sieve[5] = 5
>
> > >     for i in range( 5, int(ceil(sqrt(val))), 2) :
> > >        if not sieve[i]:  continue
> > >        #  p1= 5*i,  k = 6*i,  p2 = 7*i,
> > >        p1 = 5*i;  k = p1+i;  p2 = k+i
> > >        while p2 <= val:
> > >           sieve[p1] = False;  sieve[p2] = False;  p1 += k;  p2 += k
> > >        if p1 <= val:  sieve[p1] = False
>
> > >     primes = [2,3]
> > >     if val < 3 : return [2]
> > >     primes.extend( i for i in range(5, val+(val&1), 2)  if sieve[i] )
>
> > >     return primes
>
> > > Now to generate an array of the primes up to some N just do:
>
> > > Ruby:      10000001.primesP3a
> > > Python:   primesP3a(10000001)
>
> > > The paper presents benchmarks with Ruby 1.9.0-1 (YARV).  I would love
> > > to see my various prime generators benchmarked with optimized
> > > implementations in other languages.  I'm hoping C/C++ gurus will do
> > > good implementations.  The methodology is very simple, since all I do
> > > is additions, multiplications, and array reads/writes.
>
> > > I would also like to the C implementations benchmarked against the
> > > versions create by Daniel J Bernstein of the Sieve of Atkin (SoA). The
> > > C code is here:
>
> > >http://cr.yp.to/primesgen.html
>
> > I might give it a go.  Bernstein's primegen is not your best
> > competition.  This is:http://www.primzahlen.de/files/referent/kw/sieb.htm
>
> > Here is output using a single 3 GHz CPU:
>
> > Microsoft Windows [Version 6.0.6001]
> > Copyright (c) 2006 Microsoft Corporation.  All rights reserved.
>
> > C:\math\sieve\ecprime\x64\Release 64>ecprime 100000000000
>
> > 100%
> > primes: 4118054813
> > time: 60.153 sec
> > C:\math\sieve\ecprime\x64\Release 64>
>
> > So your mark is to beat sieving through 100 billion in a minute.
>
> > I might give your algorithm a crack (no promises though).  I set
> > follow-ups to news:comp.programming
>
> Now that I have read your paper, it seems to me that this is quite an
> old idea.
>
> If we use a product of some small primes, we can reduce our storage
> and operations by a small, constant factor.
>
> Consider this thread, from 1998:http://groups.google.com/group/comp.lang.asm.x86/browse_thread/thread...
>
> Here is a 20 year old implementation by Terje M. in C that uses the
> idea with 2*3*5:
>
> //
> // Terje Mathisen wrote this beautiful sieve algorithm in Feb, 1998
> //
> #include <stdio.h>
> #include <stdlib.h>
> #include <time.h>
> #include <math.h>
>
> unsigned char  *buf;
>
> unsigned char   mask_tab[30] = {
>     0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0,
>     0, 0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0, 128
>
> };
>
> unsigned        max_prime_number;
>
> int      test_dm(unsigned d, unsigned m)
> {
>     return mask_tab[m] && ((buf[d] & mask_tab[m]) == 0);
>
> }
>
> void     set_dm(unsigned d, unsigned m)
> {
>     buf[d] |= mask_tab[m];
>
> }
>
> unsigned divmod30(unsigned index, unsigned *pm)
> {
> //  d = index / 30;
>     unsigned        d;
>     __int64         t = (__int64) index * 2290649225;
>     d = *(((unsigned *) &t) + 1);
>     d >>= 4;
>     *pm = index - d * 30;
>     return d;
>
> }
>
> unsigned next_prime(unsigned p)
> {
>     unsigned        d,
>                     m;
>
>     if (p & 1)
>         p++;
>     p++;
>     d = divmod30(p, &m);
>     while (m < 30) {
>         if (p > max_prime_number)
>             return 0;
>         if (test_dm(d, m))
>             return p;
>         m += 2;
>         p += 2;
>     }
>     m = 1;
>     d++;
>     while (buf[d] == 0xff) {
>         p += 30;
>         d++;
>         if (p > max_prime_number)
>             return 0;
>     }
>     do {
>         if (test_dm(d, m))
>             return p;
>         m += 2;
>         p += 2;
>     } while (p <= max_prime_number);
>     return 0;
>
> }
>
> unsigned        sieve(unsigned alloc, unsigned max)
> {
>     unsigned        maxtest = (unsigned) sqrt((double) max);
>     unsigned        curr,
>                     step,
>                     dstep[30],
>                     next_m[30],
>                     dcurr,
>                     mcurr;
>     unsigned        prime,
>                     i;
>     unsigned        count;
>     unsigned        zero_bits[256];
>
>     max_prime_number = max;
>
>     prime = 5;
>     for (;;) {
>         prime = next_prime(prime);
>         if (prime > maxtest)
>             break;
>
>         step = prime * 2;
>         curr = prime * prime;
>         dcurr = divmod30(curr, &mcurr);
>
>         for (i = 1; i < 30; i += 2) {
>             unsigned        d,
>                             m,
>                             s = i;
>             do {
>                 s += step;
>                 d = divmod30(s, &m);
>             } while (mask_tab[m] == 0);
>             dstep[i] = d;
>             next_m[i] = m;
>         }
>
>         do {
>             if ((buf[dcurr] & mask_tab[mcurr]) == 0)
>                 set_dm(dcurr, mcurr);
>
>             dcurr += dstep[mcurr];
>             mcurr = next_m[mcurr];
>         } while (dcurr < alloc);
>     }
>     // Count the primes!
>     printf("...Counting!\n");
>
>     for (curr = 0; curr < 256; curr++) {
>         unsigned        c = curr;
>         for (count = 8; c; count--) {
>             c &= c - 1;
>         }
>         zero_bits[curr] = count;
>     }
>
>     for (count = 2, curr = 0; curr < (max / 30); curr++) {
>         count += zero_bits[buf[curr]];
>     }
>     for (prime = curr * 30; prime < max;) {
>         prime = next_prime(prime);
>         if (!prime)
>             break;
>         count++;
>     }
>
>     return count;
>
> }
>
> int             main(int argc, char *argv[])
> {
>     unsigned        count,
>                     alloc,
>                     s,
>                     max,
>                     show;
>     time_t          begin;
>     if (argc > 1)
>         max = (unsigned) atof(argv[1]);
>     else
>         max = 1000L;
>
>     if (argc > 2)
>         show = (unsigned) atol(argv[2]);
>     else
>         show = 100;
>
>     alloc = (max + 29) / 30 + 4;
>     buf = (unsigned char *) calloc(alloc, 1);
>     if (buf == NULL) {
>         printf("Cannot handle this many numbers!\n");
>         return EXIT_FAILURE;
>     }
>     printf("Searching prime numbers to : %ld", max);
>     begin = time(NULL);
>
>     count = sieve(alloc, max);
>
>     printf(" %ld prime numbers found in %ld secs\n",
>            count, time(NULL) - begin);
>
>     s = 1;
>     for (;;) {
>         if (s >= 5) {
>             s = next_prime(s);
>             if (!s || (s >= show))
>                 break;
>         } else if (s >= 3)
>             s = 5;
>         else if (s >= 2)
>             s = 3;
>         else
>             s = 2;
>         printf("%8d", s);
>     }
>     return 0;
>
> }
>
> P.S.
> Miller-Rabin shows this to be composite immediately:
> 921133374485363450332059416899657

Here is Pari/GP output:
(20:41) gp > c = isprime(921133374485363450332059416899657);
(20:45) gp > print (c);
0
(20:45) gp > c = isprime(7);
(20:45) gp > print(c)
1
(20:45) gp >

> How long will your algorithm take?
>
> It is not very important to have a fast prime/note-prime test for
> small numbers (e.g. 32 bit numbers).  Any method will do, including
> brute force.  It when the numbers are hundreds of digits that it
> becomes very important how you try to find the answers.- Hide quoted text -
>
> - Show quoted text -

From: Juha Nieminen on 14 Jun 2008 11:50
user923005 wrote:
> Here is Pari/GP output:
> (20:41) gp > c = isprime(921133374485363450332059416899657);
> (20:45) gp > print (c);
> 0

How come Pari/GP is so slow there? Here it's enormously faster.
isprime() gives the answer immediately. Even factoring that number is
incredibly fast:

? factor(921133374485363450332059416899657)
time = 56 ms.
%1 =
[2024035091655457 1]

[455097531798210601 1]
From: user923005 on 17 Jun 2008 20:28
On Jun 14, 8:50 am, Juha Nieminen <nos...(a)thanks.invalid> wrote:
> user923005 wrote:
> > Here is Pari/GP output:
> > (20:41) gp > c = isprime(921133374485363450332059416899657);
> > (20:45) gp > print (c);
> > 0
>
>   How come Pari/GP is so slow there? Here it's enormously faster.
> isprime() gives the answer immediately. Even factoring that number is
> incredibly fast:
>
> ? factor(921133374485363450332059416899657)
> time = 56 ms.
> %1 =
> [2024035091655457 1]
>
> [455097531798210601 1]

Probably I had a lot of things in memory and so page faulting slowed
things down.
The method in the published paper would take weeks just to find out if
the above number was prime or not {operation count would be
proportional to 4,597,105,011,348,277}.
 | 
Pages: 1
Prev: B+ tree
Next: Arrays and like objects