DZone Snippets: number code

Sieve of Zakiya -- Number Theory Sieves (NTS)

jzakiya (Jabari Zakiya) — Fri, 06 Jun 2008 22:47:09 GMT

// 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

First n primes in Ruby

ntk () — Thu, 17 Apr 2008 20:32:45 GMT





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



require 'benchmark' 



class Integer    

   def prime?         # cf. http://snippets.dzone.com/posts/show/4636

     n = self.abs()

     return true if n == 2

     return false if n == 1 || n & 1 == 0

     d = n-1

     d >>= 1 while d & 1 == 0

     20.times do                               # 20 = k from above

       a = rand(n-2) + 1

       t = d

       y = ModMath.pow(a,t,n)                  # implemented below

       while t != n-1 && y != 1 && y != n-1

         y = (y * y) % n

         t <<= 1

       end

       return false if y != n-1 && t & 1 == 0

     end

     return true

   end

end

 

module ModMath

   def ModMath.pow(base, power, mod)

     result = 1

     while power > 0

       result = (result * base) % mod if power & 1 == 1

       base = (base * base) % mod

       power >>= 1;

     end

     result

   end

end







class Integer



   def primes   # cf. http://snippets.dzone.com/posts/show/3734



      sieve = []

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

      sieve[1] = nil

      sieve[2] = 2



      3.step(Math.sqrt(self).floor, 2) do |i| 

         next unless sieve[i]

         (i*i).step(self, i) do |j|

            sieve[j] = nil

         end

      end



      sieve.compact!



   end 





   def primes2       # cf. http://snippets.dzone.com/posts/show/3734



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



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

         next unless primes[i]

            (i*i).step(self, i) do |j|

               primes[j] = nil

            end

      end

	

      primes.compact!



   end



end







class Integer



   @primes_cache ||= 100.primes

   #@primes_cache ||= [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]

   class << self; attr_accessor :primes_cache; end



   @nthprime_limit ||= 5_000_000

   class << self; attr_reader :nthprime_limit; end

   #class << self; attr_accessor :nthprime_limit; end





   def sieve_size

      num = self.to_i.abs

      logn = Math.log(num)

      if num < 15985     # cf. http://primes.utm.edu/howmany.shtml

         ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor

      else

         ( num * (logn + Math.log(logn) - 0.9427) ).floor

      end

   end





   def nthprime



      num = self.to_i.abs



      # set a limit; cf. http://primes.utm.edu/nthprime/: The 5,000,000th prime is 86,028,121.

      raise "#{num}: number too big for Integer#nthprime" if num > Integer.nthprime_limit      



      primes_cache_size = Integer.primes_cache.size



      if num > primes_cache_size



         # optional; cf. Integer#nthprime_add and Integer#nthprime_add_mr below

         if num >= 50_000 && num < 100_000 && (num - primes_cache_size) < 5000 then return num.nthprime_add end

         if num >= 100_000 && num < 200_000 && (num - primes_cache_size) < 15000 then return num.nthprime_add end

         if num >= 200_000 && (num - primes_cache_size) < 35000 then return num.nthprime_add end

         

         logn = Math.log(num)



         if num < 15985       # cf. http://primes.utm.edu/howmany.shtml

            limit = ( num * (logn + Math.log(logn) - 1 + 1.8 * Math.log(logn) / logn) ).floor

            Integer.primes_cache = limit.primes

         else

            limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor

            Integer.primes_cache = limit.primes

         end



=begin

         elsif num >= 15985 && num <= 1_000_000

            limit = ( num * (logn + Math.log(logn) - 0.9427) ).floor

            Integer.primes_cache = limit.primes

         elsif num > 1_000_000

            Integer.primes_cache = 20_000_000.primes

         end

