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.

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