1
$\begingroup$

The prime sequence starts:

s = {  11,   41,   61,   83,  113,  101,  151,  181,  233,  223,  263,
      293,  353,  383,  419,  401,  479,  467,  541, 1009,  599,  631,
      661,  691,  727,  751,  787,  797,  809,  877,  907,  919,  967,
      991, 9001, 1031}. 

The digit-decomposed list starts:

t = {1, 1, 4, 1, 6, 1, 8, 3, 1, 1, 3, 1, 0, 1, 1, 5, 1, 1, 8, 1, 2, 3, 
     3, 2, 2, 3, 2, 6, 3, 2, 9, 3, 3, 5, 3, 3, 8, 3, 4, 1, 9, 4, 0, 1,
     4, 7, 9, 4, 6, 7, 5, 4, 1, 1, 0, 0, 9, 5, 9, 9, 6, 3, 1, 6, 6, 1,
     6, 9, 1, 7, 2, 7, 7, 5, 1, 7, 8, 7, 7, 9, 7, 8, 0, 9, 8, 7, 7, 9, 
     0, 7, 9, 1, 9, 9, 6, 7, 9, 9, 1, 9, 0, 0, 1, 1, 0, 3, 1}.

The initial $11$ in $s$ is given. Thereafter, each new term must be the smallest possible prime not already used.

In other words, we want 11 to say: "In position 1 is a 1," 41 to say: "In position 4 is a 1," 61 to say: "In position 6 is a 1," 83 to say: "In position 8 is a 3," 113 to say: "In position 11 is a 3," 101 to say: "In position 10 is a 1," and so on. All of our statements must be truthful.

As per the $s$ and $t$ above, one may verify that the following yields all True:

Table[
  Floor[s[[i]] / 10] > Length[t] || t[[ Floor[s[[i]] / 10] ]] == Mod[s[[i]], 10], 
  {i, Length[s]}
]

Some of our primes (in the above example 9001) may refer to future positions which we have yet to calculate, which is why I had to Or the Floor[s[[i]] / 10] > Length[t] in the code. The existence of closely-spaced future digits leads to situations where, potentially, very large primes are required to comply with all of our demands.

For example, term #1446 is 190901, taking up positions 7006-7011. Positions 7012-7020 and 7022-7024 are already assigned with digits: 191737191?371. That would make term #1447 19173719153.

For the background origin of this challenge, see here. I've struggled with some very inefficient code to generate a few thousand terms over the course of many days. I'm hoping someone can do better.

I will add here my just revamped code to show that even I can kludge together a working program that is significantly faster than my original attempt. Unfortunately it introduces seemingly arbitrary parameters that — although they speed up the calculation — are grounded in poorly understood (by me) logical concepts underlying the sequence:

s = {11}; t = {1, 1}; u = {t}; v = {1}; p = 13; While[Length[s] < 101, 
 While[f = Floor[p/10]; m = Mod[p, 10]; w = Join[t, IntegerDigits[p]]; 
  a = Length[t] - 32; b = Length[w]; c = Select[v, a < # <= b &]; 
  If[c != {}, 
   q = Table[
     w[[v[[k]]]] == u[[k]][[-1]], {k, Position[v, c[[1]]][[1, 1]], 
      Position[v, c[[-1]]][[1, 1]]}]]; 
  MemberQ[v, f] || (f <= Length[w] && w[[f]] != m) || Union[q] != {True}, 
  p = NextPrime[p]]; AppendTo[s, p]; t = w; u = Union[u, {{f, m}}]; 
 v = Union[v, {f}]; p = Max[13, NextPrime[9*a]]]; s

Additionally, the program will not intelligently handle the occurrence of the very large primes that appear at #1447 and #3868.

$\endgroup$
  • $\begingroup$ Is this a puzzle or something? What is the context of this? At a minimum, share your current code. $\endgroup$ – MarcoB Mar 9 at 2:44
  • $\begingroup$ More of a challenge than a puzzle, although I am specifically intrigued by how large a number will eventually appear. I've added links to the background origin, as well as my current code. $\endgroup$ – Hans Havermann Mar 9 at 9:11
1
$\begingroup$

I present a re-coded solution involving two functions. I have avoided using AppendTo and procedural constructs where possible.

  • update generates a new addition to an existing sequence; it can be applied repeatedly in a Nest construct to grow the sequence, avoiding having to calculate known values from scratch;
  • validQ checks whether the sequence is still valid after the addition of a new number; it is used as a termination requirement in the update function.

Here are their implementations:

ClearAll[validQ]
validQ[s_] := Module[{t, tLen, validPairs},
  t = Flatten@IntegerDigits@s;
  tLen = Length[t];
  validPairs = Select[Transpose@{Floor[s/10], Mod[s, 10]}, #[[1]] <= tLen &];
  And @@ (t[[#1]] == #2 & @@@ validPairs)
]

ClearAll[update]
update[{s_, t_}] := Module[{candidate, sNew},
  candidate = 41;
  While[
    sNew = Join[s, {candidate}];
    MemberQ[s, candidate] || ! validQ[sNew],
    candidate = NextPrime[candidate]
  ];
  {sNew, Join[t, IntegerDigits@candidate]}
]

You can then calculate new elements using update repeatedly on an existing sequence. For instance, the following generates 100 new elements, starting from the first given one:

Nest[update, {{11}, {1, 1}}, 100]

(*Out:
{{  11,   41,  61,    83,  113,  101,  151,  181,  233,  223,
   263,  293,  353,  383,  419,  401,  479,  467,  541, 1009,
   599,  631,  661,  691,  727,  751,  787,  797,  809,  877,
   907,  919,  967, 991,  9001, 1031, 1063, 1151, 1171, 1187,
  1201, 1237, 1303, 1321, 1361, 1373, 1453, 1481, 1597, 1601,
  1693, 2309, 1669, 1721, 1759, 1777, 1801, 1867, 1877, 1987, 
  1997, 2069, 2131, 2143, 2237, 2273, 2393, 2267, 2293, 2357, 
  2389, 2437, 2503, 2591, 2633, 2671, 2711, 2753, 2777, 2797, 
  2837, 2917, 2939, 2999, 3037, 3023, 3083, 3163, 3191, 3203, 
  3271, 3313, 3301, 3323, 3373, 3413, 3433, 3461, 3593, 3631, 
  3623}, 
{<omitted list of list of digits>} 
*)

Your implementation takes roughly 7s to calculate the first 100 elements of the sequence; this one takes roughly 4.7s to do the same. It's not a huge difference, but it's something.


One can "checkpoint" these time-consuming calculations in the notebook by subdividing them in chunks that return intermediate results, which could be saved in named variable or exported for safe-keeping. For instance, the following steps add a total of 60 new elements, calculating and reporting them 20 at a time.

Nest[update, {{11}, {1, 1}}, 20]
Nest[update, %, 20]
Nest[update, %, 20]

A comparison with the results obtained "from scratch" shows that, of course, the two are equivalent:

% == Nest[update, {{11}, {1, 1}}, 60]     (* Out: True *)
| improve this answer | |
$\endgroup$
1
$\begingroup$

I'm not fluent in Mathematica (and don't have access to it so I don't know how my code compares to those given here) but I added a link to an efficient PARI/GP program in oeis.org/A333085, which computes 10^4 terms in a few seconds.

As I saw it, the main bottleneck were the large jumps required when the first few digits of a term are already predicted by earlier terms, as for a(1447) = 19173719153 and a(3868) = 371379371929. So I maintain a list of such "future predictions", and when a digit of a trial candidate is not compatible with one of these, I jump to the smallest larger candidate which has a different digit in that place (i.e., I let p += 10^k - p % 10^k, where 10^k is the weight of that digit, of course followed by nextprime()).

In particular, when the first digit of a term must be 1, then when that of the candidate reaches 2, it will be increased to 3, 4, ..., 9 and '10' (i.e., the smallest number with one more digit) in 8 steps. Yes, could be done in a single step but that's not essential.(*)

(*)UPDATE: in the mean time, I changed my code as to "jump" directly to the correct digit (if it is larger), using p += max(needed - current digit, 1)*10^k - p % 10^k. [I still increment by 1 if needed < current, don't know yet why it bugs when I use (needed-current)%10 instead of max(...): e.g., for a(187) = 10093 it would give something a little larger (around 101xx) instead; the "active" constraints being a(186) = 7151 and a(185) = 7193, saying that the 1st & 5th digit must be 1 and 3.]

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.