=end



      else

         return Integer.primes_cache.at(num-1)

      end





      if num <= Integer.primes_cache.size

         return Integer.primes_cache.at(num-1)

      end



      if Integer.primes_cache.size > 2_500_000

         num.nthprime_add_mr 

      else

         num.nthprime_add   # faster for a (relatively) small prime cache

      end



   end





   def nthprime_add       #  add primes to Integer.primes_cache up to the nth prime



      num = self.to_i.abs



      if num <= Integer.primes_cache.size

         return Integer.primes_cache.at(num-1)

      end



      last_prime = Integer.primes_cache.last

      last_prime_divmod = last_prime.divmod(6)



      if last_prime_divmod.last == 1

         i = last_prime_divmod.first

         Integer.primes_cache.pop      # avoid a duplicate prime

      else

         i = last_prime_divmod.first + 1

      end





      while Integer.primes_cache.size < num



         n1 = 6*i+1       # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and

         n2 = n1+4        # http://everything2.com/index.pl?node_id=1176369

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0



            next_num_sqrt = Math.sqrt(p).floor



            Integer.primes_cache.each do |d| 

               if d > next_num_sqrt 

                  Integer.primes_cache << p

                  #print "\r\e[0K#{Integer.primes_cache.size}"

                  #$stdout.flush

                  break

               elsif p % d == 0

                  break

               end  

            end

         end   # each  



      end  # while



      return Integer.primes_cache.at(num-1)



   end





   def nthprime_add_mr       #  add Miller-Rabin primes to Integer.primes_cache up to the nth prime



      num = self.to_i.abs



      if num <= Integer.primes_cache.size

         return Integer.primes_cache.at(num-1)

      end



      last_prime = Integer.primes_cache.last

      last_prime_divmod = last_prime.divmod(6)



      if last_prime_divmod.last == 1

         i = last_prime_divmod.first

         Integer.primes_cache.pop

      else

         i = last_prime_divmod.first + 1

      end





      while Integer.primes_cache.size < num



         search_next_prime = true



         while search_next_prime



            n1 = 6*i+1       

            n2 = n1+4        

            i += 1



            [n1, n2].each do |p| 



               next if p % 5 == 0 || p % 7 == 0



               if p.prime?

                  Integer.primes_cache << p

                  search_next_prime = false

                  #print "\r\e[0K#{Integer.primes_cache.size}"

                  #$stdout.flush

               end

            end

         end



      end  #  first while loop



      return Integer.primes_cache.at(num-1)



   end





#------------------------------------------- some additional prime methods





   def next_primes_in_cache    # next primes in Integer.primes_cache



      search_next_primes = self.to_i.abs



      last_prime = Integer.primes_cache.last

      last_prime_divmod = last_prime.divmod(6)



      if last_prime_divmod.last == 1

         i = last_prime_divmod.first

         Integer.primes_cache.pop

      else

         i = last_prime_divmod.first + 1

      end



      while search_next_primes > 0



         n1 = 6*i+1       

         n2 = n1+4        

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0

            next_num_sqrt = Math.sqrt(p).floor



            Integer.primes_cache.each do |d| 

               if d > next_num_sqrt 

                  Integer.primes_cache << p

                  search_next_primes -= 1

                  #print "\r\e[0K#{Integer.primes_cache.size}"

                  #$stdout.flush

                  break

               elsif p % d == 0

                  break

               end  

            end

         end 

      end   # while



      Integer.primes_cache.last(self.to_i.abs)



   end





   def next_mr_prime     # next Miller-Rabin prime

     

      last_num = self.to_i.abs

      last_num_divmod = last_num.divmod(6)



      if last_num_divmod.last == 1

         i = last_num_divmod.first

      else

         i = last_num_divmod.first + 1

      end



      next_prime = nil

      search_next_prime = true



      while search_next_prime



         n1 = 6*i+1       

         n2 = n1+4        

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0



            if p > last_num && p.prime?

               next_prime = p

               search_next_prime = false 

               break

            end

         end

      end



         next_prime



   end







   def next_mr_primes(n)     # next Miller-Rabin primes



      last_num = self.to_i.abs

      last_num_divmod = last_num.divmod(6)



      if last_num_divmod.last == 1

         i = last_num_divmod.first

      else

         i = last_num_divmod.first + 1

      end



      next_primes = []

      search_next_primes = n.to_i.abs



      while search_next_primes > 0



         n1 = 6*i+1      

         n2 = n1+4       

         i += 1



         [n1, n2].each do |p| 



            next if p % 5 == 0 || p % 7 == 0



            if p > last_num && p.prime?

               next_primes << p

               search_next_primes -= 1 

            end

         end

      end



      next_primes



   end



end







num1 = 10_000

num1 = 1_000

num1 = 5_000



num2 = 210_349

num2 = 100_125

num2 = 55_000



ret1 = nil

ret2 = nil

ret3 = nil

ret4 = nil





