# Extracting an asymptotic from a sequence defined by a recurrence relation

Suppose I have a sequence defined via its first term and a recurrence relation involving summation over all previous values with some coefficients. Here is the sequence I am interested in right now (although I would like to find a general method applicable to other sequences as well, if possible): $$a_0=1,\quad a_n=\sum_{k=0}^{n-1}\frac{a_k}{(n-k+1)!\,(2^n-1)}.\tag1$$ This sequence is positive rational, monotone decreasing, decaying faster than exponentially, with a few initial terms: $$1,\,\frac12,\,\frac5{36},\,\frac1{36},\,\frac{143}{32400},\,\frac{19}{32400},\,\frac{1153}{17146080},\,\frac{583}{857 30400},\,\dots\tag2$$

It is related to values of the Fabius function. I'm interested in its asymptotic expansion for $$n\to\infty$$. Using empirical numeric methods, and with some luck, I came to this conjecture: $$\log_2a_n\stackrel{\color{#aaaaaa}?}=-n \log_2n+\frac{n}{\log2}-\frac{\log _2^2n}{2} +\left(1-\frac{1}{2 \log ^22}-\frac{2}{\log2}\right)+\mathcal O\!\left(\frac{\log^2n}{n}\right)\!.\tag3$$ Could you suggest a way to prove it and find a few next terms of this asymptotic expansion?

How should I approach problems like this in general?

Update: The sequence in question is directly connected to the values of the Fabius function at negative integer powers of two: $$F\!\left(2^{-n}\right)=2^{-\binom n2}a_n.\tag4$$

