Solution 1:

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.

Solution 2:

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}