A closed form for an integral expressed as a finite series of $\zeta(2k+1)$, $\pi^m$ and a rational? In this paper the following beautiful integral expression for $\zeta(3)$ is derived: 
$$\zeta(3)=\frac{1}{7}\,\int_0^{\pi} x\,(\pi-x)\csc(x)\, dx$$
In a comment at the end of this question, I mentioned the following tweak: 
$$f(n):= \int_0^{\pi} x^n\,(\pi-x)\csc(x)\,dx$$
Maple and Mathematica both return closed forms for $f(n)$ as a finite series of $\zeta(2k+1)$, weighted by $π^m$ and a rational ($n,k,m \in \mathbb{N}$). This is what I get for $n=1..12$ : 

I then tried to find a closed form for the above and got as far as:
For $n=$ odd:
$$f(n)=\scriptsize \left(\Gamma(n+1)\,\sum_{k=1}^{n}(-1)^{k+1}\,\frac{(4-2^{1-2k})\,k\,\pi^{n+1-2k}}{\Gamma(2+n-2k)}\,\zeta(2k+1)\right)-\left(\frac{i}{2}\right)^{n+1}\,(2^{n+2}-1)\,\Gamma(n+2)\,\zeta(n+2)$$
For $n=$ even:
$$f(n)=\scriptsize \left(\Gamma(n+1)\,\sum_{k=1}^{n}(-1)^{k+1}\,\frac{(4-2^{1-2k})\,k\,\pi^{n+1-2k}}{\Gamma(2+n-2k)}\,\zeta(2k+1)\right)+\left(\frac{i}{2}\right)^{n}\,(2^{n+1}-1)\,\pi\,\Gamma(n+1)\,\zeta(n+1)$$
Note that only the last part differs, which encourages me to believe that a simpler formula might exist. Could this be simplified? Since the CAS-tools evaluate the integral very quickly into the series, I am also keen to better understand which mathematical steps they take.
Thanks. 
 A: Not an answer, but... Mathematica evaluates the indefinite integral
$$\int x^k \csc x d x$$ in terms of polylogs, from which your integral follows readily, so the obvious question is - how does it do the indefinite integral. The short answer is: The Risch algorithm (or, more precisely, an extension, since the original algorithm only deals with elementary antiderivatives). However, if you want to do it by hand, the two methods which come to mind are:


*

*Make the substitution $u = \exp(i x),$ which transforms your integrand into a power of log times a rational function - this is plausbible, since the arguments of the polylogs have the form $\exp(i x),$ so Mathematica might, indeed, be doing this.

*Expand $\csc x$ in its Laurent series around $0,$ and integrate term by term. The Laurent series is explicit (in terms of Bernoulli numbers), so it may be that a simple manipulation gets it to your linear combination of zetas (I haven't tried).

A: I think the following articles can give a clue:
http://www.tandfonline.com/doi/abs/10.1080/10652460701688125?journalCode=gitr20 (Closed-form evaluation of some families of cotangent and cosecant integrals, by D. Cvijović) and http://www.sciencedirect.com/science/article/pii/S0377042702003588 (Integral representations of the Riemann zeta function for odd-integer arguments, by D. Cvijović and J. Klinowski).
Update. In particular, using $2x^n=E_n(x)+E_n(1+x)$ and $E_n(1+x)=\sum\limits_{k=0}^n\binom{n}{k}E_k(x)$, we get
$$f(n)=\frac{\pi^{n+2}}{2}\left[I_n+\sum\limits_{k=0}^n\binom{n}{k}I_k-I_{n+1}-\sum\limits_{k=0}^{n+1}\binom{n+1}{k}I_k\right]=\frac{\pi^{n+2}}{2}\left[I_n-2I_{n+1}-\sum\limits_{k=1}^n\binom{n}{k-1}I_k\right], \tag{1}$$
where $$I_n=\int\limits_0^1E_n(x)\csc{(\pi x)}\,dx.$$
Odd indices in (1) doesn't contribute because $E_n(1-x)=(-1)^nE_n(x)$ leads to the integrand which is antisymmetric with regard to $x\to 1-x$. For even indices, $I_n$ was calculated in the second paper indicated above as $$I_{2n}=(-1)^n\frac{(4-2^{1-2n})(2n)!}{\pi^{2n+1}}\zeta(2n+1),$$ and we obtain, for even $n=2m$:
$$f(2m)=(-1)^m\frac{\pi}{2}(4-2^{1-2m})(2m)!\,\zeta(2m+1)- \sum\limits_{p=1}^m(-1)^p\;\frac{\pi^{2(m-p)+1}}{2}\;\binom{2m}{2p-1}(4-2^{1-2p})(2p)!\,\zeta(2p+1),$$
and for odd $n=2m-1, m>1$:
$$f(2m-1)=(-1)^{m+1}(4-2^{1-2m})(2m)!\,\zeta(2m+1)-\sum\limits_{p=1}^{m-1}(-1)^p\;\frac{\pi^{2(m-p)}}{2}\;\binom{2m-1}{2p-1}(4-2^{1-2p})(2p)!\,\zeta(2p+1).$$