Benchmark.bm(16) do |t| 



   t.report("first #{num1} primes: ") do

      ret1 = num1.nthprime_add

   end 



   Integer.primes_cache.clear

   Integer.primes_cache.concat(100.primes)



   t.report("first #{num1} primes: ") do

      ret2 = num1.nthprime_add_mr

   end 



   Integer.primes_cache.clear

   Integer.primes_cache.concat(100.primes)



   t.report("first #{num1} primes: ") do

      ret3 = num1.nthprime

   end 



   Integer.primes_cache.clear

   Integer.primes_cache.concat(100.primes)



   t.report("first #{num2} primes: ") do

      ret4 = num2.nthprime

   end 



end





puts

puts "the #{num1}th prime number: #{ret1}"

puts "the #{num1}th prime number: #{ret2}"

puts "the #{num1}th prime number: #{ret3}"

puts "the #{num2}th prime number: #{ret4}"

puts

puts "the #{num1}th prime: #{Integer.primes_cache.at(num1-1)}"

puts "the #{num2}th prime: #{Integer.primes_cache.at(num2-1)}"

puts

p Integer.primes_cache.first(10)

p Integer.primes_cache.last(10)

p Integer.primes_cache.size

puts





p 15.next_primes_in_cache



puts 594_213.next_mr_prime



p 149_137.next_mr_primes(5)



puts

p 1_000.sieve_size

p 1_000_000.sieve_size

p 2_500_000.sieve_size

p 5_000_000.sieve_size

p 10_000_000.sieve_size

p 100_000_000.sieve_size





=begin



# check with primegen.c, http://cr.yp.to/primegen.html

primes 2 48611 | nl | tail -n 1

primes 2 104729 | nl | tail -n 1

primes 2 679277 | nl | tail -n 1

primes 2 1301423 | nl | tail -n 1

primes 2 2903603 | nl | tail -n 1

primes 1303129 1303283 | nl

primes 594213 $((594213 + 500)) | head -n 1

primes 149137 $((149137 + 1000)) | head -n 5



# further prime tools from http://libtom.org

- LibTomMath

   bn_mp_prime_is_prime.c

   bn_mp_prime_miller_rabin.c

- TomsFastMath

   fp_isprime.c

   fp_prime_miller_rabin.c



=end

Print a binary number in C

scvalex (Alexandru Scvortov) — Tue, 30 Oct 2007 11:53:18 GMT

These are two functions that print the binary representation of an integer. The first simply prints it out, while the second only prints out the relevant digits (i.e. cuts the first x 0 digits) in groups of four.

Explanation here.



#include 



/* Print n as a binary number */

void printbitssimple(int n) {

	unsigned int i;

	i = 1<<(sizeof(n) * 8 - 1);



	while (i > 0) {

		if (n & i)

			printf("1");

		else

			printf("0");

		i >>= 1;

	}

}



/* Print n as a binary number */

void printbits(int n) {

	unsigned int i, step;



	if (0 == n) { /* For simplicity's sake, I treat 0 as a special case*/

		printf("0000");

		return;

	}



	i = 1<<(sizeof(n) * 8 - 1);



	step = -1; /* Only print the relevant digits */

	step >>= 4; /* In groups of 4 */

	while (step >= n) {

		i >>= 4;

		step >>= 4;

	}



	/* At this point, i is the smallest power of two larger or equal to n */

	while (i > 0) {

		if (n & i)

			printf("1");

		else

			printf("0");

		i >>= 1;

	}

}



int main(int argc, char *argv[]) {

	int i;

	for (i = 0; i < 16; ++i) {

		printf("%d = ", i);

		printbitssimple(i);

		printf("\n");

	}



	return 0;

}

Random numbers in Ruby

ntk () — Wed, 24 Oct 2007 19:01:38 GMT





# Pseudo random number generator

# From: http://gnuvince.wordpress.com/2005/11/24/pseudo-random-number-generator

# Author: Vincent Foley-Bourgon



def random_number

  t = Time.now.to_f / (Time.now.to_f % Time.now.to_i)

  random_seed = t * 1103515245 + 12345;

  (random_seed / 65536) % 32768;

end



10.times { puts random_number }



cnt = Array.new(10,0)

20000.times { cnt[(random_number % 10).to_i] += 1 }

p cnt

puts





# Ruby Gaussian Random Number Generator

# Author: Glenn

# http://webhost101.net/rails/typo/articles/2007/07/31/ruby-gaussian-random-number-generator



def gaussian_rand 

   u1 = u2 = w = g1 = g2 = 0  # declare

   begin

     u1 = 2 * rand - 1

     u2 = 2 * rand - 1

     w = u1 * u1 + u2 * u2

   end while w >= 1

   

   w = Math::sqrt( ( -2 * Math::log(w)) / w )

   g2 = u1 * w;

   g1 = u2 * w;

   # g1 is returned  

