DZone Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world

Snippets has posted 5883 posts at DZone. View Full User Profile

Sieve Of Zakiya -- Number Theory Sieves (NTS)

06.06.2008
| 3770 views |
  • submit to reddit
        // description of your code here
I present with this Ruby code a new class of Number Theory Sieves (NTS)
to generate primes numbers upto some number N, which I believe now
represent the fastest way to generate prime numbers. The general number
theory can utilize a class of prime generator functions, which have the
same form, and which can be implemented in the same manner. These NTS
can be inherently implemented in parallel with multiple processors.

The classical brute-force method to generate primes upto some N is the
Sieve of Eratothenes (SoE). In 2002/3 Atkin & Bernstein proposed what has
become known as the Sieve of Atkin (SoA), which was known to be the fastest
method to generate primes. My general number theory based Sieve of Zakiya (SoZ)
represents a multiple times increase in performance over the SoA, while being
easier to understand, code, and extend with newer prime generator functions.

I also used the number theory to created a primality tester, whose code
is short and elegant, which I also believe is the fastest method to test
for primality while being deterministic. Thus, it should replace other
non-deterministic or probabilistic methods, such as Miller-Rabin, et al.

I started this work this first week in May 2008, and present this code
herein on Friday, June 6, 2008.

The development of this work is presented in the paper:
"Ultimate Prime Sieve -- Sieve of Zakiya"
the pdf can be download from here:

http://www.4shared.com/dir/7467736/97bd7b71/sharing.html

#!/usr/local/bin/ruby -w

require 'benchmark' 