There is also a direct non-recursive formula for it (conjectured): $$a_n\stackrel{\color{#aaaaaa}?}=\frac{(-1)^n}{(2;2)_n}\sum_{k=0}^n\frac{\binom n k_2}{2^{(n-1)k}(n+k)!}\sum_{\ell=0}^{2^k-1}(-1)^{\sigma_2(\ell)}\left(\ell-2^k+\tfrac12\right)^{n+k},\tag5$$ where $$\left(q;q\right)_n=\prod_{k=1}^n\!\left(1-q^k\right)$$ is the q‑Pochhammer symbol, $${\binom n k}_q=\frac{\left(q;q\right)_n}{\left(q;q\right)_k\,\left(q;q\right)_{n-k}}$$ is the q‑binomial coefficient, and $$\sigma_2(\ell)$$ is the sum of binary digits of $$\ell$$; note that $$(-1)^{\sigma_2\left(\ell\right)}$$ is just the signed version of the Thue–Morse sequence.

• Commented Jan 13, 2022 at 4:46

The following might help, but it's not a complete solution.

The coefficients $$a_n$$ have an alternative expression, where the $$B_n$$ are the Bernoulli numbers: $$\tag{1} a_n = \sum_{k=0}^n 2^k a_k \ \frac{B_{n-k}}{(n-k)!}$$ To prove $$(1)$$, start with the recursion given in the problem statement, flip the direction of summation, and pull the $$a_n$$ term into the sum to get $$2^n a_n = \sum_{k=0}^n \frac{1}{k+1} \binom{n}{k} a_{n-k} \frac{(n-k)!}{n!}$$ Set $$c_{n} = n!a_n$$ so that $$2^n a_n = \sum_{k=0}^n \frac{1}{k+1} \binom{n}{k} c_{n-k}$$ This equation is invertible, e.g. Binomial Transform of a Particular Binomial Sum

Reverting back to $$a_n$$ gives the formula $$(1)$$. Define a generating function $$f(x) =\sum_{n \ge 0} a_n x^n$$. Then $$(1)$$ implies the functional equation

$$\tag{2} \quad f(x) = f(2x) \frac{x}{e^x -1}$$

Iterate on $$(2)$$, and formally the solution is

$$\tag{3} \quad f(x) = \prod_{k=0}^\infty \frac{ \exp{\big(x \ 2^{-k-1}\big)} - 1}{x \ 2^{-k-1} }$$

I've calculated many of the $$a_n$$ and checked that $$(3)$$ is satisfied numerically when the product is truncated and $$x$$ is sufficiently small.

$$(3)$$ looks like a formidable equation to do asymptotic analysis on. $$(2)$$ looks like a trick with a Mellin transform might be advantageous. The closest thing I know of that looks like the proposed asymptotic expansion has to do with Rice integrals, and instead of Bernoulli numbers, binomials appear in the sum. See example $$4$$ in $$[1]$$.

Edit 1/13/2022. Summary: the asymptotics of $$a_n$$ can be approximated by a saddle point analysis of the large $$x$$ asymptotics of $$f(x)$$ given above. Here is a short algorithm for calculation:

$$\tag{4} x_0=-\frac{1}{\log{2}}W_{-1}\big( -\log{2} \cdot 2^{-n} \big)$$ $$\tag{5} C=2^{-1/12}\exp{\Big(\big( \gamma_1 + \gamma^2/2-\pi^2/12 \big) }/\log{2} \Big)$$ $$\tag{6} F(x)=C\sqrt{x}\exp{\big(x -\frac{1}{2 \log{2}}\log^2 x \big)}$$ $$b(x)=x \frac{d}{dx}(x \frac{d}{dx} \log{(F(x)} )$$ $$\tag{7} a_n^{A} \sim (2 \pi b(x_0))^{-1/2} F(x_0) x_0^{-n}$$ Notation: $$W_{-1}(x)$$ is the Lambert W function on the '-1 branch' and is implemented in Mathematica as ProductLog[-1,x]. The $$\gamma$$ is Euler's constant and $$\gamma_1$$ is the first Stieltjes constant, $$\gamma_1 \approx$$ -0.0728158.

Here is a table comparing the true $$a_n,$$ the approximation (above) $$a_n^{A},$$ and the proposed approximation, $$a_n^{P} = 2^{ (eq.3)}$$ where eq. 3 corresponds to that in the proposers write up.

$$\begin{array}{c|lcr} n & a_n & a_n^{A} & a_n^{P} \quad \\ \hline 20 & 6.08905 \,\cdot\, 10^{-22} & 6.07191\,\cdot\, 10^{-22} & 9.3967\,\cdot\, 10^{-22} \\ 40 & 1.00137 \,\cdot\, 10^{-52}& 1.00047\,\cdot\, 10^{-52} & 1.3980\,\cdot\, 10^{-52} \\ 60 & 1.30647\,\cdot\, 10^{-87} & 1.30565\,\cdot\, 10^{-87} & 1.7222\,\cdot\, 10^{-87} \\ 80 & 3.13549\,\cdot\, 10^{-125} & 3.13398\,\cdot\, 10^{-125} & 3.9804\,\cdot\, 10^{-125} \end{array}$$ The approximation $$a_n^{A}$$ is excellent, about 3 significant digits for $$n=20,$$ and about 4 digits for $$n=80.$$ A comment before sketching the proof. The proposer's solution seems to have missed the constant, and otherwise the expansion seems to be an asymptotic expansion of the Lambert W function, $$W(x) \sim \log{x} - \log{\log{x}} + \cdots$$

The first part of the proof is to show that, for large $$x$$ $$\tag{4} \quad f(x) = \prod_{k=1}^\infty \frac{ \exp{\big(x \ 2^{-k}\big)} - 1}{x \ 2^{-k} } \sim F(x)\exp{\big(P(x) \big)}$$ where $$F(x)$$ is given by eqs. (5) & (6) and $$P(x)$$ is given by

$$P(x)=\sum_{k=-\infty, k \ne 0}^\infty\ \frac{1}{2\pi i k}\Gamma\big(1-\frac{2 \pi i k}{\log2} \big) \zeta\big(1-\frac{2 \pi i k}{\log2} \big) \exp{\big(x \, \frac{2 \pi i k}{\log2} \big)}.$$

The formulas presented are a slight modification of that presented in Example 9.3 in reference [2]. We are going to ignore $$P(x),$$ in subsequent analysis, mainly because the amplitudes of the oscillations are tiny, with magnitudes of $$7.4\,\cdot\,10^{-7}, 4.4\,\cdot\,10^{-13}\cdots$$ for $$k=1,2...$$ However, if you want even more precise asymptotics (3 or 4 digits is not enough?) then they have to be included.

Reference [2] deals with a problem that has the functional equation $$g(x) = g(2x)\frac{x}{1-e^{-x}}$$ Note how close that is to the functional equation in (3). On taking logarithms we have $$\log{\big(f(x)\big)} = \sum_{k=1}^\infty \log{\Big( \frac{ \exp{\big(x \ 2^{-k}\big)} - 1}{x \ 2^{-k} } \Big) }$$ $$= \sum_{k=1}^\infty \log{\Big(\exp{\big(x \ 2^{-k}\big)}\Big) } + \log{\Big( \frac{ 1-\exp{\big(-x \ 2^{-k}\big)} }{x \ 2^{-k} } \Big) }$$ $$=x+ \sum_{k=1}^\infty \log{\Big( \frac{ 1-\exp{\big(-x \ 2^{-k}\big)} }{x \ 2^{-k} } \Big) }$$ The infinite sum is precisely that given in ref. [2] for $$g(x)$$, so its asymptotics can be copied verbatim, which I have done. Very briefly, the technique is to find the Mellin transform for the sum, then find the asymptotics by moving the contour of integration over the function's singularities. Later on, the reference extracts the asymptotics for the coefficients of $$g(x)$$ by using a process called dePoissionization, but that requires a growth restriction on the function. The addition of $$e^x$$ for my $$f(x)$$ destroys that possibility. However, traditional saddle point works. In that situation, one can use eq. (7) for a generic function $$F(x).$$ How does one determine the saddle point? Look at the Cauchy representation of $$a_n,$$ $$a_n = \frac{1}{2 \pi i} \oint \frac{F(x)}{x^{n+1}} dx$$ $$\rightarrow \frac{1}{2 \pi i} \oint \exp{\big( x - n\log{x} - \frac{1}{2 \log{2}}\log^2{x} \big)} \frac{dx}{\sqrt{x}}$$ The saddle point $$x0$$ is that in which the derivative of the argument of the exponential is zero,

$$\frac{d}{dx} \big(x - n\log{x} - \frac{1}{2 \log{2}}\log^2{x} \big)\overset{!}{=} 0$$ The solution to this equation for $$x=x_0$$ is given in eq. (4).

$$[1]$$ Mellin transform and asymptotics: Finite differences and Rice's integrals, P. Flajolet & R. Sedgewick, Ther. Comp. Sci 144 (1995) 101-124.

$$[2]$$ Average Case Analysis of Algorithms on Sequences W. Szpankowski, Wiley-Interscience, 2001.

• Thanks. Yes, my constant term was pure guesswork. I did not expect Stieltjes constants to appear there, this is a big surprise for me. Commented Jan 13, 2022 at 22:16
• This gives us an asymptotic for the Fabius functions for $x\to 0^+$, right? $F(x)\sim \frac{1}{2^{7/12} \sqrt{\pi x}} \exp \left(\frac{1}{\log 2} \left(\gamma _1+\frac{\gamma ^2}{2}-\frac{\pi ^2}{12}-W_{-1}(-x \log2)-\frac{1}{2} W_{-1}^2(-x \log 2)\right)\right)$. Commented Jan 13, 2022 at 22:19
• @VladimirReshetnikov Hello, I've enjoyed reading some of your posts in the past. I'm not familiar with the Fabius function, but you've added some links that I haven't gone through yet. I'll check them out. It'll be tomorrow at the earliest before I can respond. Commented Jan 13, 2022 at 22:39
• By the way, $\small W_{-1}(x)\sim\log (-x)-\frac{\log (-x) \log (-\log (-x))}{1+\log (-x)}-\sum _{n=2}^{\infty } \sum _{k=1}^n \frac{(-1)^n S_n^{(n-k)} \log ^{k+1}(-\log (-x))}{(k+1)! \log ^n(-x)}$ for $x\to0^-$. Commented Jan 14, 2022 at 6:06
• @VladimirReshetnikov Thanks for pointing out the correct asymptotic form to be used; I meant the asymptotic form to explain the appearance of a $\log{\log}$ term in the empirically discovered formula. I looked at your links on the Fabius function and do not see the immediate link. The functional equation relates a derivative of a function to the same function at double its argument. However, I have no experience with it, so take my comment with a grain of salt. Commented Jan 14, 2022 at 15:34

This is not a proper answer since it involves several heuristics. I will use the generating function $$f(x)$$ found by @skbmoore. For convenience, let $$g(x)=f(2x)$$. I would like to use Hayman's formula for the asymptotics of the Maclaurin coefficients of certain entire functions. The first non-trivial assumption is that $$g(x)$$ is admissible in the sense of Hayman, i.e., we can use his formula. Note that the Maclaurin coefficients of $$g(x)$$ are $$2^n a_n$$.

First, we estimate the growth of the logarithmic derivative. It can be shown that \begin{align*} a(r) :\!&= r\frac{{g'(r)}}{{g(r)}} = 2r + \sum\limits_{k = 0}^\infty {\left( {\frac{{r/2^k }}{{e^{r/2^k } - 1}} - 1} \right)} \\ & = 2r - \frac{{\log r}}{{\log 2}} - \frac{1}{2} - \sum\limits_{k = 1}^\infty {\frac{{2^k r}}{{e^{2^k r} - 1}}} +c(r), \end{align*} where the function $$c(r)$$ is bounded and staisfies $$c(r)=c(2r)$$. Numerics suggests that $$|c(r)|<2\times 10^{-5}$$. Hayman defines the sequence $$r_n$$ via $$a(r_n ) = n$$. Using the above, I found $$r_n = \frac{n}{2} + \frac{{\log n}}{{2\log 2}} - \frac{1}{4} + o(1)$$ as $$n\to +\infty$$. We also have $$b(r) := ra'(r) = 2r + \mathcal{O}(1)$$ for large $$r$$. Now \begin{align*} \log g(r) & = \int_1^r {\frac{{a(t)}}{t}dt} + \log g(1) \\ & = 2r - \frac{{\log ^2 r}}{{2\log 2}} - \frac{{\log r}}{2} + C_0 - \sum\limits_{k = 1}^\infty \log (1 - e^{ - 2^k r} ) + d(r), \end{align*} with some real constant $$C_0$$ and with a bounded function $$d(r)$$ which staisfies $$d(r)=d(2r)$$. Therefore, by Hayman's formula, \begin{align*} \log (2^n a_n ) & = - n\log r_n + \log g(r_n ) - \frac{1}{2}\log (2\pi ) - \frac{1}{2}\log b(r_n ) + o(1) \\ & = - n\log \left( {\frac{n}{2}} \right) + n - \frac{{\log ^2 n}}{{2\log 2}} + C_1 +d(r_n)+ o(1) \end{align*} as $$n\to +\infty$$ with some real constant $$C_1$$. Equivalently, $$\log _2 a_n = - n\log _2 n + \frac{n}{{\log 2}} - \frac{{\log _2^2 n}}{2} + C_2 +d(r_n)+ o(1)$$ as $$n\to +\infty$$ with some real constant $$C_2$$.

• Ah, Hayman;s theorem; very good. It should get the first few terms right. The only thing that bothers me is that terms even smaller, like the OP wants, might exhibit a small amplitude oscillation, which is what happens when the normalize Bernoulli number is replaced by a binomial. Commented Jan 12, 2022 at 20:16
• @Gary Thanks, Dr. Nemes. This is very valuable input. I always appreciate reading your analysis of asymptotics. Commented Jan 13, 2022 at 2:14
• @VladimirReshetnikov Thank you! If the observations in my comment above are correct, my asymptotics becomes the weaker $$\log _2 a_n = - n\log _2 n + \frac{n}{{\log 2}} - \frac{{\log _2^2 n}}{2} + \mathcal{O}(\log _2 n),$$ where the implied constant is small (somthing of the size $10^{-5}$). So if there is indeed a term of size $\log_2 n$, it is difficult to detect it unless $n$ is extremely large.
– Gary
Commented Jan 13, 2022 at 2:25
• Initially, I also expected $\mathcal{O}(\log _2 n)$ term to be present as the next natural order of growth. But all numerical experiments strongly suggest that it has $0$ coefficient, and the next is a constant term, followed by $\mathcal O\!\left(\frac{\log^2n}{n}\right)$. I tentatively conjecture the latter is $-\frac1{2\log^3 2}\frac{\log^2n}{n}$. Commented Jan 13, 2022 at 2:46
• One would need to compute $a_n$ for $n$s of size $n \approx 10^{10^5 }$ to see the $\log_2 n$ term. It is not exactly $\log_2 n$ though, it is some more complicated function of that size coming from the term $d(r)$ above.
– Gary
Commented Jan 13, 2022 at 3:46