# OEIS A144311 Generating function

I'm looking for a way to use calculate OEIS A144311 efficiently in Mathematica. First, let's define the series.

In one sense or another, this series considers the number between "relative" twin primes. What do I mean by this?

Well, the number $77$ is relatively prime to {$2,3,5$} even tho it's not actually prime. We're looking for relative twin primes to this set tho: consecutive odds relatively prime to {$2,3,5$}. For instance, $(77,79)$ would be an example of such. Neither number is divisible by {$2,3,5$}.

$78$ then is the number between these. The next smallest number between a similarly relative twin prime pair would be $90$. $(89,91)$ is the next relative (and real) twin prime pair. That gives us a relative twin prime gap of 12 (being $90-78$), and 11 fails between successes. The question is, what is the largest possible twin prime gap relative to just {$2,3,5$}? This so happens to be it. 11 is the max fails for {$2,3,5$}. The series says that 29 is the max fails for {$2,3,5,7$}, and 41 for {$2,3,5,7,11$}, the first 5 primes.

Now that we've well defined these, what would be the best way to calculate them? Please, as well, show your answer calculating the 17th number in this series with the time for the calculation. All methods I've used to try to calculate these have resulted in numbers beyond what can be handled well time/memory wise even by the 10th number in the series.

• Were there a "Generating function", it would be in the OEIS entry. Find one, you're probably up for something like an Abel prize...
– ciao
May 10 '16 at 22:11

This function calculates the maximum fails for the given list upto n. It also shows the numbers which give these gaps.

 twinPrime[n_] :=
