# Find a formula for the recurrent sequence $q_{n+1}=q_n(q_n+1)+1$

Find an analytic formula for the recurrent sequence $$q_{n+1}=q_n(q_n+1)+1,\;\;q_0\in\mathbb N.$$

(The question was asked on 03.05.2018 by M. Pratsovytyi, see page 109 of Volume 1 of the Lviv Scottish Book).

• The question arises: what for? Isn't it art for art's sake? – user64494 Jun 16 '18 at 11:08
• Iterations of $f(x)=ax^2+bx+c$ have explicit formulae when $b^2-4ac=2b+8$, when $f$ is conjugate to Chebyshev polynomial $2x^2-1$ (iterations of which are Chebyshev polynomials of degree $2^n$, they all have explicit formulae). Another situation is when $f$ is similar to $cx^2$. Neither holds here, so I doubt that a formula exists. – Fedor Petrov Jun 16 '18 at 11:25
• @Qfwfq Why did you change "analytic formula" to "explicit formula"? They are not the same. An explicit formula would be something like Binet's formula for the $n$th Fibonacci number, but an analytic formula would suggest something like a convergent power series giving the $n$th term. So I think that you may have changed the meaning of the question, and unless the OP indicates that they agree with your edit, I'd suggest that you revert back to the previous version. – Joe Silverman Jun 16 '18 at 19:23
• @user64494 if complex quadratic dynamics produces a lot of interesting mathematics, is not it natural to ask what happens on a formal level? – Fedor Petrov Jun 17 '18 at 8:17
• @JoeSilverman in the Lviv book original question it is simply "formula". – Fedor Petrov Jun 17 '18 at 8:19

I'm not quite sure what you mean by an "analytic formula." As Fedor Petrov indicated, there is unlikely to be a closed formula. However, there is a convergent power series. More precisely, consider iteration of the function $f(x)=x^2+x+1$ on the attracting basin at $\infty$, which includes in particular all positive integers. There is an invertible power series, called the Böttcher coordinate of $f$ at $\infty$, that conjugate $f$ to the map $z^2$, and thus it conjugates $f^n$ to the map $z^{2^n}$. This gives an analytic formula is the sense that one gets a convergent power series, but there is not a simple formula for, say, the coefficients of the Böttcher coordinate. Never-the-less, the Böttcher coordinate is a key tool in dynamics for studying iteration in a neighborhood of a superattracting fixed point.

Explicitly, there is an invertible power series \begin{aligned} \phi(x) &= x + \frac{1}{2} x^2 + \frac{7}{8} x^3 + \frac{3}{4} x^4 + \frac{183}{128} x^5+ O(x^6), \\ \phi^{-1}(x) &= x - \frac{1}{2} x^2 - \frac{3}{8} x^3 + \frac{13}{16} x^4 - \frac{77}{128} x^5+ O(x^6), \end{aligned} so that $f(x) = 1/\phi\bigl(\phi^{-1}(1/x)^2\bigr)$. Hence if $q\in\mathbb N$, then $$f^n(q) = 1/\phi\bigl(\phi^{-1}(1/q)^{2^n}\bigr).$$

• How explicit is $\varphi$? – Fedor Petrov Jun 16 '18 at 19:24
• @FedorPetrov Depends on your definition of explicit. There's no nice formula for its coefficients, as far as I'm aware. There is a recursive formula for its coefficients. It's easy to compute as many coefficients as your want. – Joe Silverman Jun 16 '18 at 20:14
• Well, but it is easy to compute as many terms $q_n$ as you want either. – Fedor Petrov Jun 16 '18 at 20:39
• If we write the functional equation of $\phi$ in terms of $g(z):=1/f(1/z)=z^2/(1+ z+ z^2)$, we get $g(\phi(z))=\phi(z^2)$, that might be used for a recurrence for the coefficients (so in fact $\phi$ conjugates $z^2$ with $g(z)$, or am I wrong?) – Pietro Majer Jun 16 '18 at 20:41
• @FedorPetrov It's true that you can compute as many terms $q_n$ as you want for a given $q_0$. The point of the Bottcher formula is that it's a universal series that you can, in principal, substitute any starting point $q_0$. – Joe Silverman Jun 16 '18 at 21:31

This is a second answer with a somewhat different viewpoint from my other answer. It is an expansion, in some sense, of Gerry Myerson's answer. There is a general theory for estimating these sorts of sequences (even for rational starting points) that goes by the name dynamical canonical heights. I'll restrict attention to this polynomial $f(x)=x^2+x+1$ and integer starting points. For a given starting point $q_0\in\mathbb Z$, the dynamical canonical height of $q_0$ for the map $f$ is given by the limit $$\hat h_f(q_0) := \lim_{n\to\infty} \frac{1}{2^n} \log \bigl| f^{\circ n}(q_0)\bigr|.\qquad(*)$$ Here $f^{\circ n}(q_0)$ is the quantity that the OP called $q_n$. It is a standard fact that the limit converges, and that $$\log|q_n| = 2^n \cdot \hat h_f(q_0) + O(1),$$ where the $O(1)$ is bounded independently of both $q_0$ and $n$. (It is also relatively easy to give an explicit bound for $O(1)$.) The quantity $\hat h_f(q_0)$ is the logarithm of the $c$ values in Gerry's answer.

