%stern.tex: a Plain TeX file by Shalosh B. Ekhad and Doron Zeilberger
%Automated Generation of Generating Functions Related to Generalized Stern's Diatomic Arrays in the footsteps of Richard Stanley
%begin macros
\def\L{{\cal L}}
\baselineskip=14pt
\parskip=10pt
\def\halmos{\hbox{\vrule height0.15cm width0.01cm\vbox{\hrule height
0.01cm width0.2cm \vskip0.15cm \hrule height 0.01cm width0.2cm}\vrule
height0.15cm width 0.01cm}}
\font\eightrm=cmr8 \font\sixrm=cmr6
\font\eighttt=cmtt8
\magnification=\magstephalf
\def\A{{\cal A}}
\def\B{{\cal B}}
\def\C{{\cal C}}
\def\S{{\cal S}}
\def\P{{\cal P}}
\def\Q{{\cal Q}}
\def\1{{\overline{1}}}
\def\2{{\overline{2}}}
\parindent=0pt
\overfullrule=0in
\def\Tilde{\char126\relax}
\def\frac#1#2{{#1 \over #2}}
%\headline={\rm \ifodd\pageno \RightHead \else \LeftHead \fi}
%\def\RightHead{\centerline{
%Title
%}}
%\def\LeftHead{ \centerline{Doron Zeilberger}}
%end macros
\centerline
{\bf
Automated Generation of Generating Functions Related to Generalized Stern's Diatomic Arrays
}
\centerline
{\bf
in the footsteps of Richard Stanley
}
\bigskip
\centerline
{\it Shalosh B. EKHAD and Doron ZEILBERGER}
\bigskip
{\bf Abstract}: Using {\it Symbolic Dynamic Programming} we describe algorithms, fully implemented in Maple,
for automatically generating generating functions introduced by Richard Stanley in his study
of generalized Stern arrays, generalized even further, to arrays defined in terms of general
sequences satisfying linear recurrences with constant coefficients, rather than just the Fibonacci
and k-bonacci sequences.
{\bf Appetizer: A Computational Challenge}
Mathematics is so sensitive to `{\it initial conditions}'. In other words the function
{\bf Mathematical Question} $\rightarrow$ {\bf Difficulty of Finding The Answer to that question}
is very `chaotic'. Just `tweaking' an easy question by a {\it proverbial} $\epsilon$ may change it from
{\it tractable} (and sometimes trivial) to (very likely) {\it intractable} (but proving the
intractability for {\it sure} may also be intractable!).
Here are a few of our favorite examples.
$\bullet$ We all know how to prove that there are infinitely many primes, but just
insert the word {\it twin} in front of `primes' and no one (yet) can prove it.
$\bullet$ We all know how to compute the number of $10000$-step simple walks in the two dimensional square lattice ($4^{10000}$),
but stick {\it self-avoiding} in front of {\it walks} and we are all stumped.
$\bullet$ We all know, since Levi Ben Gerson, (exactly!) $700$ years ago, to compute the
number of {\it permutations} of size $10000$, namely $10000!$, a certain 35660-digit integer
that Maple can compute (and display!) in $0.002$ seconds. Now stick the phrase
{\it $1234$-avoiding} in front of {\it permutations} and it would not take much longer,
(using the well-known second-order linear recurrence).
Now {\it transpose} the $2$ and the $3$ in ``$1234$-avoiding'',
and we are all stumped on computing {\it exactly} the integer that counts the number
of $1324$-avoiding permutations of length $10000$.
$\bullet$ The proof of the five-color theorem is half-page long, but the proof of the four-color theorem is
several-thousands-page long.
And we can come up with lots of other illustrations of the {\it sensitivity} of the ``difficulty'' function.
\vfill\eject
Here is yet another example.
{\bf Easy Problem}
Let
$$
F_n(x) \,:=\, \prod_{i=0}^{n-1} \left ( 1 \,+\, x^{2^i} \, + \, x^{2^{i+1}} \right ) \quad,
$$
and write
$$
F_n(x) \, = \, \sum_{k \geq 0} \, a(n,k) \, x^k \quad .
$$
Now let
$$
v(n) := \, \sum_{k \geq 0} \, a(n,k)^2 \quad .
$$
Find $v(10000)$.
Note that it is {\it hopeless} to compute this number from the definition. The degree, and hence the size, of the
polynomials $F_n(x)$ grow exponentially (in fact, it is $2(2^{n}-1)$).
Naively using the definition, we would have to add the squares of $2^{10001}-1$ numbers.
Using [S1], to be reviewed and extended later in this article, Maple can find $v(10000)$, a certain $6591$-digit number,
in a {\it split second.}
Now comes our {\it tweak}
{\bf Hard Problem (at least for us)}
Let
$$
G_n(x) \,:=\, \prod_{i=0}^{n-1} \left ( 1 \,+\, x^{2^i+1} \, + \, x^{2^{i+1}+1} \right ) \quad,
$$
and write
$$
G_n(x) \, = \, \sum_{k \geq 0} \, b(n,k) \, x^k \quad .
$$
Let
$$
w(n):=\sum_{k \geq 0} b(n,k)^2 \quad .
$$
Find the {\bf exact value} of $w(10000)$.
One of us (DZ) is pledging a donation of $100$ US dollars to the OEIS (On-Line Encyclopedia of Integer Sequences)
in honor of the first (correct) solver of this very concrete problem.
Using the definition, Maple can compute the first $26$ values, starting at $n=0$:
$$
1, 3, 13, 55, 233, 1033, 4359, 19081, 83653, 363973, 1604755, 7071677, 31361931, 139661731, 623089471,
$$
$$
2788501361, 12507807967, 56197511503, 252874682743,1139273972183, 5137458451565,
$$
$$
23186535210405, 104711215601401, 473121563716987, 2138654595620755, 9670566829508677 \quad .
$$
Unlike the previous problem, of computing $v(10000)$, for which it is easy to detect a simple ``pattern'' (see below) (and also easy to prove it
rigorously, also see below), and then easily deduce the $10000$-th term,
the modified problem, of computing $w(10000)$, seems much harder.
At any rate, the algorithms described later in the present article, that can handle everything in [S1] and [S2],
and much more, fail miserably on this innocent problem. Of course, it is possible that there
{\bf exists} another algorithm that would make computing $w(10000)$ possible, but we have no clue,
and would love to know!
{\bf The original Stern Diatomic Array}
In the delightful article [S1], Richard Stanley's starting point was the double sequence $a(n,k)$ defined by
$$
F_n(x) \, = \, \sum_{k \geq 0} \, a(n,k) \, x^k \, = \, \prod_{i=0}^{n-1} \left ( 1 \,+\, x^{2^i} \, + \, x^{2^{i+1}} \right ) \quad,
$$
and he was interested in the sequences (for positive integers $r$)
$$
u_r(n) \, := \, \sum_{k \geq 0} a(n,k)^r \quad .
$$
More generally for $\alpha=(\alpha_0, \dots, \alpha_{m-1})$ (where the $\alpha_i$ are non-negative integers), the
sequences
$$
u_\alpha(n) := \sum_{k \geq 0} a(n,k)^{\alpha_0}\, a(n,k+1)^{\alpha_1} \dots a(n,k+m-1)^{\alpha_{m-1}} \quad .
$$
He proved that all these sequences are $C$-finite, i.e., satisfy a {\it linear recurrence equation with constant coefficients},
or equivalently ([Z],[KP]) that their generating function
$$
F_\alpha(x) \, := \, \sum_{n=0}^{\infty} u_\alpha(n) \, x^n \quad,
$$
is a {\bf rational function} of $x$.
In order to illustrate the general theory, he {\it humanly} proved that
$$
F_2(x)= \frac{1-2x}{1-5x+2x^2} \quad .
$$
We will now redo, in {\it excruciating detail}, his proof, {\it our way}, in order to motivate the algorithm that will come later.
Our notation is a little different than the one in [S1], but the bottom line is the same.
Since
$$
F_n(x) \, = \left ( 1 \,+\, x^{2^{n-1}} \, + \, x^{2 \cdot 2^{n-1}} \right ) F_{n-1}(x) \quad,
$$
we have the recurrence that relates the entries in the $n$-row of Stern's array to those of the previous one.
$$
a(n,k) \, = \, a(n-1 \, , \, k) \,+ \,a(n-1 \, , \, k- 2^{n-1}) \,+ \, a(n-1 \, , \, k - 2 \cdot 2^{n-1}) \quad .
$$
For reasons to be made clear later on, let's call $u_2(n)$, $f_{0}(n)$
Incorporating this in the definition of $u_2(n)$ (alias $f_{0}(n)$) we have
$$
f_{0} (n) \, = \, \sum_{k \geq 0}\, a(n,k)^2 \, = \,
$$
$$
\sum_{k \geq 0}
$$
$$
\, \left (a(n-1 \, , \, k) \, + \, a(n-1 \, ,\, k- 2^{n-1}) \, + \,a(n-1 \, , \, k- 2 \cdot 2^{n-1}) \right ) \cdot
$$
$$
\left (a(n-1 \, , \, k) \, + \, a(n-1 \, ,\, k- 2^{n-1}) \, + \,a(n-1 \, , \, k- 2 \cdot 2^{n-1}) \right )
$$
$$
\, = \, \sum_{k \geq 0} \, a(n-1,k) \cdot a(n-1,k)
\eqno(1.1)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k) \cdot a(n-1,k-2^{n-1})
\eqno(1.2)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k) \cdot a(n-1,k- 2 \cdot 2^{n-1})
\eqno(1.3)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k-2^{n-1}) \cdot a(n-1,k)
\eqno(1.4)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k-2^{n-1}) \cdot a(n-1,k-2^{n-1})
\eqno(1.5)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k-2^{n-1}) \cdot a(n-1,k- 2 \cdot 2^{n-1})
\eqno(1.6)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k- 2\cdot 2^{n-1}) \cdot a(n-1,k)
\eqno(1.7)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k-2 \cdot 2^{n-1}) \cdot a(n-1,k-2^{n-1})
\eqno(1.8)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k- 2 \cdot 2^{n-1}) \cdot a(n-1,k- 2 \cdot 2^{n-1}) \quad .
\eqno(1.9)
$$
Note that since the degree of $F_{n-1}(x)$ is $2(2^{n-1}-1)$, $(1.3)$ and $(1.7)$ are $0$.
Also note that by a `shift of the discrete variable $k$', $(1.1)$, $(1.5)$ and $(1.9)$ are the same.
By the commutativity of multiplication, $(1.2)$ and $(1.4)$ are identical, as are $(1.6)$ and $(1.8)$,
and by shifting the variable of summation, we see that these two pairs are identical to each other. Hence
$$
f_0(n) \, = \,3\, \sum_{k \geq 0}\, a(n-1 \, , \, k)^2 \, + \, \, 4\, \sum_{k \geq 0}\, a(n-1 \, , \, k) \, a(n-1 \, , \, k-2^{n-1}) \quad .
$$
The first sum is $f_0(n-1)$ , but the second sum is a new creature, let's call if $f_1(n-1)$, where
$$
f_1(n) \, := \, \sum_{k \geq 0}\, a(n \, , \, k) \, a(n \,, \, k-2^{n}) \quad.
$$
So far we have
$$
f_0(n) \,= \, 3 f_0(n-1) + 4 f_1(n-1) \quad .
$$
We are {\bf forced} to consider $f_1(n)$.
We have
$$
f_{1} (n) \, = \, \sum_{k \geq 0}\, a(n \,, \, k)\, a(n \, , \, k - 2^n)\, =
$$
$$
\, = \, \sum_{k \geq 0} \,
\left( a(n-1 \,, \, k) + a(n-1, k-2^{n-1})+ a(n-1,k- 2 \cdot 2^{n-1} \right ) \cdot
$$
$$
\left( a(n-1 \,, \, k- 2 \cdot 2^{n-1}) + a(n-1, k- 3 \cdot 2^{n-1})+ a(n-1,k- 4 \cdot 2^{n-1}) \right )
$$
$$
\, = \, \sum_{k \geq 0} \, a(n-1,k) \cdot a(n-1,k- 2 \cdot 2^{n-1})
\eqno(2.1)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k) \cdot a(n-1,k- 3 \cdot 2^{n-1})
\eqno(2.2)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k) \cdot a(n-1,k- 4 \cdot 2^{n-1})
\eqno(2.3)
$$
$$
+ \, \sum_{k \geq 0} \, a(n-1,k-2^{n-1}) \cdot a(n-1,k- 2 \cdot 2^{n-1})
\eqno(2.4)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k-2^{n-1}) \cdot a(n-1,k- 3 \cdot 2^{n-1})
\eqno(2.5)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k-2^{n-1}) \cdot a(n-1,k- 4 \cdot 2^{n-1})
\eqno(2.6)
$$
$$
+ \, \sum_{k \geq 0} \, a(n-1,k -2 \cdot 2^{n-1} ) \cdot a(n-1,k- 2 \cdot 2^{n-1})
\eqno(2.7)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k - 2 \cdot 2^{n-1}) \cdot a(n-1,k- 3 \cdot 2^{n-1})
\eqno(2.8)
$$
$$
+\, \sum_{k \geq 0} \, a(n-1,k - 2 \cdot 2^{n-1}) \cdot a(n-1,k- 4 \cdot 2^{n-1}) \quad .
\eqno(2.9)
$$
Once again since the degree of $F_{n-1}(x)$ is $2(2^{n-1}-1)$, $(2.1)$, $(2.2)$, $(2.3)$, $(2.5)$, $(2.6)$ and $(2.9)$ vanish.
By shifting the summation variable $k$, both $(2.4)$ and $(2.8)$ equal $f_1(n-1)$ while $(2.7)$ is our old friend $f_0(n-1)$.
Hence we get, in addition to the previous equation, the following one:
$$
f_1(n) \,= \, f_0(n-1) + 2 f_1(n-1) \quad .
$$
Yea! We did not encounter any new `uninvited guests'. Defining the
generating functions
$$
F_0(x)= \sum_{n=0}^{\infty} \, f_0(n) x^n \quad ,
$$
$$
F_1(x)= \sum_{n=0}^{\infty} \, f_1(n) x^n \quad,
$$
the above two recurrences translate to a {\bf system of two linear equations} in the two {\bf unknowns}
(also using the {\it initial conditions} $f_0(0)=1$, $f_1(0)=0$)
$$
F_0(x)=1+ x(3 F_0(x)+ 4 F_1(x)) \quad, \quad F_1(x)=0+ x( F_0(x)+ 2 F_1(x)) \quad.
$$
Of course, any seventh-grader can solve this system, but why not use Maple?
Typing
{\tt latex(solve($\{$ F0=1+x*(3*F0+4*F1), F1=x*(F0+2*F1) $\}$,$\{$ F0,F1 $\}$ ));}
we get
$$
\left\{ {\it F0}=-{\frac {2\,x-1}{2\,{x}^{2}-5\,x+1}},{\it F1}={\frac {x}{2\,{x}^{2}-5\,x+1}} \right\} \quad.
$$
Confirming that indeed $\sum_{n \geq 0} u_2(n) x^n$ equals $\frac{1-2x}{1-5x+2x^2}$ as claimed in [S1].
(Two lines below equation $(3)$ there).
Let's try and understand what is going on here. We started with an {\bf object of desire}, $f_0(n)$, and we were hoping to
relate it to $f_0(n-1)$. Alas, we were {\bf forced} to consider an {\it uninvited guest}, $f_1(n)$.
Analyzing $f_1(n)$, we were able to express it in terms of $f_1(n-1)$ and $f_0(n-1)$, and luckily, there
were no new `uninvited guests'. Taking the {\it $z$-transforms}, we got a system of two equations and two unknowns
and solving them, gave us the generating function of $f_0(n)$, that we called $F_0(x)$, as well as
the generating function of $f_1(n)$, that we called $F_1(x)$. Of course we can
ungratefully disregard $F_1(x)$ at the end, if we wish, but we needed it in order to computer $F_0(x)$.
Let's now consider the general problem. Suppose that you have {\bf arbitrary} positive integers
$$
0 \leq c_1 \leq c_2 \leq \dots \leq c_r \quad .
$$
We need the sequence
$$
\sum_{k \geq 0} a(n,k+c_1) \, a(n,k+c_2) \, \cdots \, a(n,k+c_r) \quad .
$$
(Note that one can always take $c_1=0$.)
Using the recurrence
$$
a(n,k)=a(n-1,k)+a(n-1,k-2^{n-1})+a(n-1,k-2 \cdot 2^{n-1}) \quad,
$$
above, this sum can be written as a sum of $3^r$ sums, each of them having the form
$$
f[d_1, \dots, d_r; \beta_1, \dots, \beta_r] (n-1) \, =
$$
$$
\sum_{k \geq 0} a(n-1,k+d_1- \beta_1 2^{n-1}) \,\, a(n-1,k+d_2 - \beta_2 2^{n-1}) \,\, \cdots \, a(n-1,k+d_r- \beta_r 2^{n-1} ) \quad .
$$
Here $(d_1, \dots, d_r)$ is a permutation of $(c_1, \dots, c_r)$ and $\beta_1, \dots, \beta_r$ are non-negative integers.
This {\bf forces} us to consider $f[d_1, \dots, d_r; \beta_1, \dots, \beta_r] (n)$ {\it in general}, not just the
initial case of $\beta_1=0, \dots, \beta_r=0$.
So we need to be able to express the general quantity
$$
f[d_1, \dots, d_r \,;\, \beta_1, \dots, \beta_r] (n):=
$$
$$
\sum_{k \geq 0} a(n,k+d_1- \beta_1 2^{n}) \,\, a(n,k+d_2 - \beta_2 2^{n}) \, \, \cdots \,\, a(n,k+d_r- \beta_r 2^{n} ) \quad ,
$$
corresponding to the ``state'' $[d_1, \dots, d_r \,; \, \beta_1, \dots, \beta_r]$ in terms of other such ``states''.
We have to teach the computer:
$\bullet$ How to decide whether such a state is {\it dead on arrival}, i.e. identically zero (like $(1.3)$, $(1.7)$, $(2.1$,
$(2.2)$, $(2.3)$, $(2.5)$,$(2.6)$, and $(2.9)$ in the example above).
$\bullet$ How to {\it automatically} express each of these in terms of other such creatures.
Luckily, computer algebra comes to the rescue! Using the commutativity of multiplication, we can get a {\it canonical form} of each
state, by sorting the list of pairs
$$
[d_1, \beta_1; d_2, \beta_2; \dots ; d_r, \beta_r] \quad,
$$
such that $\beta_1 \leq \beta_2 \leq \dots \leq \beta_r$ and $\beta_1=0$.
Of course as the $\beta_i$ change places, they must bring with them their corresponding $d_i$.
To each such state we assign the {\bf monomial}
$$
x_1^{d_1} \cdots x_r^{d_r} X_1^{\beta_1} \cdots X_r^{\beta_r} .
$$
First one has to replace $X_i$ by $X_i^2$, since $k+d-\beta 2^{n-1}$ in the $n-1$ level becomes
$k+d-\beta 2^{n}\,=\, k+d- (2 \beta) 2^{n-1}$ at the $n$ level.
We let Maple multiply this (converted) monomial by
$$
\prod_{i=1}^{r} (1+X_i+X_i^2) \quad ,
$$
and {\it expand} it, thereby expressing it as a sum of similar-looking monomials.
We replace each monomial by its {\it canonical form} (sorting the powers of $X_i$ and permuting the corresponding
$x_i$ to move-along with their corresponding $X_i$, see Maple source-code).
Finally if the power of $X_1$ in the converted monomial is larger than $0$ we subtract
it from all the powers of $X_i$ (in other words if the power of $X_1$ is, $e$, we divide the monomial by
$X_1^{e} \cdots X_r^{e}$). This corresponds to shifting the variable $k$ in the sum corresponding to the given state.
Now each of these converted monomials corresponds to a state.
It is easy to see that there are only finitely many states (by the upper bound for the degree of $F_n(x)$),
so this process is {\it guaranteed} to {\it terminate}, and then
Maple automatically sets up a system of linear equations, that it can solve {\it all by itself}.
This approach works not just for the original Stern array, but for the more general scenario (introduced in [S1])
$$
F_n(x) \, = \, P(x) \prod_{i=0}^{n-1} Q( x^{b^i} ) \quad,
$$
for {\it any} polynomial $P(x)$ (before $P(x)$ was $1$), and {\it any} polynomial $Q(X)$, (before $Q(x)$ was $1+X+X^2$)
and for {\it any} integer $b \geq 2$ (before $b$ was $2$).
The computer finds the equation for each
`still-to-do' state, by forming its corresponding monomial,
then replacing each $X_i$ by $X_i^{b}$ (due to the transition from the $n$ level to the $n-1$ level),
and then multiplying the resulting monomial by
$$
\prod_{i=1}^{r} Q(X_i) \quad,
$$
expanding, taking the canonical forms, and discarding the monomials that
are `dead-on-arrival'.
This is fully implemented in procedure {\tt RS} in the Maple package {\tt StanleyStern.txt} available from the front of this article
{\tt https://sites.math.rutgers.edu/\~{}zeilberg/mamarim/mamarimhtml/stern.html} \quad ,
or directly from
{\tt https://sites.math.rutgers.edu/\~{}zeilberg/tokhniot/StanleyStern.txt} \quad .
The function-call is
{\tt RS(P,Q,b,A,x,t)} \quad,
where {\tt P,Q,b} are as above and {\tt A} is $[\alpha_0, \dots, \alpha_{m-1}]$ featuring in Stanley's definition of $u_\alpha(n)$
mentioned above.
For example, to get the generating function of what Stanley [S1] calls $u_5(n)$, in the variable $t$
(rather than the variable $x$ that he is using), type
{\tt latex(RS(1+x+x**2,1,2,[5],x,t));} \quad \quad ,
immediately getting
$$
{\frac {20\,{t}^{2}+11\,t-1}{47\,{t}^{2}+14\,t-1}} \quad .
$$
Procedure {\tt RS} {\it dynamically} collects all the needed quantities and by repeatedly invoking procedure {\tt Eq}
{\it dynamically} generates a system of linear equations, until it does not encounter any {\it new guys}.
Of course, we know, {\it a priori} (as proved in [S1] and also follows from
our algorithm) that this must terminate, but even if we did not know that fact, we could have
set an {\bf upper limit} to the number of equations that we are willing to solve, and return {\tt FAIL} if
it got exceeded.
Since we have a {\it theoretical guarantee} that a rational function exists, and we can {\bf bound} the
degree of the denominator, we can use an {\it empirical} approach, and collect enough data
(as Stanley [S1] did in simple cases) and then fit it into the generating function (using something similar
to Maple's {\tt gfun[listtorec]}, but we prefer our own home-made version). This is implemented in procedure
{\tt RSe(p,q,b,A,N,x,t)} \quad,
where $N$ is the number of data points used. Note that while for small cases, where the
order of the recurrences is expected to be relatively low,
this is much faster, but as noted in the above {\it appetizer}, this can't go very far,
since the polynomials $F_n(x)$ have exponential-size degree.
Nevertheless to compute generating functions for what Stanley calls $u_r(n)$ in [S1],
for small $r$ it works very fast. For example, to get the generating function of $u_{10}(n)$, type:
{\tt latex(RSe(1+x+x**2,1,2,[10],15,x,t));}
getting, in $0.16$ seconds
$$
-{\frac {4\,{t}^{4}+1852\,{t}^{3}+7945\,{t}^{2}+96\,t-1}{ \left( t+1 \right) \left( 4\,{t}^{4}-200\,{t}^{3}-9601\,{t}^{2}-100\,t+1 \right) }} \quad .
$$
Using the non-guessing approach (using the algorithm outlined above), i.e. typing
{\tt latex(RS(1+x+x**2,1,2,[10],x,t));} \quad ,
gives that same thing, but in $18.26$ seconds.
On the other hand, for the generating function of what is called $u_{11111}(n)$ in [S1], in other words, the generating function
of the sequence
$$
\sum_{k \geq 0} a(n,k)a(n,k+1)a(n,k+2)a(n,k+3)a(n,k+4) \quad,
$$
typing (using the {\it non-guessing} approach)
{\tt latex(RS(1+x+x**2,1,2,[1,1,1,1,1],x,t));} \quad ,
yields, in $4.3$ seconds
$$
-4\,{\frac { \left( 4\,{t}^{4}-55\,{t}^{3}-69\,{t}^{2}-21\,t-3 \right) {t}^{2}}{ \left( t-1 \right) ^{3} \left( 47\,{t}^{2}+14\,t-1 \right) }} \quad ,
$$
while using the {\it guessing approach}, typing
{\tt latex(RSe(1+x+x**2,1,2,[1,1,1,1,1],20,x,t));} \quad ,
yields the same thing in twice as long.
Eventually, the `guessing approach' will explode of course (as already mentioned before).
{\bf Sample Output files for the Maple package {\tt StanleyStern.txt}}
For the generating functions of [S1]'s $u_r(n)$ for $r \leq 25$, using the {\it guessing} (yet rigorous) approach see the
output file
{\tt https://sites.math.rutgers.edu/\~{}zeilberg/tokhniot/oStanleyStern1eA.txt} \quad .
For the generating functions of [S1]'s $u_{1^r}(n)$, in other words of the sequences
$$
\sum_{k \geq 0} \, a(n,k) \, a(n,k+1) \, \cdots a(n,k+r-1) \quad,
$$
for $r \leq 9$, using the {\it guessing} (yet rigorous) and non-guessing approaches respectively, see the
output files
{\tt https://sites.math.rutgers.edu/\~{}zeilberg/tokhniot/oStanleyStern2e.txt} \quad ,
and
{\tt https://sites.math.rutgers.edu/\~{}zeilberg/tokhniot/oStanleyStern2.txt} \quad .
Note that the latter took quite a big longer.
{\bf Extension to Generalized Stern Arrays Induced by C-finite sequences}
In [S2], Richard Stanley extended his study to analogs where instead of $x^{b^i}$ that feature in the
definition of $F_n(x)$ above one has powers of the form $x^{F_i}$, and more generally, $x^{F^{(k)}_i}$
where $F_i$ are the Fibonacci numbers, and $F^{(k)}_i$ are the $k$-bonacci numbers defined by
$$
F^{(k)}_{i+1} \, = \, F^{(k)}_{i} \, + \, F^{(k)}_{i-1} \, + \, \dots \, + F^{(k)}_{i-k+1} \quad,
$$
with initial conditions
$$
F^{(k)}_1 \, = \,F^{(k)}_2 \, =\, \dots \,= \,F^{(k)}_k \, = \,1 \quad .
$$
The sequences $\{b^i\}$, can be defined as a solution of the first-order recurrence
$$
f(i)\,=\, b \,f(i-1) \quad, \quad f(0)=1 \quad,
$$
and $F_i$, and $F^{(k)}_i$ are specific examples of $C$-finite sequences.
Let's first formally define a $C$-finite sequence.
{\bf Definition}: A $C$-finite sequence $f(i)$ (no relation to the Fibonacci numbers) of order $L$, is
a sequence defined by a recurrence of the form,
$$
f(i)=c_1 f(i-1) \, + \, c_2 f(i-2) \, + \, \dots \, + \, c_L f(i-L) \quad,
$$
where $c_1, \dots c_L$ are {\bf constants}, subject to some initial conditions
$$
f(0)=d_0 \, , \, f(1)=d_1 \, , \, \dots \, , \, f(L-1)= d_{L-1} \quad ,
$$
for some (other) constants $d_0, \dots d_{L-1}$.
We denote a $C$-finite sequence by the pair $[[d_0, \dots, d_{L-1}] \, ; \, [c_1, \dots, c_L]]$ and in this paper
assume that the $d's$ and $c's$ are integers, so the sequence $\{f(i)\}_{i=0}^{\infty}$ is an {\bf integer sequence}.
So in this notation, the sequence $\{2^i\}_{i=0}^{\infty}$ featuring in the original definition of the Stern array is
denoted by $[[1],[2]]$ and more generally $\{b^i\}_{i=0}^{\infty}$ is denoted by $[[1],[b]]$, while
the Fibonacci sequence is $[[0,1],[1,1]]$ and the $k$-bonacci sequence is $[[0,1^{k-1}], [1^k]]$.
Before going on, we need an important {\it observation}.
{\bf Observation}: Given a $C$-finite sequence $f(i)$ of order $L$,
{\it any} linear combination of the form
$$
\sum_{i=M_1}^{M_2} a_i f(n+i) \quad,
$$
where $M_1