Integrality of a sequence formed by sums Consider the following sequence defined as a sum
$$a_n=\sum_{k=0}^{n-1}\frac{3^{3n-3k-1}\,(7k+8)\,(3k+1)!}{2^{2n-2k}\,k!\,(2k+3)!}.$$

QUESTION. For $n\geq1$, is the sequence of rational numbers $a_n$ always integral?

 A: (More of a comment than an answer.) Maple's GuessRecurrence function finds that the following recurrence holds for at least 500 terms:
$$8(7n+8)(2n+5)(n+2)a_{n+2} - 6(252n^3+1233n^2+1930n+960)a_{n+1}
+81(7n+15)(3n+2)(3n+4)a_n=0.$$
I don't know if that would lead to an easy proof of integrality even if it can be established.
A: I've found a way around the problem that employs the Wilf-Zeilberger method. The challenging part was to come up with one side of the below "partner function" leading to an identity, as depicted in the process.
To this end, define these two discrete functions
$$F(n,k)=\frac{3^{k + 1}\binom{3n-k}{n-1-k}2^{2n}}{3^{3n}} \qquad \text{and}$$
$$G(n,k)=\frac{3^{k+1}\binom{3n-k+1}{n-k}2^{2n-1}(2k^2-18kn-19k-21n-24)}
{3^{3n+3}(2n^2+5n+3)}.$$
Next, one checks $F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k)$ and sum both sides over all $k=0$ through $k=n$, with $b_n=\sum_{k=0}^{n-1}F(n,k)$, so that the right-hand side telescopes in the manner
$$b_{n+1}-b_n=\sum_{k=0}^nG(n,k+1)-\sum_{k=0}^nG(n,k)=G(n,n+1)-G(n,0).$$
After some simplification, one finds
$$\frac{2^{2n}}{3^{3n}}\sum_{k=0}^{n-1}3^{k + 1}\binom{3n-k}{n-1-k}
=\sum_{k=0}^{n-1}\frac{2^{2k}(7k + 8)(3k + 1)!}{3^{3k+1}k!(2k+3)!}.$$
Or, equivalently, we arrive at
$$\sum_{k=0}^{n-1}\frac{3^{3n-3k-1}(7k + 8)\,(3k + 1)!}{2^{2n-2k}\,k!\,(2k+3)!}
=\sum_{k=0}^{n-1}3^{k + 1}\binom{3n-k}{n-1-k}$$
which proves the desired integrality.
UPDATE. Following the pointer by RobPratt to OEIS, we garner these identities:
\begin{align*}
\sum_{k=0}^{n-1}\frac{3^{3n-3k-1}(7k + 8)\,(3k + 1)!}{2^{2n-2k}\,k!\,(2k+3)!}
&=\sum_{k=0}^{n-1}3^{k + 1}\binom{3n-k}{n-1-k} \\
&=\frac32\sum_{k=0}^n\binom{3k}k\binom{3n-3k}{n-k}-\frac{3(3n+1)}{2(2n+1)}\binom{3n}n.\end{align*}
A: Let $A(x) = \sum_{n=1}^\infty a_n x^n$ and let
$$S(x) = \sum_{k=0}^\infty (7k+8)\frac{(3k+1)!}{k!\,(2k+3)!} x^k.$$
Then the formula for $a_n$ gives
$A(x) = R(x)S(x)$,
where
$$R(x) = \frac{1}{3}\biggl(\frac{1}{1-\frac{27}{4} x} -1\biggr).$$
A standard argument, for example by Lagrange inversion, gives
$$S\left(\frac{y}{(1+y)^3}\right)=\frac{4+y}{3(1+y)^2}.$$
A straightforward computation gives
$$R\left(\frac{y}{(1+y)^3}\right) = \frac{9y}{(4+y)(1-2y)^2}.$$
Thus
$$A\left(\frac{y}{(1+y)^3}\right)=\frac{3y}{(1-2y)^2(1+y)^2}.$$
Since the power series expansion of $y/(1+y)^3$ starts with $y$ and has integer coefficients, its compositional inverse has integer coefficients, so $A(x)$ does also.
A: Here is another proof, inspired by Tewodros Amdeberhan's. We represent the sum as a constant term in a power series.
To represent $(7k+8) \frac{(3k+1)!}{k!\,(2k+3)!}$ as a constant term, we need to express it as a linear combination of binomial coefficients. To do this we express $7k+8$ as a linear combination of $(2k+2)(2k+3)$, $k(2k+3)$, and $k(k-1)$ since each of these polynomials yields a binomial coefficient when multiplied by $\frac{(3k+1)!}{k!\,(2k+3)!}$. We find that
$$(7k+8) \frac{(3k+1)!}{k!\,(2k+3)!}=\frac13 \left[4\binom{3k+1}{k} -7\binom{3k+1}{k-2} -2\binom{3k+1}{k-2}\right].$$
Using $\binom{n}{j} =\text{CT}\, (1+x)^n/x^j$, where CT denotes the constant term in $x$, we have
$$
\begin{aligned}
(7k+8) \frac{(3k+1)!}{k!\,(2k+3)!} &= \text{CT}\, \frac13\left( 4\frac{(1+x)^{3k+1}}{x^k}
   -7\frac{(1+x)^{3k+1}}{x^{k-1}}-2\frac{(1+x)^{3k+1}}{x^{k-2}}\right)\\
   &=\text{CT}\,\frac{(1-2x)(x+4)(1+x)^{3k+1}}{3x^k}.
 \end{aligned}
$$
Multiplying by $3^{3n-3k-1}/2^{2n-2k}$, summing on $k$ from 0 to $n-1$, and simplifying gives
$$
a_n = 3\,\text{CT}\, \frac{(1+x)^{3n+1}}{x^{n-1}(1-2x)} -3\left(\frac{27}{4}\right)^n\!\text{CT}\,\frac{x(1+x)}{1-2x}.
$$
The second constant term is 0, so $a_n = 3\,[x^{n-1}]\, (1+x)^{3n+1}/(1-2x)$, which is clearly an integer.
In fact this gives the formula
$$a_n = 3 \sum_{k=0}^{n-1}2^{n-k-1}\binom{3n+1}{k}.$$
Tewodros's formula can also be derived in the same way; if we represent his sum as a constant term and simplify we also get
$3\,[x^{n-1}]\, (1+x)^{3n+1}/(1-2x)$.
