I have a sequence of polynomials $Q_k(x, y)$, $k\geq 1$ defined recursively as follows:

  1. $Q_1=x$.

  2. There is a sequence of polynomials $p_j(y)$ of degree $j$ such that $Q_{2m}$ is of the form \begin{eqnarray}\frac{p_{0}(y)}{(2m)!}x^{2m}+\frac{p_{1}(y)}{(2m-2)!}x^{2m-2}+\cdots+\frac{p_{m-1}(y)}{2!}x^2+p_{m}(y)\end{eqnarray} and $Q_{2m+1}$ \begin{eqnarray} \frac{p_{0}(y)}{(2m+1)!}x^{2m+1}+\frac{p_{1}(y)}{(2m-1)!}x^{2m-1}+\cdots+\frac{p_{m-1}(y)}{6}x^3+p_{m}(y)x \end{eqnarray}

  3. $Q_k(k+1-2i, k+1)$ as a polynomial in $i$ has roots $1, 2, \cdots, k$.

Find a general formula of $Q_k(x, y)$. The following are the first 4 polynomials in the sequence: \begin{align*} Q_1&=x\\ Q_2&=\frac{x^2}{2}-\frac{y}{6}\\ Q_3&=\frac{x^3-xy}{6}\\ Q_4&=\frac{x^4-2x^2y}{24}+\frac{y(5y+2)}{360} \end{align*}

In fact, it suffices to find the sequence of polynomials $p_j(y)$. Here's what I've got so far. By condition (2) and (3), we have \begin{eqnarray} \sum_{j=0}^m \frac{p_j(2m+1)}{(2m-2j)!}(2i-2m-1)^{2m-2j}=\frac{2^{2m}}{(2m)!}(i-1)(i-2)\cdots(i-2m) \end{eqnarray} If we let $m=j+k$, differentiate both sides with respect to $i$ $2k$ times and put $\displaystyle i=\frac{2j+2k+1}{2}$, we have \begin{eqnarray} p_j(2j+2k+1)=\left.\frac{d^{2k}}{di^{2k}}\right|_{i=\frac{2j+2k+1}{2}}\frac{2^{2j}}{(2(j+k))!}(i-1)(i-2)\cdots(i-2(j+k)) \end{eqnarray} Letting $k=0, 1, \cdots, j$, we get the $j+1$ values taken by $p_j$ at $2j+1, 2j+3, \cdots, 4j+1$. $p_j$ then can be computed using Lagrangian interpolation.

It seems to me that, though the above algorithm can be implemented on a computer to get a few polynomials in the sequence, it does not yield directly a general formula I want. Is there another better way to go about getting a general formula?

Edit: The first 6 members in the polynomial sequence $p_j(y)$ are the following: \begin{align*} p_0(y)&=1\\ p_1(y)&=-\frac{y}{6}\\ p_2(y)&=\frac{y(5y+2)}{360}\\ p_3(y)&=-\frac{y(35y^2+42y+16)}{45360}\\ p_4(y)&=\frac{y(5y+4)(35y^2+56y+36)}{5443200}\\ p_5(y)&=-\frac{y(385y^4+1540y^3+2684y^2+2288y+768)}{359251200} \end{align*}