end



sum = 0

sumsq = 0

n = 1000

0.upto(n) do 

  #r = gaussian_rand

  r = gaussian_rand * 100 + 50   # new_random_number = gaussian_rand * standard_deviation + average

  #puts r

  sum += r

  sumsq += (r*r)

end



ave = sum/n

stddev = Math::sqrt(sumsq/n - ave * ave)

puts "Average = #{ave}"

puts "StdDev = #{stddev}"

puts





# Random Number Generator

# http://scutigena.sourceforge.net/test-random.html

# http://scutigena.sourceforge.net/sources/ruby-1.7.2/random.html

# $Id: random.ruby,v 1.2 2003/12/30 01:25:05 davidw Exp $



IM = 139968

IA = 3877

IC = 29573



$last = 42.0

def gen_random (max) (max * ($last = ($last * IA + IC) % IM)) / IM end



10.times do

    puts gen_random(100000.0)

end



printf "%.5f\n", gen_random(100.0)

puts





# From: http://matt.blogs.it/entries/00002641.html

# Author: Matt Mower

# cf. http://www.taygeta.com/random/gaussian.html

# cf. http://www.bearcave.com/misl/misl_tech/wavelets/hurst/random.html



def box_mueller( mean = 0.0, stddev = 1.0 )

  x1 = 0.0, x2 = 0.0, w = 0.0



  until w > 0.0 && w < 1.0

    x1 = 2.0 * rand - 1.0

    x2 = 2.0 * rand - 1.0

    w = ( x1 * x2 ) + ( x2 * x2 )

  end



  w = Math.sqrt( -2.0 * Math.log( w ) / w )

  r = x1 * w



  mean + r * stddev

end



10.times { puts box_mueller(5.0, 1.0) }

puts





# Generating random numbers with a specified distribution

# From: http://www.cs.nyu.edu/~michaels/blog/?p=24

# Author: Michael Schidlowsky



def gen

  (x=rand)>0.5 ? x : (x < rand/2.0) ? 1.0 - x : x

end



def gen2

   (x = rand) && rand ? 1.0 - x : x

end



def compute_distribution(numSamples, &block)

   samples = []

   values = 10

   numSamples.times{samples << yield}

   dist = Array.new(values, 0)

   samples.each{ |s| dist[(s*values).floor] += 1 }

   dist

end



p compute_distribution(1000){rand}

p compute_distribution(1000){gen}

p compute_distribution(1000) { gen2 }

puts





# Random number generation for a triangular distribution

# From: http://www.brighton-webs.co.uk/distributions/triangular.asp

# cf. http://en.wikipedia.org/wiki/Triangular_distribution



def rng(m, low, high)



   # cf. the parameter info at http://www.brighton-webs.co.uk/distributions/triangular.asp

   return nil unless high > low && m > low && m < high    



   u = rand



   if u <= (m-low)/(high-low)

      r = low+ Math.sqrt(u*(high-low)*(m-low))

   else

      r = high - Math.sqrt((1.0-u)*(high-low)*(high-m))

   end



end



10.times do

  #puts rng(0.0, -25.0, 25.0)

  puts rng(50.0, 25.0, 75.0)

end

puts





# From: http://www.ruby-forum.com/topic/95931

# Author: Christoffer Lernö



class Range

  def rand

    return first if exclude_end? && last == first

    Kernel::rand(last - first + (exclude_end? ? 0 : 1)) + first

  end

end



r1 = -9..9

r2 = -90...100

p r1.rand, r2.rand



r3 = 0..9



num = []

10.times do

   num << r3.rand

   if num.first.zero? then num.shift; redo end

end



p num

puts num.to_s.to_i

puts





# See:

# http://redcorundum.blogspot.com/2008/01/randomizing-array-revisited.html

# http://redcorundum.blogspot.com/2006/12/randomizing-array-and-other-faqs.html

# http://szeryf.wordpress.com/2007/06/19/a-simple-shuffle-that-proved-not-so-simple-after-all/



class Array

  def shuffle

    array = dup

    size.downto 2 do |j|

      r = rand j

      array[j-1], array[r] = array[r], array[j-1]

    end

    array

  end

end



array = (1..50).to_a



10.times { p array.shuffle.first(10) }





#-------------------





