Constant term in $\cot^n(x)$ The Laurent expansion for $\cot x$ is
$$\cot x =\frac{1}{x} +\sum_{k=1}^\infty \frac{(-1)^k 4^k B_{2k}}{(2k)!}x^{2k-1}$$
where $B_k$ are the Bernoulli numbers.
I would like to have a formula for the constant term of $\cot^n x$. I mean something more explicit that writing down the sum of all products of $n$ terms whose exponents add to zero. From numerical experimentation, it seems that the numerators are sequence A216254 and the denominators are A225149 but I have no proof and in any case there are no formulas or recurrence relations on OEIS for those sequences.
 A: This follows from Cauchy's formula after the change of variables $
z = \arctan t$:
$$
\frac{1}{{2\pi i}}\oint_{(0 + )} {\cot ^n z\frac{{dz}}{z}}  = \frac{1}{{2\pi i}}\oint_{(0 + )} {\frac{1}{{\tan ^n z}}\frac{{dz}}{z}}  = \frac{1}{{2\pi i}}\oint_{(0 + )} {\frac{t}{{(t^2  + 1)\arctan t}}\frac{{dt}}{{t^{n+1} }}} .
$$
Thus, the constant term in the Laurent expansion of $\cot^n z$ is precisely the $n$th coefficient in the Maclaurin expansion of $$
\frac{z}{{(z^2  + 1)\arctan z}}.
$$
Now if we write
$$
\frac{z}{{(z^2  + 1)\arctan z}} = \sum\limits_{n = 0}^\infty  {a_n z^{2n} } 
$$
and multiply both sides by
$$
\frac{{\arctan z}}{z} = \sum\limits_{n = 0}^\infty  {\frac{{( - 1)^{n} }}{{2n + 1}}z^{2n} },
$$
use the geometric series on the left-hand side, and compare like powers of $z$, we obtain the recurrence relation
$$
( - 1)^n  = \sum\limits_{k = 0}^n {\frac{{( - 1)^k }}{{2k + 1}}a_{n - k} }  \Longleftrightarrow a_n  = ( - 1)^n  - \sum\limits_{k = 1}^n {\frac{{( - 1)^k }}{{2k + 1}}a_{n - k} } 
$$
with $a_0=1$.
To obtain an asymptotic formula for the coefficients $a_n$, we can proceed as follows. The singularities of the generating function are located at $z=\pm i$ and
$$
\frac{z}{{(z^2  + 1)\arctan z}} \sim  - \frac{1}{{(1 \pm iz)\log (1 \pm iz)}}
$$
as $z\to \pm i$. Using Theorem VI.2 in Analytic Combinatorics by Flajolet
and Sedgewick, we derive
$$
a_n  \sim ( - 1)^n \frac{2}{{\log (2n + 1)}}\left( {1 - \frac{\gamma }{{\log (2n + 1)}} - \frac{{\pi ^2  - 6\gamma ^2 }}{{6\log^2 (2n + 1)}} +  \cdots } \right)
$$
as $n\to +\infty$, where $\gamma$ is the Euler–Mascheroni constant.
