Generate a sequence of primes giving position and value of a digit in its digit-decomposed version

Question

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.

Answer 1

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 *)

Answer 2

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.]