Module[{list = {2, 3, 5, 7, 11}, set, twinSet, max, res},
set = DeleteCases[Table[If[Count[Divisible[i, list], False] == Length[list],
Sow[i]], {i, 1, n}], Null];
twinSet =
Mean /@ Select[Partition[set, 2, 1], #[[2]] - #[[1]] == 2 &];
max = Max@Differences[twinSet];
res = Select[Partition[twinSet, 2, 1], #[[2]] - #[[1]] == max &];
{max - 1, res}
]


For list = {2, 3, 5} we get

twinPrime[10^2] // AbsoluteTiming

(*{0.00182891, {11, {{18, 30}, {30, 42}, {48, 60}, {60, 72}, {78, 90}}}}*)

First@twinPrime[10^5] // AbsoluteTiming

(*{1.49442, 11}*)


For list = {2, 3, 5,7} we get

 First@twinPrime[10^5] // AbsoluteTiming

(*{1.47472, 29}*)

First@twinPrime[10^6] // AbsoluteTiming

(*{16.0307, 29}*)


For list = {2, 3, 5,7,11} we get

 First@twinPrime[10^6] // AbsoluteTiming

(*{15.5588, 41}*)


For list = Table[Prime[i], {i, 1, 16}] we get

 points = Table[First@twinPrime[10^n] // AbsoluteTiming, {n, 2, 7}]

(*{{0.00368573, 11}, {0.0307609, 149}, {0.246445, 197}, {2.41628,
203}, {25.5115, 317}, {262.341, 317}}*)


Here the timing is increasing by an order of magnitude and the convergence should be checked for higher n values.

## Version 2 : Slightly improved efficiency

twinPrime2[n_, m_] :=
Module[{list = Table[Prime[x], {x, 2, m}], set, twinSet, max, res},
set = Select[Range[1, n, 2], Count[Mod[#, list], 0] == 0 &];
twinSet =
Mean /@ Select[Partition[set, 2, 1], #[[2]] - #[[1]] == 2 &];
max = Max@Differences[twinSet];
res = Select[Partition[twinSet, 2, 1], #[[2]] - #[[1]] == max &];
{max - 1}]

First@twinPrime2[10^6, 3] // AbsoluteTiming

(*{6.77699, 11}*)

First@twinPrime2[10^6, 4] // AbsoluteTiming

(*{7.06732, 29}*)

First@twinPrime2[10^6, 5] // AbsoluteTiming

(*{7.37997, 41}*)

First@twinPrime2[10^6, 16] // AbsoluteTiming

(*{11.1849, 317}*)

• How about "list = Table[Prime[x], {x,1,16}]" ? Do you get 869 as expected? May 10 '16 at 6:39
• @Elem-Teach-w-Bach-n-Math-Ed See the edit. the run time is increasing by a factor of 10. I did it up-to n=7. However it did show the same answer (317) for n=6,7 but that's not enough. u should check it for more values of n. May 10 '16 at 6:57
• I just saw that the answer should be 869, which means u will surely have to run it for higher values but its time consuming. Maybe someone else will find a more efficient way to do this. May 10 '16 at 7:02
• Yeah, this still seems more efficient than my original go at it, but not quite enough. May 10 '16 at 7:03
• We can improve the efficiency by just checking for odd numbers instead of going over all integers up to n May 10 '16 at 7:35

I was able to speed up @Hubble07's code by defining

RelativePrimes[n_Integer, p_List] :=
Complement[Range[1, n - 1, 2], Apply[Sequence, Map[Range[#, n - 1, #] &, p]]]


to find the values in set. The new function becomes

twinPrime3[n_, m_] :=
Module[{set},
set = RelativePrimes[n, Prime[Range[2, m]]];
Max[Differences[Pick[Rest[set], Differences[set], 2] - 1]] - 1]


Tests of twinPrime3 with n=10^6 and m equal from 2 to 9 are about 100 times faster than twinPrime2; however, the answer is incorrect for m greater than 9. Even with twinPrime3[5*10^8,10], the answer is still incorrect.

• Since the product of the first 10 primes is $6469693230$ or about $6.5*10^9$, that may be the reason. The context value must exist in that order of magnitude's worth of difference. May 10 '16 at 23:22
• Here's what I mean code-wise: RelativePrimes[n_Integer, p_List] := Complement[Range[1, n - 1, 2], Apply[Sequence, Map[Range[#, n - 1, #] &, p]]] twinPrime3[m_] := Module[{set, plist, mp}, plist = Table[Prime[n], {n, 1, m}]; mp = Apply[Times, plist]; set = RelativePrimes[mp, Prime[Range[2, m]]]; Max[Differences[Pick[Rest[set], Differences[set], 2] - 1]] - 1] AbsoluteTiming[twinPrime3[8]] May 11 '16 at 4:39
• This replaces $n$ with the product of the first $m$ primes. Problem is, on my little raspberry pi, I get a memory error running $m=9$ or above. Perhaps it's a memory error that's causing your issues, Kenny? May 11 '16 at 4:43
• Hey, was wondering if you had any thought regarding my answer down below. I used some of your thinking in mine, but I'm a novice and not sure if there are any major improvements to make. May 12 '16 at 18:17
• Sorry to keep bothering you. Found out that my code is slower than yours by about a factor of 8. You think we could similarly adapt your code as I adapted mine? Since I didn't write it, I'm not sure all it does. Thanks! May 14 '16 at 22:59

Ok, I got one now. Took some of the ideas of other posts, but then calculated not based on the number between the twins of the twin prime pair, but based on the number $2$ less than the first. I first create a set $\equiv 3 \pmod 6$ as only these work, and I no longer need to check $2$ or $3$. Then I create a massive set of numbers 2 or 4 mod the set of primes 5 and above being tested. I complement the former with the latter, find max of differences, then voila!

tpgaps[x_] :=
Module[{plist, tplist},
plist = Times @@ Table[Prime[y], {y, 1, x}];
tplist =
Complement[Range[3, plist, 6],
Flatten[Table[
Union[Range[2, plist, Prime[n]], Range[4, plist, Prime[n]]], {n,
3, x}]]];
AbsoluteTiming[Max[Differences[tplist]] - 1]]
tpgaps[8]

(*{0.078764, 149}*)


This is very speedy, but extremely memory intensive, and just like the codes above, can't really do $9$ or above on a raspberry pi. =(

My thought is, what if we could split plist into overlapping chunks? Do $(0,10000)$, then $(5000,15000)$, $(10000,20000)$ etc., find the gaps there, then take the max of those? Anyone willing to help me adjust my code to do so?

EDIT:

Ok, so I spent much time putting this together, but finally got something that works. It still needs testing and probably some performance tuning, but here it is:

tpgaps[x_, b_] := Module[{plist, m, tplist},
plist = Times @@ Table[Prime[y], {y, 1, x}];
m = plist/b;
tplist =
Max[ParallelTable[
Max[Differences[
Complement[Range[3 + m*z/2, m (1 + z/2), 6],
Flatten[Table[Union[
Range[2 + z/2*Prime[n]*Floor[m/Prime[n]], m (1 + z/2),
Prime[n]],
Range[4 + z/2*Prime[n]*Floor[m/Prime[n]], m (1 + z/2),
Prime[n]]],
{n, 3, x}]]]]], {z, 0, 2*plist/m - 2}]] - 1;
{tplist}]

TableForm[Table[Flatten[myAbsoluteTiming[tpgaps[8, i]],{i,1,101,10}]]


$b$ is the number of pieces to split the operation into.

2ND EDIT:

After some personal performance tuning and troubleshooting, I've got a formula that runs fast, is accurate, and doesn't have memory issues. Here it is:

tpgaps1000[x_] := Module[{plist, m, c, tplist},
plist = Times @@ Table[Prime[y], {y, 1, x}];
m = 2*Floor[(plist/(200*Times @@ Table[Prime[y], {y, 10, x}]))/2];
c = z*Prime[n]*Floor[m/(2*Prime[n])];
tplist = Max[Table[Max[Differences[Complement[
Range[3 + m*z/2, m (1 + z/2), 6], Flatten[Table[Union[
Range[2 + c, m (1 + z/2), Prime[n]],
Range[4 + c, m (1 + z/2), Prime[n]]],
{n, 3, x}]]]]], {z, 0, 2*plist/m - 2}]];
{tplist - 1}]
AbsoluteTiming[TableForm[Table[Flatten[AbsoluteTiming[tpgaps1000[i]]],{i,6,9}]]]
{81.3416,
0.098396   65
0.196206  107
3.67567   149
77.3713   203
}


3rd Edit

With some help here I now have this:

list3[a_, b_, c_] := Union[
Range[ChineseRemainder[{3, a}, {6, Prime[c]}], b, 6*Prime[c]],
Range[ChineseRemainder[{3, a + 2}, {6, Prime[c]}], b, 6*Prime[c]]]

tpgaps7[x_, z_]:=Module[{plist, tplist, tuples, k, mgap},
plist = Times @@ Table[Prime[y], {y, 1, x}];
tplist =
Complement[Range[3, plist, 6],
Flatten[Table[
Union[Range[2, plist, Prime[n]], Range[4, plist, Prime[n]]], {n,
3, x}]]];
tuples=Tuples[Table[Table[{a, plist, c}, {a, 0, Prime[c - 1]}], {c, x + 1, z}]];
k = Length[tuples];
mgap = AbsoluteTiming[Max[Table[
Max[
Differences@
Complement[tplist,
Flatten[Union@@ Apply[list3, Part[tuples, j], {1}]]]], {j,1, k}]] - 1];
{mgap}]

Do[Print[tpgaps7[i, i]], {i, 3, 7}]
Do[Print[tpgaps7[7, i]], {i, 8, 10}]

(*{{0.00009,11}}
{{0.000019,29}}
{{0.000035,41}}
{{0.000186,65}}
{{0.003238,107}}
{{0.016698,149}}
{{0.501871,203}}
{{15.6921,257}}*)