QUADRUPLET PRIMES AND THEIR GAPS - Page 8

mart_r · 2020-02-09, 11:36

Quote:

Originally Posted by robert44444uk

Do bother, mart_r. What you are doing is progress.

Personally I dont have much time to devote to learning/ running new software at present. Others may have more time.

I only needed some extra space for sulking. My room wasn't big enough.

If you can't use the program right now, that's perfectly okay. All I needed was a tiny indication of how useful my code is, so thanks for the reply

Thomas11 · 2020-02-09, 12:08

Thanks for sharing your VBA code, mart_r!

I will have a look into it and give it a try when I'm back in my office next week (don't have Excel at home)...

Bobby Jacobs · 2020-02-09, 22:35

Quote:

Originally Posted by mart_r

First holes, in this case, are the gaps for which you have to search farthest to find them.

For ordinary prime gaps, this is A100180

I call these late prime gaps. I wonder if there are any patterns with the late gaps for prime quadruplets. The first few late gaps are 1, 4, 6, 20, 125, ...

mart_r · 2020-02-12, 19:21

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

Thomas11 · 2020-02-13, 16:45

Meanwhile I've found the time to test your VBA code on my workplace PC.
While I'm getting only 550e6/s (indicating that my PC is slower than yours), I can confirm that the two programs are comparable in speed.
Moreover, since both programs yield identical results it gives us further confidence that our independent and different algorithms and implementations are both working correctly.

mart_r · 2020-02-13, 17:53

Quote:

Originally Posted by Thomas11

Meanwhile I've found the time to test your VBA code on my workplace PC.
While I'm getting only 550e6/s (indicating that my PC is slower than yours), I can confirm that the two programs are comparable in speed.
Moreover, since both programs yield identical results it gives us further confidence that our independent and different algorithms and implementations are both working correctly.

Nice. Thanks!

Uncwilly · 2020-02-13, 18:54

Quote:

Originally Posted by mart_r

I finally got around comparing the speed of Thomas' and my code on my only PC that's able to run both programs.

That crashes my version of Excel (2010, corporate license issues).

Thomas11 · 2020-02-14, 08:52

Quote:

Originally Posted by Uncwilly

That crashes my version of Excel (2010, corporate license issues).

Well, occasionally a had the same effect here (Excel 2016, 32 bit):
The screen freezes and it tells me "the application doesn't respond anymore".
But the computations continue and once they are finished, the screen is refreshed and the results are shown.

When this happens, one can still minimize the application window, do some other things, and check once in a while whether the computations are finished...

Uncwilly · 2020-02-14, 20:35

Quote:

Originally Posted by Thomas11

Well, occasionally a had the same effect here (Excel 2016, 32 bit):.

Excel dies and exits.

mart_r · 2020-02-15, 13:51

Quote:

Originally Posted by Uncwilly

Excel dies and exits.

That's odd. I haven't encountered that behaviour with either Excel 2003, 2010, or 2016. Could you elaborate on the corporate license issues? That also seems new to me, even after more than 20 years of working with Excel.

Something about the non-responding window though, sometimes when the window doesn't reappear after I click on it on the task bar, I click on and off it while pressing Esc until the VBA window appears where you have the options to stop or continue the calculation.

Uncwilly · 2020-02-15, 22:49

Well, I got a new machine, it runs Win 10. But, rather than running a new version of Excel it runs the older version. Either SysAdmin does not like 365 for security or license issues.

When I tried your code, I opened Excel, Alt-F11, pasted code in and hit the run button. Poof, excel flat out is gone. Can't test it again for a few days. I won't be at the machine for a while.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
GAPS BETWEEN PRIME PAIRS (Twin Primes)	rudy235	Prime Gap Searches	227	2021-01-21 21:59
Gaps between maximal prime gaps	Bobby Jacobs	Prime Gap Searches	52	2020-08-22 15:20
Prime Quadruplet Emirps	Trejack	Puzzles	11	2016-04-22 19:27
Gaps of Primes?	PawnProver44	Miscellaneous Math	10	2016-04-10 19:32
Top-5 prime quadruplet found	gd_barnes	Riesel Prime Search	10	2007-10-27 01:54

2020-02-09, 12:08	#79
Thomas11 Feb 2003 3564₈ Posts	Thanks for sharing your VBA code, mart_r! I will have a look into it and give it a try when I'm back in my office next week (don't have Excel at home)...

2020-02-13, 16:45	#82
Thomas11 Feb 2003 2²×3²×53 Posts	Meanwhile I've found the time to test your VBA code on my workplace PC. While I'm getting only 550e6/s (indicating that my PC is slower than yours), I can confirm that the two programs are comparable in speed. Moreover, since both programs yield identical results it gives us further confidence that our independent and different algorithms and implementations are both working correctly.

2020-02-15, 22:49	#88
Uncwilly 6809 > 6502 """"""""""""""""""" Aug 2003 101×103 Posts 250A₁₆ Posts	Well, I got a new machine, it runs Win 10. But, rather than running a new version of Excel it runs the older version. Either SysAdmin does not like 365 for security or license issues. When I tried your code, I opened Excel, Alt-F11, pasted code in and hit the run button. Poof, excel flat out is gone. Can't test it again for a few days. I won't be at the machine for a while.