All of this wouldn't be so useful if we had to use the limit formula $(*)$ to compute $\hat h_f(q_0)$, since $(*)$ already requires us to compute $q_n=f^{\circ n}(q_0)$ for large values of $n$. However, there is a rapidly converging series that computes $\hat h_f(q_0)$. It is a modification of a formula originally due to Tate for elliptic curves and can be found in [1]. Computing $k$ terms of this series gives the value of $\hat h_f(q_0)$ with an error of $O(2^{-k})$, so it is quite feasible to compute $\hat h_f(q_0)$ to 100, or even 1000, decimal places.

[1] Call, Gregory S.; Silverman, Joseph H., Canonical heights on varieties with morphisms, Compos. Math. 89, No. 2, 163-205 (1993). ZBL0826.14015.

The sequence (at any rate, the case $q_0=1$) has been studied, and references are given at OEIS. The closest thing to a formula given there is $a(n) = [c^{2^n}]$ for $n > 0$, where $c = 1.385089248334672909882206535871311526236739234374149506334120193387331772\dots$

• Strictly speaking that is the formula for $q_0=0$; the formula for $q_0=1$ would have $c= 1.385089\ldots^2$, while for $q_0=2$ would be something like $c=2.753832\ldots$ and for $q_0=3$ you have $c= 1.385089\ldots^4$ – Henry Jun 16 '18 at 14:07
• @Henry Do you know how these formulas were obtained? – Shamisen Jun 17 '18 at 22:54
• @Shamisen When I have done similar calculations for sequences where each term is close to the square of the previous term, I simply used arbitrarily high precision calculations of each term and their $2^n$th roots and then checked whether rounding down (as here) or rounding up or something else worked – Henry Jun 17 '18 at 23:48
• @Shamisen See for example A. V. Aho and N. J. A. Sloane, Some doubly exponential sequences, Fib. Quart., 11 (1973), 429-437 and in particular pages 434 and 435 – Henry Jun 18 '18 at 0:03
• @Sham, it looks to me like Joe Silverman's answers and links should give you some idea how to get these formulas. – Gerry Myerson Jun 18 '18 at 0:21

If we denote $A_n=q_n+1/2$, then $$A_n=A_{n-1}^2+5/4$$ with $A_0=q_0+1/2\ge 3/2$ by $q_0\in\mathbb{N}$. Further, $$\log A_n=2\log A_{n-1}+\log\left(1+\frac{5}{4A_{n-1}^2}\right),$$ namely $$\frac{1}{2^n}\log A_n-\frac{1}{2^{n-1}}\log A_{n-1}=\frac{1}{2^n}\log\left(1+\frac{5}{4A_{n-1}^2}\right).$$ Thus $$\log A_n=2^n\left(\log A_{0}+\sum_{k=1}^{n}\frac{1}{2^k}\log\left(1+\frac{5}{4A_{k-1}^2}\right)\right).$$ Clearly, $$A_n> A_{n-1}^{2^1}> A_{n-2}^{2^2}>\cdots A_0^{2^n}\ge (3/2)^{2^n}.$$ Thus, $$0<\sum_{k\ge n+1}\frac{2^n}{2^k}\log\left(1+\frac{5}{4A_{k-1}^2}\right)<\sum_{k\ge n+1}\frac{2^n}{2^k}\frac{5}{4A_{k-1}^2}<\frac{5}{4A_{n}^2}.$$ Hence note that $A_0>3/2$ we obtain that $$1-\frac{5}{4A_n^2}<e^{-\frac{5}{4A_n^2}}\le A_n\kappa^{-2^n}<1$$ with $$\kappa=A_0\prod_{k=1}^{\infty}\left(1+\frac{5}{4A_{k-1}^2}\right)^{\frac{1}{2^k}}$$ a constant depends only on $A_0$. Thus we can prove that $$\kappa^{2^{n}}-3\kappa^{-2^{n}}<A_n<\kappa^{2^n}$$ for all $n\ge 2$ by note that $\kappa>A_0=3/2$. For the computing of $\kappa$, it follows from above that $$A_n^{1/2^n}<\kappa<A_n^{1/2^n}\left(1+\frac{3}{A_n^2}\right)^{1/2^n}<A_n^{1/2^n}\left(1+\frac{3}{2^nA_n^2}\right).$$