View Single Post
 2020-02-12, 19:21 #81 mart_r     Dec 2008 you know...around... 2×11×29 Posts I finally got around comparing the speed of Thomas' and my code on my only PC that's able to run both programs. My former code was actually slower than I thought, but the current one below is almost on par, and may be faster on other systems. In the range 6e15+, its throughput is about 780e6/s on my workplace PC. Maybe going to run it for this search, after some minor other projects... Code: 'please enter your values of choice here: k1=range start, k2=range end, mg&=mingap report '---------------------------------------------------------------------------------------- k1 = 9E+16 k2 = 9.00005E+16 mg& = 200000 '---------------------------------------------------------------------------------------- tm = Timer Cells(1, 4) = "initializing..." Dim p&(22914109), c&(700244), d&(17506124), o&(22914109, 4), e&(22914109), a%(215656441), aa%(6), ab%(76), ac%(1000) 'p: primes 3 thru 29#/15, c: temporary variable for the calculation of d, d: index of admissible quads in Z/29#Z, o: offsets for sieve, e: change in offsets per 29#, a: 0/1=prime/composite, aa/ab/ac: 0/1=skip sieve in sieving routine below n/y w = Int(k1 / 6469693230#) x = 1 + Int(k2 / 6469693230#) b& = x - w: 'length of search interval 'calculate primes p(1)=3, p(2)=5, ... p(22914109)=431312881 - for k2>186e15 the number of primes has to be increased! Cells(1, 4) = "calculate primes below 29#/15... ( 0%)" For q& = 3 To 20759 Step 2 i% = 1 For t& = 3 To Sqr(q&) Step 2 If q& Mod t& = 0 Then i% = 0: t& = Sqr(q&) Next If i% Then j& = j& + 1: p&(j&) = q&: e&(j&) = q& - 215656441 Mod q& Next Cells(1, 4) = "calculate primes below 29#/15... ( 1%)" q2& = 1 For q& = 1 To 2338 m& = p&(q&) m& = (m& - 1) / 2 For t& = m& + p&(q&) To 215656441 Step p&(q&) a%(t&) = 1 Next If q& = q2& Then q2& = q2& * 4: t1$= Str$(Int(Sqr(q&) / 0.8)): Cells(1, 4) = "calculate primes below 29#/15... (" + t1$+ "%)" Next Cells(1, 4) = "calculate primes below 29#/15... ( 60%)" For t& = 10385 To 215656441 If a%(t&) = 0 Then j& = j& + 1: p&(j&) = 2 * t& + 1: e&(j&) = p&(j&) - 215656441 Mod p&(j&) Next t& = 0 'calculate indexes of admissible quads in Z/29#Z where d=p/30 corresponds to a quadruplet p+{11,13,17,19} (mod 29#) Cells(1, 4) = "calculate admissible quads in Z/29#Z..." d&(0) = 0 d&(1) = 3 d&(2) = 6 dq& = 1 jp& = 30 For n% = 4 To 8 Step 2 q& = -1 jp& = jp& * p&(n% - 1) dq& = dq& * (p&(n% - 1) - 4) For j& = 0 To p&(n%) - 1 For k& = 0 To dq& - 1 r& = jp& Mod p&(n%) r& = (r& * j& + 30 * d&(k&) + 19) Mod p&(n%) If r& <> 0 And r& <> 2 And r& <> 6 And r& <> 8 Then q& = q& + 1: c&(q&) = jp& / 30 * j& + d&(k&) Next Next q& = -1 jp& = jp& * p&(n%) dq& = dq& * (p&(n%) - 4) For j& = 0 To p&(n% + 1) - 1 For k& = 0 To dq& - 1 r& = jp& Mod p&(n% + 1) r& = (r& * j& + 30 * c&(k&) + 19) Mod p&(n% + 1) If r& <> 0 And r& <> 2 And r& <> 6 And r& <> 8 Then q& = q& + 1: d&(q&) = jp& / 30 * j& + c&(k&) Next Next Next 'estimate time for calculating the offset values l& = Int(Sqr(x * 6469693230#)): 'upper trial division limit t0 = Int((Timer - tm) * l& / 800000000 + 0.5) * 10 If t0 < 90 Then t1$ = " seconds)" Else t0 = Int(t0 / 60 + 0.5): t1$= " minutes)" 'initialize offset values for sieve in a, where a(1,2,3...n)={0 or 1, quadruplet or at least one composite} correspond to the quad p+{11,13,17,19} where p=30*d(n) (mod 29#) Cells(1, 4) = "calculate offsets for sieve... (this may take about" + Str$(t0) + t1$j& = 10: 'start at p(10)=31, since all smaller factors are taken care of with the d array Do y = Int(((w - 1) * 510510 / p&(j&) - Int((w - 1) * 510510 / p&(j&))) * p&(j&) + 0.5) y = Int(((y * 12673 + 29) / p&(j&) - Int((y * 12673 + 29) / p&(j&))) * p&(j&) + 0.5) r& = y: 'the calculation of r=((w-1)*29#+29) mod p is split up because I only have 53 bits of precision r& = r& + p&(j&) * (r& Mod 2) pp = p&(j&) For s& = 2 To 30 Step 2 If ((p&(j&) Mod 30) * s& - r& + 18) Mod 30 = 0 Then o&(j&, 1) = 1 + (pp * s& - r& + 18) / 30 If ((p&(j&) Mod 30) * s& - r& + 16) Mod 30 = 0 Then o&(j&, 2) = 1 + (pp * s& - r& + 16) / 30 If ((p&(j&) Mod 30) * s& - r& + 12) Mod 30 = 0 Then o&(j&, 3) = 1 + (pp * s& - r& + 12) / 30 If ((p&(j&) Mod 30) * s& - r& + 10) Mod 30 = 0 Then o&(j&, 4) = 1 + (pp * s& - r& + 10) / 30 Next j& = j& + 1 Loop While p&(j&) < l& 'these are used to speed up the sieving routine below For q& = 0 To 180 Step 30 If q& Mod 7 > 2 Then aa%(q& / 30) = 1 Next For q& = 0 To 2280 Step 30 If q& Mod 7 > 2 Or q& Mod 11 = 0 Or q& Mod 11 = 2 Or q& Mod 11 = 6 Or q& Mod 11 = 8 Then ab%(q& / 30) = 1 Next For q& = 0 To 30000 Step 30 If q& Mod 7 > 2 Or q& Mod 11 = 0 Or q& Mod 11 = 2 Or q& Mod 11 = 6 Or q& Mod 11 = 8 Or q& Mod 13 = 11 Or q& Mod 13 = 0 Or q& Mod 13 = 4 Or q& Mod 13 = 6 Then ac%(q& / 30) = 1 Next sk& = mg& * 2 / 25 'this is for the number of values in d skipped after a quadruplet is found to accelerate the search 'a small overview of the max. difference between sk consecutive values of d vs. the largest possible gap that would be missed: '1000/12782, 2000/25155, 3000/37552, 4000/49958, 5000/62374, 6000/74788, 7000/87080, 8000/99404, 9000/111696, 10000/123987, 15000/185683, 20000/247370, 25000/308830 (ratio approaching 12.3189 = 29#/30/(3*7*9*13*15*19*25)) 'overhead finished, start search in intervals of 6469693230=29# Do Cells(1, 4) = "searching... (k =" + Str$((w + f&) * 6469693230#) + " at a rate of" + Str$(Int(6469.69323 / (Timer - tm))) + "e6 per second)" tm = Timer Erase a%() For m& = 10 To 2047 For nn% = 1 To 4 o&(m&, nn%) = (o&(m&, nn%) + e&(m&)) Mod p&(m&) v& = o&(m&, nn%) For n% = 0 To 1000 If ac%(v& Mod 1001) Then GoTo 1: 'skip sieve when 7, 11 or 13 divides one member of 30d+{11,13,17,19} For k& = v& To 215656441 Step p&(m&) * 1001 a%(k&) = 1 Next 1 v& = v& + p&(m&) Next Next Next For m& = 2048 To 49151 For nn% = 1 To 4 o&(m&, nn%) = (o&(m&, nn%) + e&(m&)) Mod p&(m&) v& = o&(m&, nn%) For n% = 0 To 76 If ab%(v& Mod 77) Then GoTo 2: 'skip sieve when 7 or 11 divides one member of 30d+{11,13,17,19} For k& = v& To 215656441 Step p&(m&) * 77 a%(k&) = 1 Next 2 v& = v& + p&(m&) Next Next Next For m& = 49152 To 786431 For nn% = 1 To 4 o&(m&, nn%) = (o&(m&, nn%) + e&(m&)) Mod p&(m&) v& = o&(m&, nn%) For n% = 0 To 6 If aa%(v& Mod 7) Then GoTo 3: 'skip sieve when 7 divides one member of 30d+{11,13,17,19} For k& = v& To 215656441 Step p&(m&) * 7 a%(k&) = 1 Next 3 v& = v& + p&(m&) Next Next Next For m& = 786432 To j& - 1: 'for larger values p(m) the skip (mod 7) as above has no longer an effect w.r.t. speed and can even be counterproductive For nn% = 1 To 4 o&(m&, nn%) = (o&(m&, nn%) + e&(m&)) Mod p&(m&) For k& = o&(m&, nn%) To 215656441 Step p&(m&) a%(k&) = 1 Next Next Next 'sieving done, now looking for gaps For m& = 0 To 17506124 If a%(d&(m&) + 1) = 0 Then g& = d&(m&) - u&: If g& >= mg& Then GoSub 5 Else i% = 0: GoSub 6: 'see below Next f& = f& + 1 If f& Mod 256 = 0 Then ActiveWorkbook.Save: 'saves every 920e6/(throughput in k's per second) minutes Loop While f& <= b& Cells(1, 4) = "finished search in the interval [" + Str$(k1) + "; " + Str$(k2) + " ]" End 'what to do after a quad gap is found 5 If i% = 0 Then GoSub 7: GoTo 6: 'i=0: no jump in m, so no interval left unchecked 'when i=1, there's an unchecked interval t+[1..sk] after a quad 30*d(t)+{11,13,17,19} that has to be examined i% = 0 For q& = t& + sk& To t& + 1 Step -1 If a%(d&(q&) + 1) = 0 Then u& = d&(q&): q& = t& + 1 Next g& = d&(m&) - u& If g& >= mg& Then GoSub 7 6 u& = d&(m&) If m& < 17301504 Then t& = m&: m& = m& + sk&: i% = 1: 'skip some values (sk=mg/12.5, see above) after a quadruplet is found to accelerate the search, but only if an interval overlap is out of question Return 'output gap size and initial quadruplet member p where p=(w+f)*29#+30*(d-g)+11 is transformed to a string, else it would show as "1.2345...E+15" and a digit or two would be missed 7 h& = h& + 1: 'line number in the spreadsheet Cells(h&, 1) = g& k0 = w + f& kk$ = Str$(Int(((k0 * 215656441 + d&(m&) - g&) * 30 + 11) / 10000)): 'the four rightmost digits are separated 'error checking routine b/c of possible rounding errors in k0 (compare 5th digit from the right of output value p): kc$ = Str$(Int(((Int((k0 / 100000 - Int(k0 / 100000)) * 100000 + 0.5) * 56441 + d&(m&) - g&) * 30 + 11) / 10000)) If Right$(kc$, 1) <> Right$(kk$, 1) Then kk$ = Str$(Val(kk$) + 1) If Right$(kc$, 1) <> Right$(kk$, 1) Then kk$= Str$(Val(kk$) - 2): 'if there's a difference of 1 then it could be either way If Right$(kc$, 1) <> Right$(kk$, 1) Then Cells(h&, 3) = "roundoff error - actual value is ± 10000*n": 'according to my calculations, this should never happen, but please notify me if it does k0 = (Int((k0 / 10000 - Int(k0 / 10000)) * 10000 + 0.5) * 6441 + d&(m&) - g&) * 30 + 11 kk$ = kk$+ Right$(Str$(k0), 4): 'attach the four rightmost digits of p Cells(h&, 2) = "'" + kk$ Return Attached Thumbnails   Last fiddled with by mart_r on 2020-02-12 at 19:48