View Single Post
Old 2020-02-12, 19:21   #81
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

2×11×29 Posts
Default

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
Click image for larger version

Name:	quadgaps_output.JPG
Views:	74
Size:	52.2 KB
ID:	21767  

Last fiddled with by mart_r on 2020-02-12 at 19:48
mart_r is offline   Reply With Quote