# How can I calculate OEIS A144311 efficiently?

I'm looking for a way to calculate OEIS A144311 efficiently.

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? Also, what is the 17th number in this series? All methods I've used to try to calculate these have resulted in numbers beyond what can be handled well even in something like Mathematica even by the 10th number in the series.

• +1, nice question. I removed the generating-functions tag. When in doubt, please avail yourself of the tag summaries in selecting tags. – joriki May 10 '16 at 7:18

I see basically two approaches, and a hybrid of them that's more efficient than either.

First, you could try out all combinations of residues other than $-1,1$ that the primes could have at the one end of the sequence, and then see how far you get in each case. The number of combinations would be

$$\prod_{i=3}^n(p_i-2)\;,$$

i.e. $3\cdot5\cdot9\cdot11\cdot15\cdot17\cdot21\cdot27\cdot29\cdot35\cdot39\cdot41\cdot45\cdot51\cdot57\approx4.6\cdot10^{19}$ to reach $p_{17}$, so this is not feasible.

Second, you could start out without fixing residues for the primes and start moving to the right, and whenever you hit a slot that's not covered yet, you pick a prime whose residue hasn't been fixed yet and try the two possible residues to cover the slot. That would require trying out $2^nn!$ different choices, or $2^{17}\cdot17!\approx4.7\cdot10^{19}$ to reach $p_{17}$, funny enough almost the same number, so this isn't feasible, either.

But note that the first approach is costly for large primes and the second approach is costly for many primes, so we can advantageously combine them by adopting the first approach for the small primes and the second approach for the large primes. That is, we try out all combinations of residues for the primes up to, say, $p_7=17$, and then we move right and fix one of the remaining residues whenever we hit an uncovered slot. To reach $p_{17}$, this requires trying out

$$3\cdot5\cdot9\cdot11\cdot15\cdot2^{10}\cdot10!\approx8\cdot10^{13}$$

combinations, a whole lot better and within reach of our electronic friends.

Here's java code that applies this approach. I used it to replicate the sequence up to $p_{14}$, which took $1000$ seconds on my MacBook Pro, so you can expect to find the result for $p_{17}$ in about $5760000$ seconds, or about two months. It's probably not a coincidence that this is on the boundary of what's feasible with today's computers; the calculation of the values up to $p_{16}$ (which would take about $2$ days this way) dates from $2009$, and I suspect that they were calculated in precisely this manner.

• @Elem-Teach-w-Bach-n-Math-Ed: I just figured out how to make SageMathCloud files public, so here it is directly: cloud.sagemath.com/projects/…. You can make a copy of it there to run it yourself. – joriki May 18 '16 at 6:45
• @Elem-Teach-w-Bach-n-Math-Ed: I don't mean finding all such numbers in the form of a list -- that list couldn't be held in memory. I mean iterating over them and in each case finding the distance to the next one. Both approaches are based on choosing residues at the left-hand "failure" of a sequence and traversing the sequence to find the distance to the right-hand "failure"; the difference is just that in the first approach all residues are chosen at the beginning, whereas in the second approach they're chosen as you go along the sequence and need new primes to defer the right-hand "failure". – joriki May 24 '16 at 0:07
• @Elem-Teach-w-Bach-n-Math-Ed: It sounds like you got most of it right now, but two things that may be causing the problem at the end: a) I don't think of it as starting at $30$. In the hybrid, you start at $x$, like in the second approach, but you fix the residues of $x$ up to $p_7=17$ up front. So if you like, you can think of it as starting at $30+x$, where $x$ is some multiple of $\prod_{i=1}^7p_i$; but I just think of it as $x$ with $7$ residues fixed. – joriki May 24 '16 at 6:26
• @Elem-Teach-w-Bach-n-Math-Ed: b) It sounds as if you're trying the larger primes in a fixed order, from $p_10$ to $p_8$. The idea is to try them in all possible orders, with all possible choices of $\pm1$. When you've tried something, say $+1$ for $p_{10}$, then $+1$ for $p_9$, then $+1$ for $p_8$, you keep going along the sequence until you hit the next failure; since you have no primes left to cover it, you can't go any further, so you update the maximum, backtrack and try $-1$ for $p_8$ at the point where you tried $+1$ before. Once you're done trying that, you backtrack even further ... – joriki May 24 '16 at 6:32
• @Elem-Teach-w-Bach-n-Math-Ed: ... to the point where you tried $+1$ for $p_9$ and try $-1$ for $p_9$ instead at that point. That leads to some more recursion for the two choices $\pm1$ for $p_8$ again, and then when you backtrack after that you now try $p_8$ before $p_9$, again all $4$ combinations of $\pm1$. The code uses recursion to make the backtracking easy; at each decision, you just make one recursive call for each of the options, and when the call returns you're back where you were, ready to try the next option. ("$p_10$" above was meant to be "$p_{10}$"; too late to edit now. :-) – joriki May 24 '16 at 6:36

