Calculating "factorial sequence" of a rational function

Question

Concise statement

For a reciprocal of a polynomial, $f = \frac{1}{p}$, I want to construct a sequence $(c_n)_{n=0}^\infty$ such that for all $N\ge 0$ $$f(k)k! = \sum_{n=0}^{N-1} c_n(k-n)! + O((k-N)!). $$ How can I rigorously calculate $(c_n)$? Bonus points for options which generalize to rational functions $f$, or can practically implemented as computer program.

Verbose statement + context

In a paper I am writing, I prove a lower bound of the form: $$ L_k\ge\left(1+\frac{1}{k}+\frac{1}{k(k-1)}+\frac{1}{(k-1)(k(k-1)(k-2)-k)}\right)k!+O(1).$$ In previous literature, lower bounds have always been expressed as a sum of factorials. (e.g. $L_k \ge k!+(k-1)!+(k-2)!+O(1)$ and $L_k \ge k!+(k-1)!+O(1)$ were the previously recorded bounds)

So obviously, I wished to express my result in this format, for comparison. Of course, that begs the question, for a reciprocal of a polynomial, $f = \frac{1}{p}$, can I construct a sequence $(c_n)_{n=0}^\infty$ such that for all $N\ge 0$ $$f(k)k! = \sum_{n=0}^{N-1} c_n(k-n)! + O((k-N)!). $$

Through some dirty approximations, I got the following sequence in my case$$c_0,c_1,\dots = 0,0,0,0,1,-2,7,-32,179,-1182,8993,-77440,744425,-7901410,\dots$$ where $f = \frac{1}{(k-1)(k(k-1)(k-2)-k)}$. Checking Laurent expansions at infinity with Wolfram Alpha, this seems to be correct. I saw $f -\sum_{n=0}^N c_n\frac{1}{\prod_{i=0}^n (k-n)}$ always had the right order of magnitude.

Nevertheless, my methods were quite unsound. I was using a Python program to do polynomial division on polynomials of the form $P_m(x):=\prod_{i=0}^m (x-i)$, and eyeballing convergence taking $m$ to be large. Rinse and repeat to get each term.

Presumably $(c_n)$ can be calculated in a better, more rigorous way. My question is how. Bonus points if the method easily extends to handle general rational functions, where $f = p/q$ with both the numerator and denominator being polynomials. I am also interested in implementing such a method into a program.

Answer 1

A comment on the issue of determining the coefficients of the expansion in function series in Fedor Petrov's answer. Recall the generating function of the Stirling numbers of the second kind $S(n,r)$ : $$\phi_r(x):=\frac{x^r}{(1-x)\dots(1-rx)}=\sum_{n\ge0}S(n,r)x^n,$$ and the "Stirling inversion" (i.e. the fact that the infinite triangular matrices of the Stirling numbers of the first and second kind are inverses of each other), that can be expressed as $$x^n=\sum_{n\ge0} s(r,n)\phi_r(x).$$

As a consequence, any formal $\phi_r$ series can be expanded in a formal power series and conversely,
$$\sum_{r\ge0} c_r\phi_r=\sum_{n\ge0}b_nx^n,$$ the correspondence being: $$ b_n=\sum_{r=0}^n S(n,r)c_r$$ $$ c_r=\sum_{n=0}^r s(r,n)b_n.$$

Warning. These notations differ slightly from your expansion, that is by a shift of indices and by a factor $x$, which can be easily adjusted. In fact, if you do not have strong objections to change notation, it could be worth to adapt it to that of Stirling numbers.

Answer 2

Dividing by $k!$ we write $$ f(k)=c_0+c_1\cdot \frac1k+c_2\cdot \frac1{k(k-1)}+\ldots+c_{N-1}\frac1{k(k-1)\ldots(k-N+2)}+O(k^{-N}). $$ Denoting $k=1/x$ this reads as $$ f(1/x)=c_0+c_1x+c_2\frac{x^2}{1-x}+c_3\frac{x^3}{(1-x)(1-2x)}+\ldots+c_{N-1}\frac{x^{N-1}}{(1-x)(1-2x)\ldots(1-(N-2)x)}+O(x^N), \quad x\to 0. $$ If $f(k)=p(k)/q(k)$ and $\deg p\leqslant \deg q$, we may expand $f(1/x)$ as a power series in $x$ and find $c_0,c_1,\ldots$ consequently equalizing the coefficients of $1,x,x^2,\ldots$ in both sides.

Answer 3

Let's assume that $p$ is monic of degree $d$ and has distinct roots $r_1,\dots,r_d$. Consider the partial fraction decomposition: $$f(x) = \frac{1}{p(x)} = \sum_{i=1}^d \frac{a_i}{x-r_i},$$ where $a_i:=\frac{1}{p'(r_i)}$.

Then $$\frac{1}{x}f(\frac{1}{x}) = \sum_{i=1}^d \frac{a_i}{1-r_ix}.$$ From answers of Fedor and Pietro, and the generating function for Stirling numbers of first kind, it follows that $$\sum_{n\geq 0} c_n \frac{z^n}{n!} = \sum_{i=1}^d a_i (1+z)^{r_i},$$ and thus $$c_n = n! \sum_{i=1}^d a_i \binom{r_i}n = \sum_{i=1}^d a_i (r_i)_n.$$

It is straightforward to further generalize this to the case of repeated roots and/or $f$ being a ratio of two polynomials.