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

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]


• Possibly better suited to Mathematics? Commented Mar 12, 2017 at 14:12
• A minor point, but instead of reversing {First, Length} why not just do {Length, First} ? Commented Mar 12, 2017 at 14:32
• If you're able to get $\lambda$ as a Root object, why not use N to get the digits? Commented Mar 13, 2017 at 14:00
• You could use the continued fraction form. The sequence can be found here oeis.org/A014715. The following (Javascript) code works correctly up to 190 digits: c=[1,3,3,2,2...] // about 200 elements b=r=BigInt(10)**BigInt(200); for (i=c.length-1; i>=0; i--) { r = BigInt(c[i])*b + b*b/r; } // r now contains 10**200 * Conway'c constant Commented Aug 3, 2018 at 10:19

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 = InternalBag[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++
];
InternalStuffBag[b, count];
InternalStuffBag[b, value];
];
InternalBagPart[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.

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.