Sieve of Zakiya -- Number Theory Sieves (NTS)
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
1 2 #!/usr/local/bin/ruby -w 3 4 require 'benchmark' 5 6 class Integer 7 8 def primes1 9 # classical brute-force Sieve of Eratosthenes (SoE) 10 # initialze sieve array with all integers starting at i=2 11 12 sieve = [nil, nil].concat((2..self).to_a) 13 14 (2 .. Math.sqrt(self).to_i).each do |i| 15 16 next unless sieve[i] 17 (i*i).step(self, i) { |j| sieve[j] = nil } 18 end 19 sieve.compact! 20 end 21 22 def primes2 23 # classical brute-force Sieve of Eratosthenes (SoE) 24 # initialize sieve array with odd integers for odd indexes, then 2 25 26 sieve = []; 3.step(self, 2) { |i| sieve[i] = i }; sieve[2] = 2 27 28 3.step(Math.sqrt(self).to_i, 2) do |i| 29 next unless sieve[i] 30 (i*i).step(self, i<<1) { |j| sieve[j] = nil } 31 end 32 sieve.compact! 33 end 34 35 def primes3 36 # http://everything2.com/index.pl?node_id=1176369 37 # all prime candidates > 3 are of form 6*k+1 and 6*k+5 38 # initialize sieve array with only these candidate values 39 n1, n2 = 1, 5; sieve = [] 40 while n2 < self 41 n1 += 6; n2 += 6; sieve[n1] = n1; sieve[n2] = n2 42 end 43 # now initialize sieve array with first 3 primes, resize array 44 sieve[2]=2; sieve[3]=3; sieve[5]=5; sieve=sieve[0..self] 45 46 5.step(Math.sqrt(self).to_i, 2) do |i| 47 next unless sieve[i] 48 (i*i).step(self, i<<1) {|j| sieve[j] = nil } 49 end 50 sieve.compact! 51 end 52 53 def primesP3 54 # all prime candidates > 3 are of form 6*k+1 and 6*k+5 55 # initialize sieve array with only these candidate values 56 n1, n2 = 1, 5; sieve = [] 57 while n2 < self 58 n1 += 6; n2 += 6; sieve[n1] = n1; sieve[n2] = n2 59 end 60 # now initialize sieve array with primes < 6, resize array 61 sieve[2]=2; sieve[3]=3; sieve[5]=5; sieve=sieve[0..self] 62 63 5.step(Math.sqrt(self).to_i, 2) do |i| 64 next unless sieve[i] 65 # p5= 5*i, k = 6*i, p7 = 7*i 66 p1 = 5*i; k = p1+i; p2 = k+i 67 while p1 <= self 68 sieve[p1] = nil; sieve[p2] = nil; p1 += k; p2 += k 69 end 70 end 71 sieve.compact! 72 end 73 74 75 def primesP3a 76 # all prime candidates > 3 are of form 6*k+1 and 6*k+5 77 # initialize sieve array with only these candidate values 78 # where sieve contains the odd integers representatives 79 # convert integers to array indices/vals by i = (n-3)>>1 = (n>>1)-1 80 n1, n2 = -1, 1; lndx= (self-1) >>1; sieve = [] 81 while n2 < lndx 82 n1 +=3; n2 += 3; sieve[n1] = n1; sieve[n2] = n2 83 end 84 #now initialize sieve array with (odd) primes < 6, resize array 85 sieve[0] =0; sieve[1]=1; sieve=sieve[0..lndx-1] 86 87 5.step(Math.sqrt(self).to_i, 2) do |i| 88 next unless sieve[(i>>1) - 1] 89 # p5= 5*i, k = 6*i, p7 = 7*i 90 # p1 = (5*i-3)>>1; p2 = (7*i-3)>>1; k = (6*i)>>1 91 i6 = 6*i; p1 = (i6-i-3)>>1; p2 = (i6+i-3)>>1; k = i6>>1 92 while p1 < lndx 93 sieve[p1] = nil; sieve[p2] = nil; p1 += k; p2 += k 94 end 95 end 96 return [2] if self < 3 97 [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 } 98 end 99 100 def primesP5 101 # all prime candidates > 5 are of form 30*k+(1,7,11,13,17,19,23,29) 102 # initialize sieve array with only these candidate values 103 n1, n2, n3, n4, n5, n6, n7, n8 = 1, 7, 11, 13, 17, 19, 23, 29; sieve = [] 104 while n8 < self 105 n1 +=30; n2 += 30; n3 += 30; n4 += 30 106 n5 +=30; n6 += 30; n7 += 30; n8 += 30 107 sieve[n1] = n1; sieve[n2] = n2; sieve[n3] = n3; sieve[n4] = n4 108 sieve[n5] = n5; sieve[n6] = n6; sieve[n7] = n7; sieve[n8] = n8 109 end 110 # now initialize sieve with the primes < 30, resize array 111 sieve[2]=2; sieve[3]=3; sieve[5]=5; sieve[7]=7; sieve[11]=11; sieve[13]=13 112 sieve[17]=17; sieve[19]=19; sieve[23]=23; sieve[29]=29; sieve=sieve[0..self] 113 114 7.step(Math.sqrt(self).to_i, 2) do |i| 115 next unless sieve[i] 116 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i, k=30*i 117 n6 = 6*i; p1 = i+n6; p2 = n6+n6-i; p3 = p1+n6; p4 = p2+n6 118 p5 = p3+n6; p6 = p4+n6; p7 = p6+n6; k = p7+i; p8 = k+i 119 while p1 <= self 120 sieve[p1] = nil; sieve[p2] = nil; sieve[p3] = nil; sieve[p4] = nil 121 sieve[p5] = nil; sieve[p6] = nil; sieve[p7] = nil; sieve[p8] = nil 122 p1 += k; p2 += k; p3 += k; p4 += k; p5 += k; p6 += k; p7 += k; p8 += k 123 end 124 end 125 sieve.compact! 126 end 127 128 def primesP5a 129 # all prime candidates > 5 are of form 30*k+(1,7,11,13,17,19,23,29) 130 # initialize sieve array with only these candidate values 131 # where sieve contains the odd integers representatives 132 # convert integers to array indices/vals by i = (n-3)>>1 = n>>1 - 1 133 n1, n2, n3, n4, n5, n6, n7, n8 = -1, 2, 4, 5, 7, 8, 10, 13 134 lndx= (self-1) >>1; sieve = [] 135 while n8 < lndx 136 n1 +=15; n2 += 15; n3 += 15; n4 += 15 137 n5 +=15; n6 += 15; n7 += 15; n8 += 15 138 sieve[n1] = n1; sieve[n2] = n2; sieve[n3] = n3; sieve[n4] = n4 139 sieve[n5] = n5; sieve[n6] = n6; sieve[n7] = n7; sieve[n8] = n8 140 end 141 # now initialize sieve with the (odd) primes < 30, resize array 142 sieve[0]=0; sieve[1]=1; sieve[2]=2; sieve[4]=4; sieve[5]=5 143 sieve[7]=7; sieve[8]=8; sieve[10]=10; sieve[13]=13 144 sieve = sieve[0..lndx-1] 145 146 7.step(Math.sqrt(self).to_i, 2) do |i| 147 next unless sieve[(i>>1) - 1] 148 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i, k=30*i 149 # maps to: p1= (7*i-3)>>1, --- p8=(31*i-3)>>1, k=30*i>>1 150 n2=i<<1; n4=i<<2; n8=i<<3; n16=i<<4; n32=i<<5; n12=n8+n4 151 p1 = (n8-i-3)>>1; p2 = (n12-i-3)>>1; p3 = (n12+i-3)>>1; p4 = (n16+i-3)>>1 152 p5 = (n16+n4-i-3)>>1; p6 = (n16+n8-i-3)>>1; p7 = (n32-n4+i-3)>>1 153 p8 = (n32-i-3)>>1; k = (n32-n2)>>1 154 while p1 < lndx 155 sieve[p1] = nil; sieve[p2] = nil; sieve[p3] = nil; sieve[p4] = nil 156 sieve[p5] = nil; sieve[p6] = nil; sieve[p7] = nil; sieve[p8] = nil 157 p1 += k; p2 += k; p3 += k; p4 += k; p5 += k; p6 += k; p7 += k; p8 += k 158 end 159 end 160 return [2] if self < 3 161 [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 } 162 end 163 164 def primesP7 165 # all prime candidates > 7 are of form 210*k+(1,11,13,17,19,23,29,31,37 166 # 41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121,127,131 167 # 137,139,143,149,151,157,163,167,169,173,179,181,187,191,193,197,199,209) 168 # initialize sieve array with only these candidate values 169 n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 = 1, 11, 13, 17, 19, 23, 29, 31, 37, 41 170 n11, n12, n13, n14, n15, n16, n17, n18 = 43, 47, 53, 59, 61, 67, 71, 73 171 n19, n20, n21, n22, n23, n24, n25, n26 = 79, 83, 89, 97, 101, 103, 107, 109 172 n27, n28, n29, n30, n31, n32, n33, n34 = 113, 127, 131, 137, 139, 149, 151, 157 173 n35, n36, n37, n38, n39, n40, n41, n42, n43 = 163, 167, 173, 179, 181, 191, 193, 197, 199 174 n44, n45, n46, n47, n48 = 121, 143, 169, 187, 209 175 num = self - self[0]^1; sieve = [] 176 while n48 < num 177 n1 +=210; n2 += 210; n3 += 210; n4 += 210; n5 += 210 178 n6 +=210; n7 += 210; n8 += 210; n9 += 210; n10 += 210 179 n11 +=210; n12 += 210; n13 += 210; n14 += 210; n15 += 210 180 n16 +=210; n17 += 210; n18 += 210; n19 += 210; n20 += 210 181 n21 +=210; n22 += 210; n23 += 210; n24 += 210; n25 += 210 182 n26 +=210; n27 += 210; n28 += 210; n29 += 210; n30 += 210 183 n31 +=210; n32 += 210; n33 += 210; n34 += 210; n35 += 210 184 n36 +=210; n37 += 210; n38 += 210; n39 += 210; n40 += 210 185 n41 +=210; n42 += 210; n43 += 210; n44 += 210; n45 += 210 186 n46 += 210; n47 += 210; n48 += 210 187 sieve[n1] = n1; sieve[n2] = n2; sieve[n3] = n3; sieve[n4] = n4 188 sieve[n5] = n5; sieve[n6] = n6; sieve[n7] = n7; sieve[n8] = n8 189 sieve[n9] = n9; sieve[n10] = n10; sieve[n11] = n11; sieve[n12] = n12 190 sieve[n13] = n13; sieve[n14] = n14; sieve[n15] = n15; sieve[n16] = n16 191 sieve[n17] = n17; sieve[n18] = n18; sieve[n19] = n19; sieve[n20] = n20 192 sieve[n21] = n21; sieve[n22] = n22; sieve[n23] = n23; sieve[n24] = n24 193 sieve[n25] = n25; sieve[n26] = n26; sieve[n27] = n27; sieve[n28] = n28 194 sieve[n29] = n29; sieve[n30] = n30; sieve[n31] = n31; sieve[n32] = n32 195 sieve[n33] = n33; sieve[n34] = n34; sieve[n35] = n35; sieve[n36] = n36 196 sieve[n37] = n37; sieve[n38] = n38; sieve[n39] = n39; sieve[n40] = n40 197 sieve[n41] = n41; sieve[n42] = n42; sieve[n43] = n43; sieve[n44] = n44 198 sieve[n45] = n44; sieve[n46] = n46; sieve[n47] = n47; sieve[n48] = n48 199 end 200 # now initialize sieve with the (odd) primes < 210, resize array 201 sieve[2]=2; sieve[3]=3; sieve[5]=5; sieve[7]=7; sieve[11]=11; sieve[13]=13 202 sieve[17]=17; sieve[19]=19; sieve[23]=23; sieve[29]=29; sieve[31]=31 203 sieve[37]=37; sieve[41]=41; sieve[43]=43; sieve[47]=47; sieve[53]=53 204 sieve[59]=59; sieve[61]=61; sieve[67]=67; sieve[71]=71; sieve[73]=73 205 sieve[79]=79; sieve[83]=83; sieve[89]=89; sieve[97]=97; sieve[101]=101 206 sieve[103]=103; sieve[107]=107; sieve[109]=109; sieve[113]=113; sieve[127]=127 207 sieve[131]=131; sieve[137]=137; sieve[139]=139; sieve[149]=149; sieve[151]=151 208 sieve[157]=157; sieve[163]=163; sieve[167]=167; sieve[173]=173; sieve[179]=179 209 sieve[181]=181; sieve[191]=191; sieve[193]=193; sieve[197]=197; sieve[199]=199 210 sieve = sieve[0..self] 211 11.step(Math.sqrt(self).to_i, 2) do |i| 212 next unless sieve[i] 213 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i, k=30*i 214 # maps to: p1= (7*i-3)>>1, --- p8=(31*i)>>1, k=30*i>>1 215 #n6=6*i; n7=n6+i; n11=n6+n6-i; n13=n6+n7; n17=n11+n6 216 #n19=n13+n6; n23=n17+n6; n29=n23+n6; n31=n29+(i<<1) 217 #n4=i<<2; n8=i<<3; n16=i<<4; n7=n8-i; n11=n7+n4; n13=n11+(i<<1) 218 #n17=n16+i; n19=n11+n8; n23=n16+n7; n29=n16+n13; n31=(i<<5)-i 219 p1 = (11*i); p2 = (13*i); p3 = (17*i); p4 = (19*i) 220 p5 = (23*i); p6 = (29*i); p7 = (31*i); p8 = (37*i) 221 p9 = (41*i); p10 = (43*i); p11 = (47*i); p12 = (53*i) 222 p13 = (59*i); p14 = (61*i); p15 = (67*i); p16 = (71*i) 223 p17 = (73*i); p18 = (79*i); p19 = (83*i); p20 = (89*i) 224 p21 = (97*i); p22 = (101*i); p23 = (103*i); p24 = (107*i) 225 p25 = (109*i); p26 = (113*i); p27 = (127*i); p28 = (131*i) 226 p29 = (137*i); p30 = (139*i); p31 = (149*i); p32 = (151*i) 227 p33 = (157*i); p34 = (163*i); p35 = (167*i); p36 = (173*i) 228 p37 = (179*i); p38 = (181*i); p39 = (191*i); p40 = (193*i) 229 p41 = (197*i); p42 = (199*i); p43 = (211*i) 230 p44 = (121*i); p45 = 143*i; p46 = 169*i; p47 = 187*i; p48 = 209*i 231 k = (210*i) 232 while p1 <= num 233 sieve[p1] = nil; sieve[p2] = nil; sieve[p3] = nil; sieve[p4] = nil 234 sieve[p5] = nil; sieve[p6] = nil; sieve[p7] = nil; sieve[p8] = nil 235 sieve[p9] = nil; sieve[p10] = nil; sieve[p11] = nil; sieve[p12] = nil 236 sieve[p13] = nil; sieve[p14] = nil; sieve[p15] = nil; sieve[p16] = nil 237 sieve[p17] = nil; sieve[p18] = nil; sieve[p19] = nil; sieve[p20] = nil 238 sieve[p21] = nil; sieve[p22] = nil; sieve[p23] = nil; sieve[p24] = nil 239 sieve[p25] = nil; sieve[p26] = nil; sieve[p27] = nil; sieve[p28] = nil 240 sieve[p29] = nil; sieve[p30] = nil; sieve[p31] = nil; sieve[p32] = nil 241 sieve[p33] = nil; sieve[p34] = nil; sieve[p35] = nil; sieve[p36] = nil 242 sieve[p37] = nil; sieve[p38] = nil; sieve[p39] = nil; sieve[p40] = nil 243 sieve[p41] = nil; sieve[p42] = nil; sieve[p43] = nil; sieve[p44] = nil 244 sieve[p45] = nil; sieve[p46] = nil; sieve[p47] = nil; sieve[p48] = nil 245 p1 += k; p2 += k; p3 += k; p4 += k; p5 +=k; p6 += k; p7 += k 246 p8 += k; p9 += k; p10 +=k; p11 += k; p12 +=k; p13 +=k; p14 +=k 247 p15 +=k; p16 +=k; p17 +=k; p18 += k; p19 +=k; p20 +=k; p21 +=k 248 p22 +=k; p23 +=k; p24 +=k; p25 += k; p26 +=k; p27 +=k; p28 +=k 249 p29 +=k; p30 +=k; p31 +=k; p32 += k; p33 +=k; p34 +=k; p35 +=k 250 p36 +=k; p37 +=k; p38 +=k; p39 += k; p40 +=k; p41 +=k; p42 +=k 251 p43 +=k; p44 +=k; p45 +=k; p46 +=k; p47 +=k; p48 +=k 252 end 253 end 254 sieve.compact! 255 end 256 257 def primesP7a 258 # all prime candidates > 7 are of form 210*k+(1,11,13,17,19,23,29,31,37 259 # 41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121,127,131 260 # 137,139,143,149,151,157,163,167,169173,179,181,187,191,193,197,199,209) 261 # initialize sieve array with only these candidate values 262 # where sieve contains the odd integers representatives 263 # convert integers to array indices/vals by i = (n )>>1 = (n>>1)-1 264 n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 = -1, 4, 5, 7, 8, 10, 13, 14, 17, 19 265 n11, n12, n13, n14, n15, n16, n17, n18 = 20, 22, 25, 28, 29, 32, 34, 35 266 n19, n20, n21, n22, n23, n24, n25, n26 = 38, 40, 43, 47, 49, 50, 52, 53 267 n27, n28, n29, n30, n31, n32, n33, n34 = 55, 62, 64, 67, 68, 73, 74, 77 268 n35, n36, n37, n38, n39, n40, n41, n42, n43 = 80, 82, 85, 88, 89, 94, 95, 97, 98 269 n44, n45, n46, n47, n48 = 59, 70, 83, 92, 103 270 lndx= (self-1)>>1; sieve = [] 271 while n48 < lndx 272 n1 +=105; n2 += 105; n3 += 105; n4 += 105; n5 += 105 273 n6 +=105; n7 += 105; n8 += 105; n9 += 105; n10 += 105 274 n11 +=105; n12 += 105; n13 += 105; n14 += 105; n15 += 105 275 n16 +=105; n17 += 105; n18 += 105; n19 += 105; n20 += 105 276 n21 +=105; n22 += 105; n23 += 105; n24 += 105; n25 += 105 277 n26 +=105; n27 += 105; n28 += 105; n29 += 105; n30 += 105 278 n31 +=105; n32 += 105; n33 += 105; n34 += 105; n35 += 105 279 n36 +=105; n37 += 105; n38 += 105; n39 += 105; n40 += 105 280 n41 +=105; n42 += 105; n43 += 105; n44 += 105; n45 += 105 281 n46 +=105; n47 += 105; n48 += 105 282 sieve[n1] = n1; sieve[n2] = n2; sieve[n3] = n3; sieve[n4] = n4 283 sieve[n5] = n5; sieve[n6] = n6; sieve[n7] = n7; sieve[n8] = n8 284 sieve[n9] = n9; sieve[n10] = n10; sieve[n11] = n11; sieve[n12] = n12 285 sieve[n13] = n13; sieve[n14] = n14; sieve[n15] = n15; sieve[n16] = n16 286 sieve[n17] = n17; sieve[n18] = n18; sieve[n19] = n19; sieve[n20] = n20 287 sieve[n21] = n21; sieve[n22] = n22; sieve[n23] = n23; sieve[n24] = n24 288 sieve[n25] = n25; sieve[n26] = n26; sieve[n27] = n27; sieve[n28] = n28 289 sieve[n29] = n29; sieve[n30] = n30; sieve[n31] = n31; sieve[n32] = n32 290 sieve[n33] = n33; sieve[n34] = n34; sieve[n35] = n35; sieve[n36] = n36 291 sieve[n37] = n37; sieve[n38] = n38; sieve[n39] = n39; sieve[n40] = n40 292 sieve[n41] = n41; sieve[n42] = n42; sieve[n43] = n43; sieve[n44] = n44 293 sieve[n45] = n45; sieve[n46] = n46; sieve[n47] = n47; sieve[n48] = n48 294 end 295 # now initialize sieve with the (odd) primes < 210, resize array 296 sieve[0]=0; sieve[1]=1; sieve[2]=2; sieve[4]=4; sieve[5]=5 297 sieve[7]=7; sieve[8]=8; sieve[10]=10; sieve[13]=13; sieve[14]=14 298 sieve[17]=17; sieve[19]=19; sieve[20]=20; sieve[22]=22; sieve[25]=25 299 sieve[28]=28; sieve[29]=29; sieve[32]=32; sieve[34]=34; sieve[35]=35 300 sieve[38]=38; sieve[40]=40; sieve[43]=43; sieve[47]=47; sieve[49]=49 301 sieve[50]=50; sieve[52]=52; sieve[53]=53; sieve[55]=55; sieve[62]=62 302 sieve[64]=64; sieve[67]=67; sieve[68]=68; sieve[73]=73; sieve[74]=74 303 sieve[77]=77; sieve[80]=80; sieve[82]=82; sieve[85]=85; sieve[88]=88 304 sieve[89]=89; sieve[94]=94; sieve[95]=95; sieve[97]=97; sieve[98]=98; 305 sieve = sieve[0..lndx-1] 306 307 11.step(Math.sqrt(self).to_i, 2) do |i| 308 next unless sieve[(i>>1) - 1] 309 # p1=7*i, p2=11*i, p3=13*i ---- p6=23*i, p7=29*i, p8=31*i, k=30*i 310 # maps to: p1= (7*i-3)>>1, --- p8=(31*i)>>1, k=30*i>>1 311 #n6=6*i; n7=n6+i; n11=n6+n6-i; n13=n6+n7; n17=n11+n6 312 #n19=n13+n6; n23=n17+n6; n29=n23+n6; n31=n29+(i<<1) 313 #n4=i<<2; n8=i<<3; n16=i<<4; n7=n8-i; n11=n7+n4; n13=n11+(i<<1) 314 #n17=n16+i; n19=n11+n8; n23=n16+n7; n29=n16+n13; n31=(i<<5)-i 315 p1 = (11*i-3)>>1; p2 = (13*i-3)>>1; p3 = (17*i-3)>>1; p4 = (19*i-3)>>1 316 p5 = (23*i-3)>>1; p6 = (29*i-3)>>1; p7 = (31*i-3)>>1; p8 = (37*i-3)>>1 317 p9 = (41*i-3)>>1; p10 = (43*i-3)>>1; p11 = (47*i-3)>>1; p12 = (53*i-3)>>1 318 p13 = (59*i-3)>>1; p14 = (61*i-3)>>1; p15 = (67*i-3)>>1; p16 = (71*i-3)>>1 319 p17 = (73*i-3)>>1; p18 = (79*i-3)>>1; p19 = (83*i-3)>>1; p20 = (89*i-3)>>1 320 p21 = (97*i-3)>>1; p22 = (101*i-3)>>1; p23 = (103*i-3)>>1; p24 = (107*i-3)>>1 321 p25 = (109*i-3)>>1; p26 = (113*i-3)>>1; p27 = (127*i-3)>>1; p28 = (131*i-3)>>1 322 p29 = (137*i-3)>>1; p30 = (139*i-3)>>1; p31 = (149*i-3)