This can be considered as a kind of covering problem, which can be formulated as an integer linear programming problem. Suppose you want to see if $A144311(n) \ge m$.
We ask to cover the set $S = \{1, \ldots, m\}$, where each prime $p$ of the first $n$ primes will cover those members of $S$ congruent to $a_p \pm 1 \mod p$ for some $a_p \in \{0, \ldots, p-1\}$. Given such a covering, if $T$ is an integer $\equiv -a_p \mod p$ for each of these $p$, then $T+1, \ldots, T+m$ is a sequence of $m$ consecutive integers that satisfies the requirement (the existence of $T$ is guaranteed by the Chinese Remainder Theorem).

In the integer linear programming formulation, we have $0-1$ variables $x_{p,a}$ for each of our primes $p$ and $a = 0 \ldots p-1$, with the constraints

\eqalign{\sum_p \sum_{a \equiv i \pm 1 \mod p} x_{p,a} &\ge 1 \ \text{for}\ i=1 \ldots m\cr \sum_{a=0}^{p-1} x_{p,a} &= 1 \ \text{for each}\ p \cr \text{all} & x_{p,a} \in \{0,1\}}

For $n = 12$, it took Cplex on my computer about 14 seconds to show that $m=527$ is feasible, and about 22 seconds to show that $m=528$ is infeasible. However, I expect computing times to rise rather rapidly as $n$ increases, although maybe not as rapidly as for exhaustive search methods.

Another possibility is to attack this with a SAT solver.

Feasible solutions might be found using heuristic search methods such as tabu search or simulated annealing, providing lower bounds. Proving upper bounds is likely to be more difficult.

• So I'd imagine n=16 m=870 is beyond reasonable time limits? – Elem-Teach-w-Bach-n-Math-Ed May 26 '16 at 2:52
• Interesting -- we should compare times to see which method scales better. Times in seconds (on my MacBook Pro) for my solution are: $1,6,70,1000$ for $11\le n\le14$; they increase roughly as expected. Can you go up to $n=14$ or $15$ to compare? – joriki May 26 '16 at 5:01
• For $n=14$, it took 56.8 seconds to show $m=617$ is feasible, but 892.5 seconds to show $m=618$ is infeasible. – Robert Israel May 27 '16 at 6:34
• That gives a factor of $\frac{1000\cdot22}{892.5\cdot6}\approx4$ better scaling from $n=12$ to $n=14$ than my method; so you might be able to reach $n=17$ in a slightly more reasonable time (also @Elem-Teach-w-Bach-n-Math-Ed). – joriki May 27 '16 at 9:29
• By the way, a better command for "mod" used as a binary operator is \bmod, which doesn't produce the big space that's meant for the use of \mod at the end of equations. – joriki May 27 '16 at 9:31

Taking @joriki's approach, I've now written a version that runs in Mathematica. I'm including now only the latest and just previous revisions to keep this from getting lengthy. The last revision saw a 5x speedup do to the number between twin primes being divisible by 6, and (in the case of my code) the number $2$ and $4$ more than the twins respectively being $\equiv 3 \pmod 6$. Here's the final:

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[(*Parallel*)Table[Max[Differences@Complement[
tplist, Flatten[Union @@ Apply[list3, Part[tuples, j], {1}]]]],
{j, 1, k}]] - 1];
{mgap}]

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

(*Times when running in Parallel on my Rasp Pi*)
{{0.064885,7,107}}
{{0.097869,8,149}}
{{1.26251,9,203}}
{{35.9214,10,257}}

(*Times for non-parallel online programming lab*)
{{0.003366,7,107}}
{{0.015974,8,149}}
{{0.498724,9,203}}
{{15.4054,10,257}}


Previous DEBUGGED/OPTIMIZED VERSION

list3[a_,b_,c_]:=Union[Range[a,b,Prime[c]],Range[a+2,b,Prime[c]]]

tpgaps6[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[ParallelTable[
Max[Differences@Complement[tplist,
Flatten[Union@@Apply[list3,Part[tuples,j],{1}]]]],{j,1,k}]]-1];
{z,mgap}]

Do[Print[AbsoluteTiming[tpgaps6[7,i]]],{i,7,10}]

(*Times when running in Parallel on my Rasp Pi*)
{0.554687,{7,{0.053122,107}}}
{0.659064,{8,{0.164598,149}}}
{5.49219,{9,{5.00954,203}}}
{173.238,{10,{172.752,257}}}

(*Times for non-parallel online programming lab*)
{0.078032,{7,{0.00339,107}}}
{0.105818,{8,{0.035588,149}}}
{1.45627,{9,{1.38502,203}}}
{48.1555,{10,{48.0839,257}}}

• But $227$ and $317$ are both off by $30$? – joriki May 26 '16 at 8:02
• @joriki Yeah, I did catch that as well, but was headed to bed. I think I may know why. Need to do some testing to make sure tho. – Elem-Teach-w-Bach-n-Math-Ed May 26 '16 at 14:08
• @joriki Found the error... Part of the code read "c" instead of Prime[c]. In essence, when what was needed was a list of numbers from 1 to prime 10 (29), instead it was just calculating based off just 1-10... thus although the error is fixed, everything now takes about 3x as long. =( – Elem-Teach-w-Bach-n-Math-Ed May 27 '16 at 0:56
• @joriki Found about a 5x speedup on my code! Will be posting shortly! – Elem-Teach-w-Bach-n-Math-Ed May 28 '16 at 1:03