Calculate $140$ digits of Conway's Constant from the Look and Say Sequence

Question

The look-and-say sequence is the sequence of numbers $1, 11, 21, 1211, 111221, 312211, …,$ in which each term is constructed by “reading” the previous term in the sequence. For example, the term $1$ is read as “one $1$”, which becomes the next term: $11$. Then $11$ is read as “two ones”, which becomes the next term: $21$, and so on...

then

If $L_n$ is the number of digits in the $n$th term in the sequence, then

$$\lim_{n \to \infty} \frac{L_{n+1}}{L_n}=\lambda$$

where $λ = \color{red}{1.303577269034...}$ is a Conway's Constant. ( algebraic number of degree $71$)

I tried to calculate it , but its very slow :

RunLengthEncode[x_List] := 
RunLengthEncode[x] = (Through[{First, Length}[#1]] &) /@ Split[x]
LookAndSay[n_, d_: 1] := 
NestList[Flatten[Reverse /@ RunLengthEncode[#]] &, {d}, n - 1]

nthdigit[n_] := Length /@ LookAndSay[n] // Last

(nthdigit[60]/nthdigit[59])~N~15 // AbsoluteTiming

$\{46.6713, \color{red}{1.3035}5090959742\}$

How can I speed up this calculation to find $140$ digits of Conway's Constant ?

$140$ digits of $\lambda$ is quite enough to find polynomial

lambda = 1.3035772690342963912570991121525518907307025046594048757548613906285508
878524615571268157668644252255534713930470949026839628498935515543473758248566910891;


 MinimalPolynomial[RootApproximant[lambda], x]

Answer 1

I call this ambiguous and I think you haven't thought this through. Let us assume it would be possible by such a direct method of calculating $\lambda$. First, we need to speed up the high-level function of yours that you probably have from Rosetta code. If the numbers grow larger and larger, Split becomes a bottleneck, and I don't now how to speed it up except writing in in a low-level compiled function.

Therefore, here is an implementation calculating the run-length encoding in an ugly combination of two loops:

rl = Compile[{{l, _Integer, 1}},
  Module[{b = Internal`Bag[Most[{0}]], n = Length[l], count, value, 
    i = 1},
   While[i <= n,
    count = 1;
    value = l[[i++]];
    While[i <= n && value == l[[i]],
     count++;
     i++
     ];
    Internal`StuffBag[b, count];
    Internal`StuffBag[b, value];
    ];
   Internal`BagPart[b, All]
   ], CompilationTarget -> "C", RuntimeOptions -> "Speed"
  ]

As comparison let's take a slightly adapted version of the core of your function

RunLengthEncode[x_List] := Flatten[Through[{Length, First}[#1]] & /@ Split[x]]

While yours takes

res1 = Nest[RunLengthEncode, {1}, 60]; // AbsoluteTiming
(* {35.2265, Null} *)

the compiled version only needs

res2 = Nest[rl, {1}, 60]; // AbsoluteTiming
(* {0.92388, Null} *)

and the results are the same

res1 == res2
(* True *)

Speed, unfortunately, is not our biggest problem. Let me demonstrate the byte size of the results for all sequences up to the 50th.

sizes = Table[ByteCount[Nest[rl, {1}, n]], {n, 1, 50}];

and fit an exponential model on the data

nlm = NonlinearModelFit[sizes, a*Exp[b x], {a, b}, x]

Without quantifying the quality of the model fit, let us instead just look on the combined plots:

Show[
 ListLinePlot[sizes, PlotStyle -> Thick],
 Plot[nlm[x], {x, 0, 50}, PlotStyle -> Directive[Dotted, Red]]
 ]

That looks promising if you like to eat up all your system memory, but since we are a very optimistic bunch, let's estimate at which sequence we need 16GB of precious RAM

FindRoot[Normal[nlm]/2^30 - 16, {x, 60}]
(* {x -> 78.3621} *)

This means, that for 79 steps, we will cross the 16GB line. I like to remind you, that you are using memoization to store all previous results, so your computer is going to burst into flames sooner. I have 32GB, so mine is probably dying a bit later. Let's be cautious and calculate n=75

seq75 = Nest[rl, {1}, 75]; // AbsoluteTiming
(* {44.1891, Null} *)

This monster already uses 6.6GB of RAM and has 881.752.750 elements. Let's calculate $\lambda$ using the 76th sequence

N[Length[rl[seq75]]/Length[seq75], 15]
(* 1.30358560492156 *)

Still as bad as your approximation. Therefore, I don't think that speeding your approach up will help you. I believe a completely different method is required to get $\lambda$ with more precision.

Answer 2

Carrying out Greg Hurst's suggestion, Conway's constant is the largest real root of

$$x^{71}-x^{69}-2 x^{68}-x^{67}+2 x^{66}+2 x^{65}+x^{64}-x^{63}-x^{62}-x^{61}-x^{60}-x^{59}+2 x^{58}+5 x^{57}+3 x^{56}-2 x^{55}-10 x^{54}-3 x^{53}-2 x^{52}+6 x^{51}+6 x^{50}+x^{49}+9 x^{48}-3 x^{47}-7 x^{46}-8 x^{45}-8 x^{44}+10 x^{43}+6 x^{42}+8 x^{41}-5 x^{40}-12 x^{39}+7 x^{38}-7 x^{37}+7 x^{36}+x^{35}-3 x^{34}+10 x^{33}+x^{32}-6 x^{31}-2 x^{30}-10 x^{29}-3 x^{28}+2 x^{27}+9 x^{26}-3 x^{25}+14 x^{24}-8 x^{23}-7 x^{21}+9 x^{20}+3 x^{19}-4 x^{18}-10 x^{17}-7 x^{16}+12 x^{15}+7 x^{14}+2 x^{13}-12 x^{12}-4 x^{11}-2 x^{10}+5 x^9+x^7-7 x^6+7 x^5-4 x^4+12 x^3-6 x^2+3 x-6.$$

I found that formula from the OEIS page. So the code

cc=Max[x /. Solve[-6 + 3 x - 6 x^2 + 12 x^3 - 4 x^4 + 7 x^5 - 7 x^6 + 
 x^7 + 5 x^9 - 2 x^10 - 4 x^11 - 12 x^12 + 2 x^13 + 7 x^14 + 
 12 x^15 - 7 x^16 - 10 x^17 - 4 x^18 + 3 x^19 + 9 x^20 - 7 x^21 - 
 8 x^23 + 14 x^24 - 3 x^25 + 9 x^26 + 2 x^27 - 3 x^28 - 10 x^29 - 
 2 x^30 - 6 x^31 + x^32 + 10 x^33 - 3 x^34 + x^35 + 7 x^36 - 
 7 x^37 + 7 x^38 - 12 x^39 - 5 x^40 + 8 x^41 + 6 x^42 + 10 x^43 - 
 8 x^44 - 8 x^45 - 7 x^46 - 3 x^47 + 9 x^48 + x^49 + 6 x^50 + 
 6 x^51 - 2 x^52 - 3 x^53 - 10 x^54 - 2 x^55 + 3 x^56 + 5 x^57 + 
 2 x^58 - x^59 - x^60 - x^61 - x^62 - x^63 + x^64 + 2 x^65 + 
 2 x^66 - x^67 - 2 x^68 - x^69 + x^71 == 0, x, Reals]]

produces a Root object. If you call N[cc,140], you get

1.3035772690342963912570991121525518907307025046594048757548613906285508878524615571268157668644252255534713930470949026839628498935515543474

I think 'Calculate' implies computations easily derived from the statement of the problem, as opposed to doing a lot of math to simplify the problem. So I don't expect an Accept, I just felt like doing this computation.