class RNG



   def num(min=8,max=min+5,iter=1)



      return nil unless min.is_a?(Integer) && max.is_a?(Integer) && max > min && iter.is_a?(Integer)



      ret = []

      stats = Hash.new(0)     # optional; cf. stats[random_num] += 1

      digits = Array(0..9).sort_by {rand}

      #digits = (Array(0..9) * (rand(4)+1)).sort_by {rand}

      digits_size = digits.size



      iter.times do

         count = 0

         len = min + rand(max-min+1)

         ar = []

         while count < len

            i = rand(digits_size)    # get a random array index

            random_num = digits.at(i)

            stats[random_num] += 1

            ar << random_num

            digits = digits.sort_by {rand}

            count += 1

            if count == 1 && ar.first.zero? then ar.shift; stats[0] -= 1; count = 0 end

         end

         ret << ar.to_s.to_i

      end  # iter

      #ret

      ret << stats   # optional

   end



end



puts



min = 8

max = 13

iter = 10000

result =  RNG.new.num(min, max, iter)

stats = result.pop



t = 0

stats.each_value { |v| t += v }  # get the overall frequency



n = 0

t = t * 1.0

m = t / 10.0 



stats.sort.each do |k,v| 

   i = (v / t) * 100.0

   x = v > m ? v - m : m - v

   y = (x / m) * 100.0

   n += i

   puts "#{k}  ::  #{"%.2f" % m}  ::  #{v}  ::  #{"%.2f" % i} %  ::  #{ "%.2f" % ((m-v) * -1) }  ::  #{"%.2f" % y} %"

end



puts n, m, t



puts

puts RNG.new.num(20, 25, 10)[0..-2]



puts

rn = RNG.new.num(80, 100)

puts rn[0..-2]

p rn.last

Lottery

ntk () — Thu, 18 Oct 2007 12:30:09 GMT

Author: ntk
License: The MIT License, Copyright (c) 2007 ntk
Inspired by: Lottery Quick Pick
Hints: Lottery mathematics





class Lottery



   def playing(i=1, n1=0, r1=0..0, n2=0, r2=0..0)



      return nil unless i.is_a?(Integer) && n1.is_a?(Integer) && n1 > 0 && r1.is_a?(Range) && n2.is_a?(Integer) && r2.is_a?(Range)

      stats = Hash.new(0)   # optional; stores the frequency of picked numbers, cf. stats[random_num] += 1

      ret = []



      random_procs = [ 

         lambda { |a,n| a.slice!(rand(n)) },      # a is an array, n an integer

         lambda { |a,n| a.slice!(rand(n) * -1) }, 

         lambda { |a,n| a.slice!((rand(n) * -1) + rand(n)) },

         lambda { |a,n| a.slice!((rand(n) - rand(n))) }, 

         lambda { |a,n| a.reverse!.slice!(rand(n)) },

         lambda { |a,n| a.reverse!.slice!(rand(n) * -1) }, 

         lambda { |a,n| a.reverse!.slice!((rand(n) * -1) + rand(n)) },

         lambda { |a,n| a.reverse!.slice!((rand(n) - rand(n))) } 

      ] 



      i.times do

         numbers = r1.to_a.sort_by {rand}

         num = numbers.size

         ar = []

         count = 0



         while count < n1 + n2



            count += 1



            if count > n1 && n2 > 0

               numbers2 = r2.to_a

               numbers2 = (numbers & numbers2).sort_by {rand} 

               num2 = numbers2.size

               count2 = 0

               while count2 < n2

                  count2 += 1

                  i = rand(random_procs.size)   # get a random array index for random_procs array

                  random_procs = random_procs.sort_by {rand}      

                  random_num = random_procs.at(i).call(numbers2, num2)

                  stats[random_num] += 1

                  ar << random_num

                  numbers2 = numbers2.sort_by {rand}

                  num2 -= 1

               end   # while



               if count2 > 1 then count += (count2-1) end  

               #break   # alternative

              

            else



               i = rand(random_procs.size)                               

               random_procs = random_procs.sort_by {rand}  

  

               random_num = random_procs.at(i).call(numbers, num)

               stats[random_num] += 1

               ar << random_num



=begin

              if rand(num) % (rand(num) + 1) == rand(num)               #  alternative

                  random_num = random_procs.at(i).call(numbers, num)

                  stats[random_num] += 1

                  ar << random_num

               else

                  random_num = numbers.slice!(rand(num))

                  stats[random_num] += 1

                  ar << random_num

               end