It is noteworthy that the denominators happen to be the first 6 numbers in this sequence, as pointed out by Solomonoff's Secret.

  • 649
  • 4
  • 19
  • Are you sure that $\frac {\partial Q_{2m}} {\partial x} = Q_{2m-1}$? From your formulae one can only deduce $\frac {\partial Q_{2m+1}} {\partial x} = Q_{2m}$. – Alex M. Jul 18 '15 at 12:47
  • That $\frac{\partial Q_{2m}}{\partial x}=Q_{2m-1}$ is another condition, not a consequence of the form of the polynomial. I use the word 'Besides'. – No_way Jul 19 '15 at 06:16
  • As originally written, the question suggested that the polynomials $p$ are the same for all the polynomials $Q$. In reality, they are different for each polynomial $Q$, therefore I have added a supplementary index to take this into account. – Alex M. Jul 19 '15 at 09:00
  • 1
    @No_way Is there a motivation for having condition (3) hold only for even indices? When you make a distinction between $p_{2m+1,m}(y) $ and $p_{2m,m}(y)$ then we no longer have (the implicit) $\frac{\partial}{\partial x}Q_{2m+1} = Q_{2m}$ – will Jul 20 '15 at 17:59
  • 2
    @Alex M. The polynomials $p$ for $Q_{2m}$ and $Q_{2m+1}$ are the same! You totally misunderstood what I meant, and should have consulted me before making amendment, and I find the authoratitative tone of your comment very deplorable. – No_way Jul 21 '15 at 13:24
  • 2
    @will In fact I didn't make the distinction between those polynomials originally. It is this guy Alex M. who made the modification without my consent. Please see the current version. – No_way Jul 21 '15 at 13:26
  • 1
    It's actually non-trivial that (2), (3), and (4) produce coefficients of $m$ powers $x$ that can be expressed with just $m$ polynomials in $y$. (Notice $(2m+2-2k)p_{m+1,k}(y) = p_{m,k}(y).)\ $ The fourth condition will probably yield (a combination of) the first kind of stirling numbers, that will probably be polynomials of $m$ of degree $k$. The motivation for this recurrence probably has an unnecessary parity. – will Jul 21 '15 at 15:40
  • 2
    @will I have give more clarification and information to the original question. Please take a look. – No_way Jul 24 '15 at 16:16
  • 1
    This may be relevant: https://oeis.org/search?q=1%2C+6%2C+360%2C+45360%2C+5443200&language=english&go=Search – Reinstate Monica Jul 24 '15 at 16:46
  • @Solomonoff'sSecret Thank you! This may be helpful. I thought that the denominators of the polynomials in the sequence are some random numbers... – No_way Jul 24 '15 at 16:56
  • 1
    Perhaps generating series can be used. For example if we define $f(t,y) = \sum_{n=0}^\infty p_n(y)t^n$ and $g(t,x) = \sum_{n=0}^\infty \frac{x^{2n}}{(2n)!}t^n = \cosh(x\sqrt{t})$ then it follows that $f(t,y)g(t,x) = \sum_{n=0}^\infty Q_{2n}(x,y) t^n$. Don't know if this is useful or how one would proceed from this. – Winther Jul 24 '15 at 17:06
  • 1
    Another relation that might be useful for this problem is $Q_n(x+t,y) = \sum_{k=0}^n \frac{x^{n-k}}{(n-k)!} Q_{k}(t,y)$. This follows from a Taylor expansion using $\partial_x Q_{k}(x,y)=Q_{k-1}(x,y)$. – Winther Jul 24 '15 at 18:24
  • @Winther Thank you for pointing out this generating series approach. I'll look into it. – No_way Jul 25 '15 at 07:32

2 Answers2


We are essentially looking for the coefficients of falling factorials. I am not a fan of the first kind of Stirling numbers, because they are (irreversibly?) recursive. (See the more explicit Bernoulli numbers for comparison). Never the less, begin with condition (3): $$ 1 \le m\implies\ Q_m(m+1-2i,m+1) = (-2)^m\binom{i-1}{m} $$ (Corrected by San). When we use the translation $z=m+1-2i,\ $ we can apply the structure of condition (2): $$ Q_m(z,m+1) = \sum_{k=0}^{\lfloor\frac{m}{2}\rfloor}p_k(m+1)\frac{z^{m-2k}}{(m-2k)!} = (-2)^m\binom{\frac{m-1-z}{2}}{m} $$ Notice the top of the binomial expands to $(z+1-m)(z+3-m)\cdots (z-3+m)(z-1+m),\ $ and so the parity of the powers of $z$ match the parity of $m$. This reassures us that condition (2) is well defined (for all complex $z$) even though the sum omits half of the powers of $z$. We also require roots of unity, the exact choice of roots is irrelevant as long as it is consistent: $\mu_N \in\mathbb{C}$ is a $N$th root unity when $k\in N\mathbb{Z}\Leftrightarrow \mu_N^k = 1.\ $ When $m$ is odd and $m < 2(m-2N),\ \ (m-2N)\mathbb{Z} \cap \{m,m-2,\ldots,3,1\} = \{m-2N\};\ $ which makes the discrete fourier transform of $z^{m-2k}$ in $Q_m$ useful: $$\mbox{when } m \mbox{ is odd, and } \mu_B := \exp(\frac{2\pi i}{B}): $$ $$ 4N < m\implies\ p_N(m+1)\frac{m-2N}{(m-2N)!} = (-2)^m\sum_{j=1}^{m-2N} \binom{(m-1-\mu_{m-2N}^j)/2}{m} $$ We do not see an obvious way to show that $p_N$ is a degree $N$ polynomial, other than to tediously expand the powers of $z$ in the right hand binomial. We assume that everyone is convinced that $p_N$ is polynomial. Now we too resort to Lagrange interpolation: $$\mbox{using }\ m+1=4N+2+2h,\ \ h=0\ldots N$$ $$ p_N(y) = 4^{2N+1}(-1)^{N+1}(N+1)\binom{\frac{y}{2}-2N-1}{N+1} $$$$ \sum_{h=0}^N \sum_{j=1}^{2N+1+2h} \frac{(2N+2h)!(-4)^h}{y-4N-2-2h}\binom{N}{h} \binom{2N+h-\frac{\mu_{2N+1+2h}^j}{2}}{4N+1+2h} $$

This is a finite formula, that takes no limits. It is not suitable for computation. (The case $N=4$ evaluates a degree 25 binomial). Now that Winther corrected the fourier transform, we can numerically evaluate the first few polynomials. Here is the output of a test program (by gcc):

This machine uses 128 bit long doubles:

Will claims 'No_way's polynomial, p1(y), is about:
-1.0000y^1   -0.0000y^0
over 1!6^1  =  6
with relative error > 7.12e-16   making coefficent errors > 4.27e-15

Will claims 'No_way's polynomial, p2(y), is about:
+1.0000y^2   +0.4000y^1   -0.0000y^0
over 2!6^2  =  72
with relative error > 3.20e-14   making coefficent errors > 2.31e-12

Will claims 'No_way's polynomial, p3(y), is about:
-1.0000y^3   -1.2001y^2   -0.4561y^1   -0.0053y^0
over 3!6^3  =  1296
with relative error > 7.44e-11   making coefficent errors > 1.16e-07

Will claims 'No_way's polynomial, p4(y), is about:
+0.8094y^4   +18.4170y^3   -500.4859y^2   +6988.6559y^1   -36277.1494y^0
over 4!6^4  =  31104
with relative error > 1.75e-06   making coefficent errors > 1.97e+03

The alternative recurrence relation to this "formula" would be to interpolate the correct combination of stirling numbers. If you construct a generating function, you might need to integrate something like $\ln(z)/(z^2+y^2)$, but this is a wild guess from past experience.

  • 636
  • 3
  • 6
  • 1
    I want to point out that in the first displayed equation in your answer, you missed the scalar multiple $\displaystyle \frac{2^m}{m!}$ in the RHS, which is the leading coefficient of $Q_m(m+1-2i, m+1)$. I am not sure if this will affect substantially your subsequent formulae. – No_way Jul 25 '15 at 07:30
  • I meant 'the leading coefficient of $Q_m(m+1-2i, m+1)$ as a polynomial in $i$'. – No_way Jul 25 '15 at 08:03
  • The $m!$ is provided by the binomial. I definitely missed the power of two through. – will Jul 25 '15 at 11:11
  • Thank you for your answer. I will look into it. Actually I would like to see formula involving Stirling numbers of the second kind, as I somehow saw them lurking in my computation. Though you are not a fan of Stirling number of the first kind, could you tell if it is possible to writing down a general formula in terms of the second kind? – No_way Jul 26 '15 at 17:06
  • 1
    It is definitely possible to use the [first kind](http://mathworld.wolfram.com/StirlingNumberoftheFirstKind.html) of Stirling numbers. You found the points of $p(y)$ by taking derivatives of the falling factorial, $4^j(i-1)(i-2)\cdots(i-2(j+k)),\ $ which happens to be the row generating function of the 1st stirlings. See Weisstein's (8) equation. – will Jul 26 '15 at 21:37
  • 1
    Working out a recurrence among just the $p_k(y)$, and thus moving past the parity of bivarite the $Q_m(x,y)$ recurrence, would make much life easier. – will Jul 26 '15 at 21:41
  • 2
    I tried to plot the formula for $p_N(y)$ and the results looks very strange. For example $p_2(y) \approx 800.0(x-11.0)^2$ which does not match well with the formulas given by OP. There is definately something wrong with the prefactors since $p_3(0) > 10^9$ and $p_4(0) > 10^{15}$. – Winther Jul 27 '15 at 17:57
  • @will Could you please address Winther's concern? I concur with his estimate of $p_2$ using your formula. – No_way Jul 29 '15 at 14:24
  • @Winther helpfully found an error in my discrete fourier transform. The indices on my inner sum originally ran from $j=0\ldots m-2N,\ $ which samples $\mu_{m-2N}^0 = \mu_{m-2N}^{m-2N}\ $ twice. I have corrected my answer. Is there a place to post my code (in c) that verifies the formula for $p_0\ldots p_4$? – will Jul 29 '15 at 19:52
  • 1
    My appologies. I had a bug in the code. I can confirm that your new formula gives the correct answer for $n=1,2,3,4$ using Mathematica (symbolic evaluation of the sums). – Winther Jul 29 '15 at 20:30
  • My apologies for writing a buggy answer. – will Jul 29 '15 at 20:34

Using Wolfram Alpha, I find that the first 6 members $p_j(x)$, $0\leq j\leq 5$, of the polynomial sequence happen to be the first 6 non-zero coefficients of the Maclaurin series of \begin{eqnarray}\left(\frac{t}{\sinh(t)}\right)^x.\end{eqnarray} I conjecture that \begin{eqnarray}\left(\frac{t}{\sinh(t)}\right)^x=\sum_{j=0}^\infty p_j(x)t^{2j}.\end{eqnarray}

  • 649
  • 4
  • 19
  • For integer $m$ we should have $\oint_{|z|=1}\frac{e^{xz}dz}{\sinh^{m+1}(z)}\ =\ (-2)^m\binom{(m-1-x)/2}{m}$ so that the coefficients of $p_k(m+1)\frac{x^{m-2k}}{(m-2k)!}$ produce the roots in $Q(x,m+1)$ – will Sep 23 '15 at 01:22