class Integer

   def primes1
      # classical brute-force Sieve of Eratosthenes (SoE)
      # initialze sieve array with all integers starting at i=2

      sieve = [nil, nil].concat((2..self).to_a)

      (2 .. Math.sqrt(self).to_i).each do |i|

         next unless sieve[i]
         (i*i).step(self, i) { |j| sieve[j] = nil }
      end
      sieve.compact!
   end  

   def primes2
      # classical brute-force Sieve of Eratosthenes (SoE)
      # initialize sieve array with odd integers for odd indexes, then 2

      sieve = [];  3.step(self, 2) { |i| sieve[i] = i };  sieve[2] = 2

      3.step(Math.sqrt(self).to_i, 2) do |i| 
         next unless sieve[i]
         (i*i).step(self, i<<1) { |j| sieve[j] = nil }
      end
      sieve.compact!
   end 

   def primes3
      # http://everything2.com/index.pl?node_id=1176369
      # all prime candidates > 3 are of form  6*k+1 and 6*k+5
      # initialize sieve array with only these candidate values
      n1, n2 = 1, 5;  sieve = [] 
      while n2 < self
         n1 += 6;   n2 += 6;   sieve[n1] = n1;   sieve[n2] = n2
      end
      # now initialize sieve array with first 3 primes, resize array
      sieve[2]=2;  sieve[3]=3;  sieve[5]=5;  sieve=sieve[0..self]

      5.step(Math.sqrt(self).to_i, 2) do |i| 
         next unless sieve[i]
         (i*i).step(self, i<<1) {|j| sieve[j] = nil }
      end   	
      sieve.compact!
   end

   def primesP3
      # all prime candidates > 3 are of form  6*k+1 and 6*k+5
      # initialize sieve array with only these candidate values
      n1, n2 = 1, 5;  sieve = [] 
      while n2 < self
         n1 += 6;   n2 += 6;   sieve[n1] = n1;  sieve[n2] = n2
      end
      # now initialize sieve array with primes < 6, resize array
      sieve[2]=2;  sieve[3]=3;  sieve[5]=5;  sieve=sieve[0..self]

      5.step(Math.sqrt(self).to_i, 2) do |i| 
         next unless sieve[i]
	 #  p5= 5*i, k = 6*i,  p7 = 7*i  
	 p1 = 5*i;  k = p1+i;  p2 = k+i
	 while p1 <= self
	     sieve[p1] = nil;  sieve[p2] = nil;  p1 += k;  p2 += k
	 end    
      end   	
      sieve.compact!
   end


   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

   def primesP5
      # all prime candidates > 5 are of form  30*k+(1,7,11,13,17,19,23,29)
      # initialize sieve array with only these candidate values
      n1, n2, n3, n4, n5, n6, n7, n8 = 1, 7, 11, 13, 17, 19, 23, 29;  sieve = [] 
      while n8 < self
         n1 +=30;   n2 += 30;   n3 += 30;   n4 += 30
	 n5 +=30;   n6 += 30;   n7 += 30;   n8 += 30
	 sieve[n1] = n1;  sieve[n2] = n2;  sieve[n3] = n3;  sieve[n4] = n4
	 sieve[n5] = n5;  sieve[n6] = n6;  sieve[n7] = n7;  sieve[n8] = n8
      end
      # now initialize sieve with the primes < 30, resize array
      sieve[2]=2;  sieve[3]=3;  sieve[5]=5;  sieve[7]=7;  sieve[11]=11;  sieve[13]=13
      sieve[17]=17; sieve[19]=19;  sieve[23]=23;  sieve[29]=29;   sieve=sieve[0..self]

      7.step(Math.sqrt(self).to_i, 2) do |i| 
         next unless sieve[i]
	 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i,  k=30*i
	 n6 = 6*i;  p1 = i+n6;  p2 = n6+n6-i;  p3 = p1+n6;  p4 = p2+n6
	 p5 = p3+n6;  p6 = p4+n6;  p7 = p6+n6;  k = p7+i;  p8 = k+i
	 while p1 <= self
	     sieve[p1] = nil;  sieve[p2] = nil;  sieve[p3] = nil;  sieve[p4] = nil
	     sieve[p5] = nil;  sieve[p6] = nil;  sieve[p7] = nil;  sieve[p8] = nil
	     p1 += k; p2 += k; p3 += k; p4 += k; p5 += k; p6 += k; p7 += k; p8 += k
	 end
      end
      sieve.compact!
   end

   def primesP5a
      # all prime candidates > 5 are of form  30*k+(1,7,11,13,17,19,23,29)
      # 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, n3, n4, n5, n6, n7, n8 = -1, 2, 4, 5, 7, 8, 10, 13
      lndx= (self-1) >>1;  sieve = []
      while n8 < lndx
         n1 +=15;   n2 += 15;   n3 += 15;   n4 += 15
	 n5 +=15;   n6 += 15;   n7 += 15;   n8 += 15
	 sieve[n1] = n1;  sieve[n2] = n2;  sieve[n3] = n3;  sieve[n4] = n4
	 sieve[n5] = n5;  sieve[n6] = n6;  sieve[n7] = n7;  sieve[n8] = n8
      end
      # now initialize sieve with the (odd) primes < 30, resize array
      sieve[0]=0;  sieve[1]=1;  sieve[2]=2;  sieve[4]=4;  sieve[5]=5
      sieve[7]=7;  sieve[8]=8;  sieve[10]=10;  sieve[13]=13
      sieve = sieve[0..lndx-1]

      7.step(Math.sqrt(self).to_i, 2) do |i|
         next unless sieve[(i>>1) - 1]
	 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i,  k=30*i
	 # maps to: p1= (7*i-3)>>1, --- p8=(31*i-3)>>1, k=30*i>>1
	 n2=i<<1;  n4=i<<2;  n8=i<<3;  n16=i<<4;  n32=i<<5;  n12=n8+n4
	 p1 = (n8-i-3)>>1;  p2 = (n12-i-3)>>1;  p3 = (n12+i-3)>>1;  p4 = (n16+i-3)>>1
	 p5 = (n16+n4-i-3)>>1;  p6 = (n16+n8-i-3)>>1;  p7 = (n32-n4+i-3)>>1
	 p8 = (n32-i-3)>>1;  k = (n32-n2)>>1
	 while p1 < lndx
	     sieve[p1] = nil;  sieve[p2] = nil;  sieve[p3] = nil;  sieve[p4] = nil
	     sieve[p5] = nil;  sieve[p6] = nil;  sieve[p7] = nil;  sieve[p8] = nil
	     p1 += k;  p2 += k;  p3 += k;  p4 += k;  p5 += k;  p6 += k;  p7 += k;  p8 += k
	 end
      end
      return [2] if self < 3
      [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
   end

   def primesP7
      # all prime candidates > 7 are of form  210*k+(1,11,13,17,19,23,29,31,37
      # 41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121,127,131
      # 137,139,143,149,151,157,163,167,169,173,179,181,187,191,193,197,199,209)
      # initialize sieve array with only these candidate values
      n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 = 1, 11, 13, 17, 19, 23, 29, 31, 37, 41
      n11, n12, n13, n14, n15, n16, n17, n18 = 43, 47, 53, 59, 61, 67, 71, 73
      n19, n20, n21, n22, n23, n24, n25, n26 = 79, 83, 89, 97, 101, 103, 107, 109
      n27, n28, n29, n30, n31, n32, n33, n34 = 113, 127, 131, 137, 139, 149, 151, 157
      n35, n36, n37, n38, n39, n40, n41, n42, n43 = 163, 167, 173, 179, 181, 191, 193, 197, 199
      n44, n45, n46, n47, n48 = 121, 143, 169, 187,  209
      num = self - self[0]^1;  sieve = []
      while n48 < num
         n1  +=210;  n2  += 210;  n3  += 210;  n4  += 210;  n5  += 210
	 n6  +=210;  n7  += 210;  n8  += 210;  n9  += 210;  n10 += 210
         n11 +=210;  n12 += 210;  n13 += 210;  n14 += 210;  n15 += 210
	 n16 +=210;  n17 += 210;  n18 += 210;  n19 += 210;  n20 += 210
         n21 +=210;  n22 += 210;  n23 += 210;  n24 += 210;  n25 += 210
	 n26 +=210;  n27 += 210;  n28 += 210;  n29 += 210;  n30 += 210
         n31 +=210;  n32 += 210;  n33 += 210;  n34 += 210;  n35 += 210
	 n36 +=210;  n37 += 210;  n38 += 210;  n39 += 210;  n40 += 210
	 n41 +=210;  n42 += 210; n43 += 210; n44 += 210; n45 += 210 
	 n46 += 210; n47 += 210; n48 += 210
	 sieve[n1] = n1;    sieve[n2]  = n2;   sieve[n3] = n3;    sieve[n4] = n4
	 sieve[n5] = n5;    sieve[n6]  = n6;   sieve[n7] = n7;    sieve[n8] = n8
	 sieve[n9] = n9;    sieve[n10] = n10;  sieve[n11] = n11;  sieve[n12] = n12
	 sieve[n13] = n13;  sieve[n14] = n14;  sieve[n15] = n15;  sieve[n16] = n16
	 sieve[n17] = n17;  sieve[n18] = n18;  sieve[n19] = n19;  sieve[n20] = n20
	 sieve[n21] = n21;  sieve[n22] = n22;  sieve[n23] = n23;  sieve[n24] = n24
	 sieve[n25] = n25;  sieve[n26] = n26;  sieve[n27] = n27;  sieve[n28] = n28
	 sieve[n29] = n29;  sieve[n30] = n30;  sieve[n31] = n31;  sieve[n32] = n32
	 sieve[n33] = n33;  sieve[n34] = n34;  sieve[n35] = n35;  sieve[n36] = n36
	 sieve[n37] = n37;  sieve[n38] = n38;  sieve[n39] = n39;  sieve[n40] = n40
	 sieve[n41] = n41;  sieve[n42] = n42;  sieve[n43] = n43;  sieve[n44] = n44
         sieve[n45] = n44;  sieve[n46] = n46;  sieve[n47] = n47;  sieve[n48] = n48
      end
      # now initialize sieve with the (odd) primes < 210, resize array
      sieve[2]=2;  sieve[3]=3;  sieve[5]=5;  sieve[7]=7;  sieve[11]=11;  sieve[13]=13
      sieve[17]=17;   sieve[19]=19;   sieve[23]=23;   sieve[29]=29;   sieve[31]=31
      sieve[37]=37;   sieve[41]=41;   sieve[43]=43;   sieve[47]=47;   sieve[53]=53
      sieve[59]=59;   sieve[61]=61;   sieve[67]=67;   sieve[71]=71;   sieve[73]=73
      sieve[79]=79;   sieve[83]=83;   sieve[89]=89;   sieve[97]=97;   sieve[101]=101
      sieve[103]=103; sieve[107]=107; sieve[109]=109; sieve[113]=113; sieve[127]=127
      sieve[131]=131; sieve[137]=137; sieve[139]=139; sieve[149]=149; sieve[151]=151
      sieve[157]=157; sieve[163]=163; sieve[167]=167; sieve[173]=173; sieve[179]=179
      sieve[181]=181; sieve[191]=191; sieve[193]=193; sieve[197]=197; sieve[199]=199
      sieve = sieve[0..self]
      11.step(Math.sqrt(self).to_i, 2) do |i|
         next unless sieve[i]
	 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i,  k=30*i
	 # maps to: p1= (7*i-3)>>1, --- p8=(31*i)>>1, k=30*i>>1
	 #n6=6*i;  n7=n6+i;  n11=n6+n6-i;  n13=n6+n7;  n17=n11+n6
	 #n19=n13+n6;  n23=n17+n6;  n29=n23+n6;  n31=n29+(i<<1)
	 #n4=i<<2;  n8=i<<3;  n16=i<<4;  n7=n8-i;  n11=n7+n4;  n13=n11+(i<<1)
	 #n17=n16+i;  n19=n11+n8;  n23=n16+n7;  n29=n16+n13;  n31=(i<<5)-i
	 p1 = (11*i);   p2 = (13*i);   p3 = (17*i);   p4 = (19*i)  
	 p5 = (23*i);   p6 = (29*i);   p7 = (31*i);   p8 = (37*i) 
	 p9 = (41*i);   p10 = (43*i);  p11 = (47*i);  p12 = (53*i) 
	 p13 = (59*i);  p14 = (61*i);  p15 = (67*i);  p16 = (71*i) 
	 p17 = (73*i);  p18 = (79*i);  p19 = (83*i);  p20 = (89*i) 
         p21 = (97*i);  p22 = (101*i); p23 = (103*i); p24 = (107*i) 
	 p25 = (109*i); p26 = (113*i); p27 = (127*i); p28 = (131*i) 
	 p29 = (137*i); p30 = (139*i); p31 = (149*i); p32 = (151*i) 
	 p33 = (157*i); p34 = (163*i); p35 = (167*i); p36 = (173*i) 
	 p37 = (179*i); p38 = (181*i); p39 = (191*i); p40 = (193*i) 
	 p41 = (197*i); p42 = (199*i); p43 = (211*i)
	 p44 = (121*i); p45 = 143*i;  p46 = 169*i; p47 = 187*i; p48 = 209*i
	 k = (210*i) 
	 while p1 <= num
	     sieve[p1] = nil;   sieve[p2] = nil;   sieve[p3] = nil;   sieve[p4] = nil
	     sieve[p5] = nil;   sieve[p6] = nil;   sieve[p7] = nil;   sieve[p8] = nil
	     sieve[p9] = nil;   sieve[p10] = nil;  sieve[p11] = nil;  sieve[p12] = nil
	     sieve[p13] = nil;  sieve[p14] = nil;  sieve[p15] = nil;  sieve[p16] = nil
	     sieve[p17] = nil;  sieve[p18] = nil;  sieve[p19] = nil;  sieve[p20] = nil
	     sieve[p21] = nil;  sieve[p22] = nil;  sieve[p23] = nil;  sieve[p24] = nil
	     sieve[p25] = nil;  sieve[p26] = nil;  sieve[p27] = nil;  sieve[p28] = nil
	     sieve[p29] = nil;  sieve[p30] = nil;  sieve[p31] = nil;  sieve[p32] = nil
	     sieve[p33] = nil;  sieve[p34] = nil;  sieve[p35] = nil;  sieve[p36] = nil
	     sieve[p37] = nil;  sieve[p38] = nil;  sieve[p39] = nil;  sieve[p40] = nil
	     sieve[p41] = nil;  sieve[p42] = nil;  sieve[p43] = nil;  sieve[p44] = nil
             sieve[p45] = nil;  sieve[p46] = nil;  sieve[p47] = nil;  sieve[p48] = nil
	     p1 += k;  p2 += k; p3 += k; p4  += k; p5  +=k; p6 += k; p7 += k
	     p8 += k;  p9 += k; p10 +=k; p11 += k; p12 +=k; p13 +=k; p14 +=k
	     p15 +=k;  p16 +=k; p17 +=k; p18 += k; p19 +=k; p20 +=k; p21 +=k
	     p22 +=k;  p23 +=k; p24 +=k; p25 += k; p26 +=k; p27 +=k; p28 +=k
	     p29 +=k;  p30 +=k; p31 +=k; p32 += k; p33 +=k; p34 +=k; p35 +=k
	     p36 +=k;  p37 +=k; p38 +=k; p39 += k; p40 +=k; p41 +=k; p42 +=k
             p43 +=k;  p44 +=k; p45 +=k; p46 +=k;  p47 +=k; p48 +=k
	 end
      end
      sieve.compact! 
   end

   def primesP7a
      # all prime candidates > 7 are of form  210*k+(1,11,13,17,19,23,29,31,37
      # 41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121,127,131
      # 137,139,143,149,151,157,163,167,169173,179,181,187,191,193,197,199,209)
      # initialize sieve array with only these candidate values
      # where sieve contains the odd integers representatives
      # convert integers to array indices/vals by  i = (n )>>1 = (n>>1)-1
      n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 = -1, 4, 5, 7, 8, 10, 13, 14, 17, 19
      n11, n12, n13, n14, n15, n16, n17, n18 = 20, 22, 25, 28, 29, 32, 34, 35
      n19, n20, n21, n22, n23, n24, n25, n26 = 38, 40, 43, 47, 49, 50, 52, 53
      n27, n28, n29, n30, n31, n32, n33, n34 = 55, 62, 64, 67, 68, 73, 74, 77
      n35, n36, n37, n38, n39, n40, n41, n42, n43 = 80, 82, 85, 88, 89, 94, 95, 97, 98
      n44, n45, n46, n47, n48 = 59, 70, 83, 92, 103
      lndx= (self-1)>>1;  sieve = []
      while n48 < lndx
         n1 +=105;   n2 += 105;   n3 += 105;   n4 += 105;   n5 += 105
	 n6 +=105;   n7 += 105;   n8 += 105;   n9 += 105;   n10 += 105
         n11 +=105;  n12 += 105;  n13 += 105;  n14 += 105;  n15 += 105
	 n16 +=105;  n17 += 105;  n18 += 105;  n19 += 105;  n20 += 105
         n21 +=105;  n22 += 105;  n23 += 105;  n24 += 105;  n25 += 105
	 n26 +=105;  n27 += 105;  n28 += 105;  n29 += 105;  n30 += 105
         n31 +=105;  n32 += 105;  n33 += 105;  n34 += 105;  n35 += 105
	 n36 +=105;  n37 += 105;  n38 += 105;  n39 += 105;  n40 += 105
	 n41 +=105;  n42 += 105;  n43 += 105;  n44 += 105;  n45 += 105
	 n46 +=105;  n47 += 105;  n48 += 105
	 sieve[n1] = n1;    sieve[n2] = n2;    sieve[n3] = n3;    sieve[n4] = n4
	 sieve[n5] = n5;    sieve[n6] = n6;    sieve[n7] = n7;    sieve[n8] = n8
	 sieve[n9] = n9;    sieve[n10] = n10;  sieve[n11] = n11;  sieve[n12] = n12
	 sieve[n13] = n13;  sieve[n14] = n14;  sieve[n15] = n15;  sieve[n16] = n16
	 sieve[n17] = n17;  sieve[n18] = n18;  sieve[n19] = n19;  sieve[n20] = n20
	 sieve[n21] = n21;  sieve[n22] = n22;  sieve[n23] = n23;  sieve[n24] = n24
	 sieve[n25] = n25;  sieve[n26] = n26;  sieve[n27] = n27;  sieve[n28] = n28
	 sieve[n29] = n29;  sieve[n30] = n30;  sieve[n31] = n31;  sieve[n32] = n32
	 sieve[n33] = n33;  sieve[n34] = n34;  sieve[n35] = n35;  sieve[n36] = n36
	 sieve[n37] = n37;  sieve[n38] = n38;  sieve[n39] = n39;  sieve[n40] = n40
	 sieve[n41] = n41;  sieve[n42] = n42;  sieve[n43] = n43;  sieve[n44] = n44 
	 sieve[n45] = n45;  sieve[n46] = n46;  sieve[n47] = n47;  sieve[n48] = n48
      end
      # now initialize sieve with the (odd) primes < 210, resize array
      sieve[0]=0;    sieve[1]=1;    sieve[2]=2;    sieve[4]=4;    sieve[5]=5
      sieve[7]=7;    sieve[8]=8;    sieve[10]=10;  sieve[13]=13;  sieve[14]=14
      sieve[17]=17;  sieve[19]=19;  sieve[20]=20;  sieve[22]=22;  sieve[25]=25
      sieve[28]=28;  sieve[29]=29;  sieve[32]=32;  sieve[34]=34;  sieve[35]=35
      sieve[38]=38;  sieve[40]=40;  sieve[43]=43;  sieve[47]=47;  sieve[49]=49
      sieve[50]=50;  sieve[52]=52;  sieve[53]=53;  sieve[55]=55;  sieve[62]=62
      sieve[64]=64;  sieve[67]=67;  sieve[68]=68;  sieve[73]=73;  sieve[74]=74
      sieve[77]=77;  sieve[80]=80;  sieve[82]=82;  sieve[85]=85;  sieve[88]=88
      sieve[89]=89;  sieve[94]=94;  sieve[95]=95;  sieve[97]=97;  sieve[98]=98; 
      sieve = sieve[0..lndx-1]

      11.step(Math.sqrt(self).to_i, 2) do |i|
         next unless sieve[(i>>1) - 1]
	 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i,  k=30*i
	 # maps to: p1= (7*i-3)>>1, --- p8=(31*i)>>1, k=30*i>>1
	 #n6=6*i;  n7=n6+i;  n11=n6+n6-i;  n13=n6+n7;  n17=n11+n6
	 #n19=n13+n6;  n23=n17+n6;  n29=n23+n6;  n31=n29+(i<<1)
	 #n4=i<<2;  n8=i<<3;  n16=i<<4;  n7=n8-i;  n11=n7+n4;  n13=n11+(i<<1)
	 #n17=n16+i;  n19=n11+n8;  n23=n16+n7;  n29=n16+n13;  n31=(i<<5)-i
	 p1 = (11*i-3)>>1;   p2 = (13*i-3)>>1;   p3 = (17*i-3)>>1;   p4 = (19*i-3)>>1
	 p5 = (23*i-3)>>1;   p6 = (29*i-3)>>1;   p7 = (31*i-3)>>1;   p8 = (37*i-3)>>1
	 p9 = (41*i-3)>>1;   p10 = (43*i-3)>>1;  p11 = (47*i-3)>>1;  p12 = (53*i-3)>>1
	 p13 = (59*i-3)>>1;  p14 = (61*i-3)>>1;  p15 = (67*i-3)>>1;  p16 = (71*i-3)>>1
	 p17 = (73*i-3)>>1;  p18 = (79*i-3)>>1;  p19 = (83*i-3)>>1;  p20 = (89*i-3)>>1
         p21 = (97*i-3)>>1;  p22 = (101*i-3)>>1; p23 = (103*i-3)>>1; p24 = (107*i-3)>>1
	 p25 = (109*i-3)>>1; p26 = (113*i-3)>>1; p27 = (127*i-3)>>1; p28 = (131*i-3)>>1
	 p29 = (137*i-3)>>1; p30 = (139*i-3)>>1; p31 = (149*i-3)>>1; p32 = (151*i-3)>>1
	 p33 = (157*i-3)>>1; p34 = (163*i-3)>>1; p35 = (167*i-3)>>1; p36 = (173*i-3)>>1
	 p37 = (179*i-3)>>1; p38 = (181*i-3)>>1; p39 = (191*i-3)>>1; p40 = (193*i-3)>>1
	 p41 = (197*i-3)>>1; p42 = (199*i-3)>>1; p43 = (211*i-3)>>1; p44 = (121*i-3)>>1
	 p45 = (143*i-3)>>1; p46 = (169*i-3)>>1; p47 = (187*i-3)>>1; p48 = (209*i-3)>>1
	 k = (210*i)>>1
	 while p1 < lndx
	     sieve[p1] = nil;   sieve[p2] = nil;   sieve[p3] = nil;   sieve[p4] = nil
	     sieve[p5] = nil;   sieve[p6] = nil;   sieve[p7] = nil;   sieve[p8] = nil
	     sieve[p9] = nil;   sieve[p10] = nil;  sieve[p11] = nil;  sieve[p12] = nil
	     sieve[p13] = nil;  sieve[p14] = nil;  sieve[p15] = nil;  sieve[p16] = nil
	     sieve[p17] = nil;  sieve[p18] = nil;  sieve[p19] = nil;  sieve[p20] = nil
	     sieve[p21] = nil;  sieve[p22] = nil;  sieve[p23] = nil;  sieve[p24] = nil
	     sieve[p25] = nil;  sieve[p26] = nil;  sieve[p27] = nil;  sieve[p28] = nil
	     sieve[p29] = nil;  sieve[p30] = nil;  sieve[p31] = nil;  sieve[p32] = nil
	     sieve[p33] = nil;  sieve[p34] = nil;  sieve[p35] = nil;  sieve[p36] = nil
	     sieve[p37] = nil;  sieve[p38] = nil;  sieve[p39] = nil;  sieve[p40] = nil
	     sieve[p41] = nil;  sieve[p42] = nil;  sieve[p43] = nil;  sieve[p44] =nil
	     sieve[p45] = nil;  sieve[p46] = nil;  sieve[p47] = nil;  sieve[p48] =nil
	     p1 += k; p2 += k; p3 += k; p4  += k; p5 += k; p6 += k; p7 += k
	     p8 += k; p9 += k; p10 +=k; p11 += k; p12 +=k; p13 +=k; p14 +=k
	     p15 +=k; p16 +=k; p17 +=k; p18 += k; p19 +=k; p20 +=k; p21 +=k
	     p22 +=k; p23 +=k; p24 +=k; p25 += k; p26 +=k; p27 +=k; p28 +=k
	     p29 +=k; p30 +=k; p31 +=k; p32 += k; p33 +=k; p34 +=k; p35 +=k
	     p36 +=k; p37 +=k; p38 +=k; p39 += k; p40 +=k; p41 +=k; p42 +=k
	     p43 +=k; p44 +=k; p45 +=k; p46 += k; p47 +=k; p48 +=k
	 end
      end
      return [2] if self < 3
      [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
   end
   
   def primesP60
      # all prime candidates > 5 are of form  60*k+(1,7,11,13,17,19,23,29,
      # 31,37,41,43,47,49,53,59
      # initialize sieve array with only these candidate values
      n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 = 1, 7, 11, 13, 17, 19, 23, 29, 31, 37
      n11, n12, n13, n14, n15, n16 = 41, 43, 47, 49, 53, 59
      num = self - self[0]^1;  sieve = []
      while n16 < num
         n1 += 60;   n2 += 60;   n3 += 60;   n4 += 60;   n5 += 60
	 n6 += 60;   n7 += 60;   n8 += 60;   n9 += 60;   n10 += 60
         n11 += 60;  n12 += 60;  n13 += 60;  n14 += 60;  n15 += 60;  n16 += 60
	 sieve[n1] = n1;  sieve[n2] = n2;  sieve[n3] = n3;  sieve[n4] = n4
	 sieve[n5] = n5;  sieve[n6] = n6;  sieve[n7] = n7;  sieve[n8] = n8
	 sieve[n9] = n9;  sieve[n10] = n10;  sieve[n11] = n11;  sieve[n12] = n12
	 sieve[n13] = n13; sieve[n14] = n14; sieve[n15] = n15;  sieve[n16] = n16
      end
      # now initialize sieve with the (odd) primes < 210, resize array
      sieve[2]=2;  sieve[3]=3;  sieve[5]=5;  sieve[7]=7;  sieve[11]=11;  sieve[13]=13
      sieve[17]=17;  sieve[19]=19;  sieve[23]=23;  sieve[29]=29;  sieve[31]=31
      sieve[37]=37; sieve[41]=41; sieve[43]=43; sieve[47]=47; sieve[53]=53; sieve[59]=59
      sieve = sieve[0..self]
      7.step(Math.sqrt(self).to_i, 2) do |i|
         next unless sieve[i]
	 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i,  k=30*i
	 # maps to: p1= (7*i-3)>>1, --- p8=(31*i)>>1, k=30*i>>1
	 #n6=6*i;  n7=n6+i;  n11=n6+n6-i;  n13=n6+n7;  n17=n11+n6
	 #n19=n13+n6;  n23=n17+n6;  n29=n23+n6;  n31=n29+(i<<1)
	 #n4=i<<2;  n8=i<<3;  n16=i<<4;  n7=n8-i;  n11=n7+n4;  n13=n11+(i<<1)
	 #n17=n16+i;  n19=n11+n8;  n23=n16+n7;  n29=n16+n13;  n31=(i<<5)-i
	 p1 = 7*i;  p2 = 11*i;  p3 = 13*i;  p4 = 17*i;  p5= 19*i; p6 = 23*i;  p7 = 29*i
	 p8 = 31*i;  p9 = 37*i; p10 = 41*i;  p11 = 43*i; p12 = 47*i; p13 = 53*i; p14 = 59*i
	 p15 = 49*i;  p16 = 61*i; k = 60*i
	 
	 #n2 = i<<1;  n4 = i<<2;  n8 = i<<3;  n16 = i<<4;  n32 = i<<5;  n64 = i<<6
	 #p1 = n8-i;  p2 = n8+n4-i;  p3 = n8+n4+i;  p4 = n16+i;  p5 = n16+n2+i;  p6 = n16+n8-i
	 #p7 = n32-n2-i;  p8 = n32-i;  p9=n32+n4+i;  p10=n32+n8+i;  p11=p2+n32;  p12=n32+n16-i
	 #p13 = p9+n16;  p14 = n64-n4-i;  p15 = n64-n16+i;  p16 = n64-n2-i;  k = n64-n4
	 while p1 < num
	     sieve[p1] = nil;  sieve[p2] = nil;   sieve[p3] = nil;  sieve[p4] = nil
	     sieve[p5] = nil;  sieve[p6] = nil;   sieve[p7] = nil;  sieve[p8] = nil
	     sieve[p9] = nil;  sieve[p10] = nil;  sieve[p11] = nil; sieve[p12] = nil
	     sieve[p13] = nil;  sieve[p14] = nil; sieve[p15] = nil; sieve[p16] = nil
	     p1 += k;  p2 += k;  p3 += k;  p4 += k;  p5 += k; p6 += k; p7 += k; p16 +=k
	     p8 += k;  p9 += k; p10 += k; p11 += k; p12 += k; p13 +=k; p14 +=k; p15 +=k
	 end
      end
      sieve.compact! 
   end
  
   def primesP60a
      # all prime candidates > 5 are of form  60*k+(1,7,11,13,17,19,23,29,
      # 31,37,41,43,47,49,53,59
      # 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, n3, n4, n5, n6, n7, n8, n9, n10 = -1, 2, 4, 5, 7, 8, 10, 13, 14, 17
      n11, n12, n13, n14, n15, n16 = 19, 20, 22, 23, 25, 28
      lndx= (self-1) >>1;  sieve = []
      while n16 < lndx
         n1 += 30;   n2 += 30;   n3 += 30;   n4 += 30;   n5 += 30
	 n6 += 30;   n7 += 30;   n8 += 30;   n9 += 30;   n10 += 30
         n11 += 30;  n12 += 30;  n13 += 30;  n14 += 30;  n15 += 30;  n16 += 30
	 sieve[n1] = n1;  sieve[n2] = n2;  sieve[n3] = n3;  sieve[n4] = n4
	 sieve[n5] = n5;  sieve[n6] = n6;  sieve[n7] = n7;  sieve[n8] = n8
	 sieve[n9] = n9;  sieve[n10] = n10;  sieve[n11] = n11;  sieve[n12] = n12
	 sieve[n13] = n13; sieve[n14] = n14; sieve[n15] = n15;  sieve[n16] = n16
      end
      # now initialize sieve with the (odd) primes < 210, resize array
      sieve[0]=0;  sieve[1]=1;  sieve[2]=2;  sieve[4]=4;  sieve[5]=5;  sieve[7]=7
      sieve[8]=8;  sieve[10]=10;  sieve[13]=13;  sieve[14]=14;  sieve[17]=17
      sieve[19]=19; sieve[20]=20; sieve[22]=22;  sieve[25]=25;  sieve[28]=28
      sieve = sieve[0..lndx-1]
      
      7.step(Math.sqrt(self).to_i, 2) do |i|
         next unless sieve[(i>>1)-1]
	 # p1=7*i, p2=11*i, p3=13*i ---- p14=53*i, p15=59*i, p16=61*i,  k=60*i
	 # maps to: p1= (7*i-3)>>1, --- p16=(61*i)>>1, k=60*i>>1 
	 n4=i<<2;  n8=i<<3;  n16=i<<4;  n32=i<<5;  n64=i<<6; n12=n16-n4
	 n19=n16+n4-i;  n37=n32+n4+i;  n41=n32+n8+i;  n43=43*i;  n48=n32+n16;  
	 n53=n64-n12+i;  n60=n64-n4
	 p1 = (n8-i-3)>>1; p2 = (n12-i-3)>>1;  p3 = (n12+i-3)>>1;  p4 = (n16+i-3)>>1
	 p5 = (n19-3)>>1;  p6 = (n16+n8-i-3)>>1; p7 = (n32-n4+i-3)>>1; p8=(n32-i-3)>>1
	 p9 = (n37-3)>>1;  p10 = (n41-3)>>1;   p11 = (n43-3)>>1;   p12 = (n48-i-3)>>1
	 p13 = (n48+i-3)>>1; p14 = (n53-3)>>1; p15 = (n60-i-3)>>1; p16 = (n60+i-3)>>1
	 k = n60>>1
	 while p1 < lndx
	     sieve[p1] = nil;  sieve[p2] = nil;  sieve[p3] = nil;  sieve[p4] = nil
	     sieve[p5] = nil;  sieve[p6] = nil;  sieve[p7] = nil;  sieve[p8] = nil
	     sieve[p9] = nil;  sieve[p10] = nil; sieve[p11] = nil; sieve[p12] = nil
	     sieve[p13] = nil; sieve[p14] = nil; sieve[p15] = nil; sieve[p16] = nil
	     p1 += k; p2 += k;  p3 += k;  p4 += k;  p5 += k; p6 += k; p7 += k; p8 += k
	     p9 += k; p10 += k; p11 += k; p12 += k; p13 +=k; p14 +=k; p15 +=k; p16 += k
	 end
      end
      return [2] if self < 3
      [2]+([nil]+sieve).compact!.map {|i| (i<<1)+3}
   end
   
   def primesP11a
      # all prime candidates > 11 are of form 2311k+(1, 13< resP11 < 2310)
      # uses generalized code with hardcoded arrays supplied, NOT OPTIMIZED!
      # initialize sieve array with only these candidate values
      # where sieve contains the odd integers representations
      # convert integers to array indices/vals by  i = (n-3)>>1 = (n>>1)-1
      # the primes and pcn array are converted values
  
      primes = [0, 1, 2, 4, 5, 7, 8, 10, 13, 14, 17, 19, 20, 22, 25, 28, 29, 32, 34, 35, 38, 40, 43, 
      47, 49, 50, 52, 53, 55, 62, 64, 67, 68, 73, 74, 77, 80, 82, 85, 88, 89, 94, 95, 97, 98, 104, 
      110, 112, 113, 115, 118, 119, 124, 127, 130, 133, 134, 137, 139, 140, 145, 152, 154, 155, 157, 
      164, 167, 172, 173, 175, 178, 182, 185, 188, 190, 193, 197, 199, 203, 208, 209, 214, 215, 218, 
      220, 223, 227, 229, 230, 232, 238, 242, 244, 248, 250, 253, 259, 260, 269, 272, 277, 280, 283, 
      284, 287, 292, 295, 298, 299, 302, 305, 307, 308, 314, 319, 320, 322, 325, 328, 329, 335, 337, 
      340, 344, 349, 353, 358, 362, 365, 368, 370, 374, 377, 379, 383, 385, 392, 397, 403, 404, 409,
      410, 412, 413, 418, 425, 427, 428, 430, 437, 439, 440, 442, 452, 454, 458, 463, 467, 469, 472, 
      475, 482, 484, 487, 490, 494, 497, 503, 505, 508, 509, 514, 515, 518, 523, 524, 529, 530, 533, 
      542, 544, 545, 547, 550, 553, 557, 560, 563, 574, 575, 580, 584, 589, 592, 595, 599, 605, 607, 
      610, 613, 614, 617, 623, 628, 637, 638, 640, 643, 644, 647, 649, 650, 652, 658, 659, 662, 679, 
      682, 685, 689, 698, 703, 710, 712, 713, 715, 718, 722, 724, 725, 728, 734, 739, 740, 742, 743, 
      745, 748, 754, 760, 764, 770, 773, 775, 778, 782, 784, 788, 790, 797, 799, 802, 803, 805, 808, 
      809, 812, 817, 827, 830, 832, 833, 845, 847, 848, 853, 859, 860, 865, 869, 872, 875, 878, 887, 
      890, 892, 893, 899, 904, 910, 914, 922, 929, 932, 934, 935, 937, 938, 943, 949, 952, 955, 964, 
      965, 973, 974, 985, 988, 992, 995, 997, 998, 1000, 1004, 1007, 1012, 1013, 1018, 1025, 1030, 
      1033, 1039, 1040, 1042, 1043, 1048, 1054, 1055, 1063, 1064, 1067, 1069, 1070, 1075, 1079, 1088, 
      1100, 1102, 1105, 1109, 1117, 1118, 1120, 1124, 1132, 1133, 1135, 1139, 1142, 1145, 1147, 1153]
      
      residues =[13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 
      103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 169, 173, 179, 181, 191, 193, 
      197, 199, 211, 221, 223, 227, 229, 233, 239, 241, 247, 251, 257, 263, 269, 271, 277, 281, 283, 
      289, 293, 299, 307, 311, 313, 317, 323, 331, 337, 347, 349, 353, 359, 361, 367, 373, 377, 379, 
      383, 389, 391, 397, 401, 403, 409, 419, 421, 431, 433, 437, 439, 443, 449, 457, 461, 463, 467, 
      479, 481, 487, 491, 493, 499, 503, 509, 521, 523, 527, 529, 533, 541, 547, 551, 557, 559, 563, 
      569, 571, 577, 587, 589, 593, 599, 601, 607, 611, 613, 617, 619, 629, 631, 641, 643, 647, 653, 
      659, 661, 667, 673, 677, 683, 689, 691, 697, 701, 703, 709, 713, 719, 727, 731, 733, 739, 743, 
      751, 757, 761, 767, 769, 773, 779, 787, 793, 797, 799, 809, 811, 817, 821, 823, 827, 829, 839, 
      841, 851, 853, 857, 859, 863, 871, 877, 881, 883, 887, 893, 899, 901, 907, 911, 919, 923, 929, 
      937, 941, 943, 947, 949, 953, 961, 967, 971, 977, 983, 989, 991, 997, 1003, 1007, 1009, 1013, 
      1019, 1021, 1027, 1031, 1033, 1037, 1039, 1049, 1051, 1061, 1063, 1069, 1073, 1079, 1081, 1087,
      1091, 1093, 1097, 1103, 1109, 1117, 1121, 1123, 1129, 1139, 1147, 1151, 1153, 1157, 1159, 1163, 
      1171, 1181, 1187, 1189, 1193, 1201, 1207, 1213, 1217, 1219, 1223, 1229, 1231, 1237, 1241, 1247, 
      1249, 1259, 1261, 1271, 1273, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1313, 1319, 
      1321, 1327, 1333, 1339, 1343, 1349, 1357, 1361, 1363, 1367, 1369, 1373, 1381, 1387, 1391, 1399, 
      1403, 1409, 1411, 1417, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1457, 1459, 1469, 1471, 
      1481, 1483, 1487, 1489, 1493, 1499, 1501, 1511, 1513, 1517, 1523, 1531, 1537, 1541, 1543, 1549, 
      1553, 1559, 1567, 1571, 1577, 1579, 1583, 1591, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 
      1633, 1637, 1643, 1649, 1651, 1657, 1663, 1667, 1669, 1679, 1681, 1691, 1693, 1697, 1699, 1703, 
      1709, 1711, 1717, 1721, 1723, 1733, 1739, 1741, 1747, 1751, 1753, 1759, 1763, 1769, 1777, 1781, 
      1783, 1787, 1789, 1801, 1807, 1811, 1817, 1819, 1823, 1829, 1831, 1843, 1847, 1849, 1853, 1861, 
      1867, 1871, 1873, 1877, 1879, 1889, 1891, 1901, 1907, 1909, 1913, 1919, 1921, 1927, 1931, 1933, 
      1937, 1943, 1949, 1951, 1957, 1961, 1963, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 
      2021, 2027, 2029, 2033, 2039, 2041, 2047, 2053, 2059, 2063, 2069, 2071, 2077, 2081, 2083, 2087, 
      2089, 2099, 2111, 2113, 2117, 2119, 2129, 2131, 2137, 2141, 2143, 2147, 2153, 2159, 2161, 2171, 
      2173, 2179, 2183, 2197, 2201, 2203, 2207, 2209, 2213, 2221, 2227, 2231, 2237, 2239, 2243, 2249, 
      2251, 2257, 2263, 2267, 2269, 2273, 2279, 2281, 2287, 2291, 2293, 2297, 2309]
      
      pcn = [-1, 5, 7, 8, 10, 13, 14, 17, 19, 20, 22, 25, 28, 29, 32, 34, 35, 38, 40, 43, 47, 49, 50, 
      52, 53, 55, 62, 64, 67, 68, 73, 74, 77, 80, 82, 83, 85, 88, 89, 94, 95, 97, 98, 104, 109, 110, 
      112, 113, 115, 118, 119, 122, 124, 127, 130, 133, 134, 137, 139, 140, 143, 145, 148, 152, 154, 
      155, 157, 160, 164, 167, 172, 173, 175, 178, 179, 182, 185, 187, 188, 190, 193, 194, 197, 199, 
      200, 203, 208, 209, 214, 215, 217, 218, 220, 223, 227, 229, 230, 232, 238, 239, 242, 244, 245, 
      248, 250, 253, 259, 260, 262, 263, 265, 269, 272, 274, 277, 278, 280, 283, 284, 287, 292, 293, 
      295, 298, 299, 302, 304, 305, 307, 308, 313, 314, 319, 320, 322, 325, 328, 329, 332, 335, 337, 
      340, 343, 344, 347, 349, 350, 353, 355, 358, 362, 364, 365, 368, 370, 374, 377, 379, 382, 383, 
      385, 388, 392, 395, 397, 398, 403, 404, 407, 409, 410, 412, 413, 418, 419, 424, 425, 427, 428, 
      430, 434, 437, 439, 440, 442, 445, 448, 449, 452, 454, 458, 460, 463, 467, 469, 470, 472, 473, 
      475, 479, 482, 484, 487, 490, 493, 494, 497, 500, 502, 503, 505, 508, 509, 512, 514, 515, 517, 
      518, 523, 524, 529, 530, 533, 535, 538, 539, 542, 544, 545, 547, 550, 553, 557, 559, 560, 563, 
      568, 572, 574, 575, 577, 578, 580, 584, 589, 592, 593, 595, 599, 602, 605, 607, 608, 610, 613, 
      614, 617, 619, 622, 623, 628, 629, 634, 635, 637, 638, 640, 643, 644, 647, 649, 650, 652, 655, 
      658, 659, 662, 665, 668, 670, 673, 677, 679, 680, 682, 683, 685, 689, 692, 694, 698, 700, 703, 
      704, 707, 710, 712, 713, 715, 718, 722, 724, 725, 727, 728, 733, 734, 739, 740, 742, 743, 745, 
      748, 749, 754, 755, 757, 760, 764, 767, 769, 770, 773, 775, 778, 782, 784, 787, 788, 790, 794, 
      797, 799, 802, 803, 805, 808, 809, 812, 815, 817, 820, 823, 824, 827, 830, 832, 833, 838, 839, 
      844, 845, 847, 848, 850, 853, 854, 857, 859, 860, 865, 868, 869, 872, 874, 875, 878, 880, 883, 
      887, 889, 890, 892, 893, 899, 902, 904, 907, 908, 910, 913, 914, 920, 922, 923, 925, 929, 932, 
      934, 935, 937, 938, 943, 944, 949, 952, 953, 955, 958, 959, 962, 964, 965, 967, 970, 973, 974, 
      977, 979, 980, 985, 988, 992, 995, 997, 998, 1000, 1004, 1007, 1009, 1012, 1013, 1015, 1018, 
      1019, 1022, 1025, 1028, 1030, 1033, 1034, 1037, 1039, 1040, 1042, 1043, 1048, 1054, 1055, 1057, 
      1058, 1063, 1064, 1067, 1069, 1070, 1072, 1075, 1078, 1079, 1084, 1085, 1088, 1090, 1097, 1099, 
      1100, 1102, 1103, 1105, 1109, 1112, 1114, 1117, 1118, 1120, 1123, 1124, 1127, 1130, 1132, 1133, 
      1135, 1138, 1139, 1142, 1144, 1145, 1147, 1153]

      # now initialize modPn, lndx and sieve then find prime candidates
      # set initial prime candidates to (converted) residue values
      lndx= (self-1)>>1;   mod = 2310;  m = mod>>1;   sieve = []
      while pcn.last < lndx
	  pcn.each_index {|j| n = pcn[j]+m;  sieve[n] = n;  pcn[j] = n }
      end
      # now initialize sieve with the (odd) primes < modPn, resize array
      primes.each {|p| sieve[p] = p };   sieve = sieve[0..lndx-1]

      # perform sieve with residue elements resPn+(modPn+1)
      residues << (mod + 1)
      13.step(Math.sqrt(self).to_i, 2) do |i|
          next unless sieve[(i>>1)-1]
	  # p1=res1*i, p2=res2*i, p3=res3*i -- pN-1=res(N-1)*i, pN=resN*i, k=mod*i
	  # maps to: p1= (res1*i-3)>>1, --- pN=(((mod+1)-3)*i)>>1, k=m*i>>1
          res = residues.map{|prm| (prm*i>>1)-1};   k = (mod*i)>>1
          while res[0] < lndx
	      res.each_index {|j| sieve[res[j]] = nil;  res[j] += k }
	  end
      end
      return [2] if self < 3
      [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
   end

   def primesPn(input=nil)
      # generalized Sieve of Zakiya, takes an input prime or modulus
      # all prime candidates > Pn are of form modPn*k+[1, res(Pn)]
      # initialize sieve array with only these candidate values
      # where sieve contains the odd integers representations
      # convert integers to array indices/vals by  i = (n-3)>>1 = (n>>1)-1

      # if no input value then just perform SoZ P7a as reference sieve
      return  primesP7a  if  input == nil

      seeds = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61]
      mcps = [289, 391, 481, 493, 529, 559, 629, 667, 793, 841, 893, 899, 923, 961, 
              1037, 1121, 1189, 1273, 1349, 1387, 1411, 1417, 1469, 1517, 1643, 1681, 
              1751, 1781, 1817, 1829, 1919, 2021]
		  
      # if input is a prime number, determine modulus = P!(n) and seed primes
      if input&1 == 1
          # find seed primes <= Pn,  compute modPn
	  return 'INVALID OR TOO LARGE PRIME' if !seeds.include? input
          seeds = seeds[0 .. seeds.index(input)];   mod = seeds.inject {|a,b| a*b}
      else  # if input a modulus (even number), determine seed primes for it
	   md,  pp, mod = 2, 3, input
	   seeds[1..-1].each {|i| md *=i;  break if mod < md;  pp = i}
	   seeds = seeds[0 .. seeds.index(pp)]
      end

      # create array of prime residues for primes Pn < pk < modPn
      primes = mod.primesP7a;    residues  = primes - seeds

      # find modular compliments of primes [1, Pn < p < modPn] if necessary
      mcp = [];  tmp = [1]+residues
      while !tmp.empty?;  b = mod - tmp.shift;  mcp << b unless tmp.delete(b)  end

      residues.concat(mcp).sort!;  residues.concat(mcps).sort! if [11, 2310].include? input
      
      # now initialize modPn, lndx and sieve then find prime candidates
      # set initial prime candidates to (converted) residue values
      lndx= (self-1)>>1;   m=mod>>1;   sieve = []
      pcn = [-1]+residues.map {|t| (t>>1)-1}
      while pcn.last < lndx
	  pcn.each_index {|j| n = pcn[j]+m;  sieve[n] = n;  pcn[j] = n }
      end
      # now initialize sieve with the (odd) primes < modPn, resize array
      primes[1..-1].each {|x| p =(x-3)>>1;  sieve[p] = p };   sieve = sieve[0..lndx-1]

      # perform sieve with residue elements resPn+(modPn+1)
      residues << (mod + 1);   res0 = residues[0]
      res0.step(Math.sqrt(self).to_i, 2) do |i|
          next unless sieve[(i>>1)-1]
	  # p1=res1*i, p2=res2*i, p3=res3*i -- pN-1=res(N-1)*i, pN=resN*i, k=mod*i
	  # maps to: p1= (res1*i-3)>>1, --- pN=(((mod+1)-3)*i)>>1, k=mod*i>>1
          res = residues.map{|prm| (prm*i-3)>>1};   k = (mod*i)>>1
          while res[0] < lndx
	      res.each_index {|j| sieve[res[j]] = nil;  res[j] += k }
	  end
      end
      return [2] if self < 3
      [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
   end

   def primz?
      # number theory based deterministic primality tester
      n = self.abs
      return true  if [2, 3, 5].include? n
      return false if n == 1 || n & 1 == 0
      return false if n > 5 && ( ! [1, 5].include?(n%6) || n%5 == 0)

      7.step(Math.sqrt(n).to_i,2) do |i|
	 #  p5= 5*i,  k = 6*i,  p7 = 7*i  
	 p1 = 5*i;  k = p1+i;  p2 = k+i
	 return false if  [(n-p1)%k , (n-p2)%k].include? 0
      end   	
      return true
   end
end


   def sieveOfAtkin(num)

      num += 1
      lng = (num>>1)-1+(num&1)
      sieve = [nil]*(lng + 1)

      x_max, x2  = Math.sqrt((num-1) / 4.0).to_i, 0
      4.step( 8*x_max + 1, 8) do |xd|
          x2 += xd;    y_max = Math.sqrt(num-x2).to_i
          n, n_diff = x2 + y_max**2, (y_max << 1) - 1
          if n&1 == 0  then  n -= n_diff;   n_diff -= 2  end
          ((n_diff - 1) << 1).step( -2, -8)  do |d|
              m = n%12
              if (m == 1 or m == 5) then m= n >> 1;  sieve[m] = !sieve[m]  end
	      #if [1,5].include? m then m = n >> 1;  sieve[m] = !sieve[m]  end
              n -= d
	  end
      end

      x_max,  x2 = Math.sqrt((num-1) / 3.0).to_i, 0
      3.step(6*x_max + 1, 6) do |xd|
          x2 += xd;    y_max = Math.sqrt(num-x2).to_i
          n,  n_diff = x2 + y_max**2, (y_max << 1) - 1
          if n&1 == 0  then  n -= n_diff;   n_diff -= 2  end
          ((n_diff - 1) << 1).step( -2, -8) do |d|
              if (n%12 == 7) then  m = n >> 1;  sieve[m] = !sieve[m]  end
              n -= d
	  end
      end

      x_max, y_min, x2, xd = ((2 + Math.sqrt(4-8*(1-num))) / 4).to_i, -1, 0, 3
      (1 ... x_max + 1).each do |x|
          x2 += xd;   xd += 6
	  y_min = (((Math.sqrt(x2 - num).ceil - 1) << 1) - 2) << 1  if x2 >= nu
          n, n_diff = ((x*x + x) << 1) - 1, (((x-1) << 1) - 2) << 1
          (n_diff).step(y_min-1, -8) do |d|
              if (n%12 == 11)  then  m = n >> 1;  sieve[m] = !sieve[m]  end
              n += d
	  end
      end

      return [2]  if num-1 < 3
      primes = [2,3]

      (2 ...  (Math.sqrt(num).to_i+ 1) >> 1).each do |n|
          next unless sieve[n]
	  i = (n << 1) +1;   j= i*i;   primes << i
          (j).step(num-1, j << 1) { |k| sieve[k>>1] = nil }
      end

      s = Math.sqrt(num).to_i + 1
      s += 1  if  s&1 == 0
      s.step(num-1, 2) {|i| primes << i  if sieve[i>>1] }

      return primes
    end


def tm(code); s=Time.now; eval code; Time.now-s end
def primesxy(x,y);  (x..y).map {|p| p.primz? ? p : nil }.compact!  end

limit= 10_000_001
ret1, ret2, ret3, ret4, ret5, ret6, ret7, ret8, ret9 = nil

Benchmark.bm(14) do |t| 
   t.report("primes up to #{limit}: ") do ret1 = sieveOfAtkin(limit)end
   t.report("primes up to #{limit}: ") do ret2 = limit.primesP60a end
   t.report("primes up to #{limit}: ") do ret3 = limit.primesPn(3) end
   t.report("primes up to #{limit}: ") do ret4 = limit.primesPn(5) end
   t.report("primes up to #{limit}: ") do ret5 = limit.primesPn(7) end
   t.report("primes up to #{limit}: ") do ret6 = limit.primesPn(11) end
   t.report("primes up to #{limit}: ") do ret7 = limit.primesP3a end
   t.report("primes up to #{limit}: ") do ret8 = limit.primesP5a end
   t.report("primes up to #{limit}: ") do ret9 = limit.primesP7a end
end

#p ret1.first(10)
puts;  p ret1.last(10);  p ret1.size

#p ret2.first(10)
puts;  p ret2.last(10);  p ret2.size

#p ret3.first(10)
puts;  p ret3.last(10);  p ret3.size

#p ret4.first(10)
puts;  p ret4.last(10);  p ret4.size

#p ret5.first(10)
puts;  p ret5.last(10);  p ret5.size

#p ret6.first(10)
puts;  p ret6.last(10);  p ret6.size

#p ret7.first(10)
puts;  p ret7.last(10);  p ret7.size

#p ret8.first(10)
puts;  p ret8.last(10);  p ret8.size

#p ret9.first(10)
puts;  p ret9.last(10);  p ret9.size