=end



               numbers = numbers.sort_by {rand}

               num -= 1



            end          # if

         end             # while

         ret << ar

      end

      #ret

      ret << stats   # optional

   end



end





puts

puts "UK National Lottery, German Lotto 6/49, ..."



ndraws = 10

result = Lottery.new.playing(ndraws, 6, 1..49)

stats = result.pop



result.each do |t| 

   puts t.sort.join('-')

end



puts

n = 0

stats.sort.each do |k,v| 

   i = (v / (ndraws * 6.0) ) * 100.0

   n += i

   puts "#{k} :: #{v} :: #{i}%"

end

p n



puts

puts "US Powerball Lottery"

Lottery.new.playing(ndraws, 5, 1..55, 1, 1..42)[0..-2].each do |t| 

   n = t.pop

   puts "#{t.sort.join('-')}  and  #{n}"

end



puts

puts "Mega Millions"

Lottery.new.playing(ndraws, 5, 1..56, 1, 1..46)[0..-2].each do |t| 

   n = t.pop

   puts "#{t.sort.join('-')}  and  #{n}"

end



puts

puts "EuroMillions"

Lottery.new.playing(ndraws, 5, 1..50, 2, 1..9)[0..-2].each do |t| 

   n1 = t.pop

   n2 = t.pop

   puts "#{t.sort.join('-')}  and  #{n1}-#{n2}"

end





puts

puts "US Powerball Lottery"



ndraws = 10

nums = [5, 23, 25, 33, 34]

pb = 31



matches = []

n = 0



Lottery.new.playing(ndraws, 5, 1..55, 1, 1..42)[0..-2].each do |t| 

   i = t.pop

   t = t.sort

   if (t & nums).size > matches.size || ( matches.size == 0 && i == pb )

      matches = t & nums 

      i == pb ? n = i : n = 0

   end

end



if pb == n && !matches.empty?

   puts "Your numbers: #{nums.join('-')} and #{pb}\nYour matches: #{matches.join('-')} and #{n}"

elsif pb == n && matches.empty?

   puts "Your numbers: #{nums.join('-')} and #{pb}\nYour matches: [] and #{n}"

else

   puts "Your numbers: #{nums.join('-')} and #{pb}\nYour matches: #{matches.join('-')}"

end

Generate prime numbers

scvalex (Alexandru Scvortov) — Fri, 22 Jun 2007 14:33:42 GMT

GenPrime(int n) generates all prime numbers smaller or equal to n, storing them in prime[]. NrPrime is the number of prime numbers calculated.

This uses the "Primes Sieve of Eratosthenes".



#include 



int prime[78498]; // 78498 is the number of prime numbers smaller than 1 000 000.

int nrPrime = 0;



void genPrime(int n) {

  char p[1000001];



  memset(p, 1, sizeof(char) * (n + 1));



  p[1] = 0;



  int i = 2;

  for (; i <= n; ++i)

    if (p[i]) {

      prime[nrPrime] = i;

      ++nrPrime;

      int j = i * 2;

      while (j <= n) {

        p[j] = 0;

        j += i;

      }

    }

}

Javascript: convert integer to word (e.g. 18 -> "eighteen")

maximile (Max Williams) — Tue, 20 Mar 2007 22:15:23 GMT

// Convert a positive integer less than 1000 to its word representation.



var units = new Array ("Zero", "One", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine", "Ten", "Eleven", "Twelve", "Thirteen", "Fourteen", "Fifteen", "Sixteen", "Seventeen", "Eighteen", "Nineteen");

var tens = new Array ("Twenty", "Thirty", "Forty", "Fifty", "Sixty", "Seventy", "Eighty", "Ninety");



function num(it) {

	var theword = "";

	var started;

	if (it>999) return "Lots";

	if (it==0) return units[0];

	for (var i = 9; i >= 1; i--){

		if (it>=i*100) {

			theword += units[i];

			started = 1;

			theword += " hundred";

			if (it!=i*100) theword += " and ";

			it -= i*100;

			i=0;

		}

	};

	

	for (var i = 9; i >= 2; i--){

		if (it>=i*10) {

			theword += (started?tens[i-2].toLowerCase():tens[i-2]);

			started = 1;

			if (it!=i*10) theword += "-";

			it -= i*10;

			i=0

		}

	};

	

	for (var i=1; i < 20; i++) {

		if (it==i) {

			theword += (started?units[i].toLowerCase():units[i]);

		}

	};

	return theword;

}

Zeropad a number

albert006 (Albert Peschar) — Mon, 01 Jan 2007 23:19:29 GMT

Pads a string with 0's until the given length limit is met



function zeropad($number, $limit) {

  return (strlen($number) >= $limit) ? $number : zeropad("0" . $number, $limit);

}

?>

Bash: checking if a variable is a number

Grzechu (Grzegorz Rumatowski) — Wed, 13 Dec 2006 18:54:32 GMT

Shellscript checking if a variable is a number



#!/bin/bash

read VARIABLE

if [ $VARIABLE -eq $VARIABLE 2> /dev/null ]; then

echo $VARIABLE is a number

else

echo $VARIABLE is not a number

fi

BigNumber //JavaScript Class

jonasraoni (Jonas Raoni Soares Silva) — Thu, 15 Jun 2006 20:19:36 GMT

Offers a extremely high precision level to make mathematical operations. For integers there is no limits and for floating point numbers, the class allows setting the maximum precision.

[UPDATED CODE AND HELP CAN BE FOUND HERE]

The class always returns new instances of BigNumber on the operations, so to the set value of a object, use the "set" method.



//+ Jonas Raoni Soares Silva

//@ http://jsfromhell.com/classes/bignumber [v1.0]



BigNumber = function(n, p, r){

    var o = this, i;

    if(n instanceof BigNumber){

        for(i in {precision: 0, roundType: 0, _sign: 0, _dec: 0}) o[i] = n[i];

        o._buffer = n._buffer.slice();

        return;

    }

    o.precision = isNaN(p = Math.abs(p)) ? BigNumber.defaultPrecision : p;

    o.roundType = isNaN(r = Math.abs(r)) ? BigNumber.defaultRoundType : r;

    o._sign = (n += "").charAt(0) == "-";

    o._dec = ((n = n.replace(/[^\d.]/g, "").split(".", 2))[0] = n[0].replace(/^0+/, "") || "0").length;

    for(i = (n = o._buffer = (n.join("") || "0").split("")).length; i; n[--i] = +n[i]);

    o.round();

};

with({$: BigNumber, o: BigNumber.prototype}){

    $.ROUND_HALF_EVEN = ($.ROUND_HALF_DOWN = ($.ROUND_HALF_UP = ($.ROUND_FLOOR = ($.ROUND_CEIL = ($.ROUND_DOWN = ($.ROUND_UP = 0) + 1) + 1) + 1) + 1) + 1) + 1;

    $.defaultPrecision = 40;

    $.defaultRoundType = $.ROUND_HALF_UP;

    o.add = function(n){

        if(this._sign != (n = new BigNumber(n))._sign) return n._sign = !n._sign, this.subtract(n);

        var o = new BigNumber(this), a = o._buffer, b = n._buffer, la = o._dec,

        lb = n._dec, n = Math.max(la, lb), i, r;

        la != lb && ((lb = la - lb) > 0 ? o._zeroes(b, lb, 1) : o._zeroes(a, -lb, 1));

        i = (la = a.length) == (lb = b.length) ? a.length : ((lb = la - lb) > 0 ? o._zeroes(b, lb) : o._zeroes(a, -lb)).length;

        for(r = 0; i; r = (a[--i] = a[i] + b[i] + r) / 10 >>> 0, a[i] %= 10);

        return r && ++n && a.unshift(r), o._dec = n, o.round();

    };

    o.subtract = function(n){

        if(this._sign != (n = new BigNumber(n))._sign) return n._sign = !n._sign, this.add(n);

        var o = new BigNumber(this), x = o.compare(n) + 1, a = x ? o : n, b = x ? n : o, la = a._dec, lb = b._dec, n = la, i, j;

        a = a._buffer, b = b._buffer, la != lb && ((lb = la - lb) > 0 ? o._zeroes(b, lb, 1) : o._zeroes(a, -lb, 1));

        for(i = (la = a.length) == (lb = b.length) ? a.length : ((lb = la - lb) > 0 ? o._zeroes(b, lb) : o._zeroes(a, -lb)).length; i;){

            if(a[--i] < b[i]){

                for(j = i; j && !a[--j]; a[j] = 9);

                --a[j], a[i] += 10;

            }

            b[i] = a[i] - b[i];

        }

        return o._sign = !x, o._dec = n, o._buffer = b, o.round();

    };

    o.multiply = function(n){

        var o = new BigNumber(this), r = o.compare(n = new BigNumber(n)) + 1, a = (r ? o : n)._buffer,

        b = (r ? n : o)._buffer, la = a.length, lb = b.length, x = new BigNumber, i, j, s;

        for(i = lb; i; x.set(x.add(new BigNumber(s.join("")))))

            for(s = (new Array(lb - --i)).join("0").split(""), r = 0, j = la; j;

            r += a[--j] * b[i], s.unshift(r % 10), r = r / 10 >>> 0);

        return r && x._buffer.unshift(r), o._dec = (o._buffer = x._buffer).length - la - lb + o._dec + n._dec, o.round();

    };

    o.divide = function(n){

        if((n = new BigNumber(n)) == "0") throw new Error("division by 0");

        var o = new BigNumber(this), a = o._buffer, b = n._buffer, la = a.length - o._dec,

        lb = b.length - n._dec, s = a._sign != b._sign, dec = 0, buffer = [], x = new BigNumber, y, i;

        la != lb && ((lb = la - lb) > 0 ? o._zeroes(b, lb) : o._zeroes(a, -lb));

        o._dec = a.length, n._dec = b.length;

        b = n, a._sign = false, b._sign = false;

        while(o != 0 && (buffer.length - dec) <= o.precision){

            x.set(0);

            for(i = 0; o.compare(y = x.add(b)) + 1 && ++i; x.set(y));

            if(!i){

                do ++i, o._buffer.push(0), ++o._dec;

                while(o.compare(b) < 0);

                !dec && (!buffer.length && buffer.push(0), dec = buffer.length);

                o._zeroes(buffer, --i);

                continue;

            }

            o.set(o.subtract(x)), buffer.push(i);

        }

        return o._buffer = buffer, o._dec = dec ? dec : buffer.length, o.round();

    };

    o.mod = function(n){

        return this.subtract(this.divide(n).intPart().multiply(n));

    };

    o.pow = function(n){

        var o = new BigNumber(this), i;

        if(n == 0) return o.set(1);

        for(i = Math.abs(n); --i; o.set(o.multiply(this)));

        return n < 0 ? o.set((new BigNumber(1)).divide(o)) : o;

    };

    o.set = function(n){

        return this.constructor(n), this;

    };

    o.compare = function(n){

        var a = this, la = this._dec, b = n, lb = n._dec, i, l;

        if(la != lb) return la > lb ? 1 : -1;

        for(la = (a = a._buffer).length, lb = (b = b._buffer).length, i = -1, l = Math.min(la, lb); ++i < l;)

            if(a[i] != b[i]) return a[i] > b[i] ? 1 : -1;

        return la != lb ? la > lb ? 1 : -1 : 0;

    }

    o.negate = function(){

        var n = new BigNumber(this); return n._sign ^= 1, n;

    };

    o.abs = function(){

        var n = new BigNumber(this); return n._sign = 0, n;

    };

    o.intPart = function(){

        return new BigNumber((this._sign ? "-" : "") + (this._buffer.slice(0, this._dec).join("") || "0"));

    };

    o.valueOf = o.toString = function(){

        var o = this;

        return (o._sign ? "-" : "") + (o._buffer.slice(0, o._dec).join("") || "0") + (o._dec != o._buffer.length ? "." + o._buffer.slice(o._dec).join("") : "");

    };

    o._zeroes = function(n, l, t){

        var s = ["push", "unshift"][t || 0];

        for(++l; --l;  n[s](0));

        return n;

    };

    o.round = function(){

        if("_rounding" in this) return this;

        var $ = BigNumber, r = this.roundType, b = this._buffer, d, p, n, x;

        for(this._rounding = true; this._dec > 1 && !b[0]; --this._dec, b.shift());

        for(d = this._dec, p = this.precision + d, n = b[p]; b.length > d && !b[b.length -1]; b.pop());

        x = (this._sign ? "-" : "") + (p - d ? "0." + this._zeroes([], p - d - 1).join("") : "") + 1;

        if(b.length > p){

            n && (r == $.DOWN ? false : r == $.UP ? true : r == $.CEIL ? !this._sign

            : r == $.FLOOR ? this._sign : r == $.HALF_UP ? n >= 5 : r == $.HALF_DOWN ? n > 5

            : r == $.HALF_EVEN ? n >= 5 && b[p - 1] & 1 : false) && this.add(x);

            b.splice(p, b.length - p);

        }

        return delete this._rounding, this;

    };

}

Usage