El uso de la reflexión de la fórmula seguida por la repetición de la fórmula y el Beta representación integral (DLMF)
\begin{align}
\frac{x(1-x)}{\sin \pi x}&=\frac{1}{\pi}x(1-x)\Gamma(x)\Gamma(1-x)\\
&=\frac{1}{\pi}\Gamma(x+1)\Gamma(2-x)\\
&=\frac{2}{\pi}B(1+x,2-x)\\
&=\frac{2}{\pi}\int_0^\infty \frac{t^x}{(1+t)^3}\,dt\\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{e^{u(2x-1)}}{\cosh^3u}\,du
\end{align}
la última expresión se obtiene con $t=e^{2u}$. Entonces, para el cálculo de
\begin{equation}
I_f=\int_0^1\frac{x(1-x)}{\sin \pi x}f(x)\,dx
\end{equation}
uno puede expresar
\begin{align}
I_f&=\frac{1}{2\pi}\int_0^1\int_{-\infty}^{\infty}\frac{e^{u(2x-1)}}{\cosh^3u}\,duf(x)\,dx\\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{e^{-u}}{\cosh^3u}\,du\int_0^1e^{2ux}f(x)\,dx\\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{e^{-u}}{\cosh^3u}F(u)\,du \label{eq:intf}
\end{align}
suponiendo que el cambio de la integración de la orden es válida y que denota
\begin{equation}
F(u)=\int_0^1e^{2ux}f(x)\,dx
\end{equation}
Si $F(u)$ es analítica en el semiplano $\Im(u)>0$ e $\left|F(u)\right|=o\left(\frac{e^{4u}}{u} \right)$ para $\left|u\right|\to \infty$,$I_f$ se evalúa mediante la integración a lo largo del eje real cerrada por el semi-círculo grande $\Im(u)>0$, utilizando el residuo método. Los polos están situados en $u_n=i(2n+1)\pi/2$ con $n=0,1,2...$. Los residuos se $1/2F''(i(2n+1)\pi/2)-F'(i(2n+1)\pi/2)$ donde $F''(z)$ e $F'(z)$ son, respectivamente, la primera y la segunda derivada de $F(z)$. Como el semi-círculo contribución se desvanece, se trata de
\begin{equation}
I_f=i\sum_{n=0}^\infty\left[\frac{1}{2}F''(i(2n+1)\frac{\pi}{2})-F'(i(2n+1)\frac{\pi}{2})\right]
\end{equation}
Al $f(x)=1$, para expresar la integral original, $F(u)=\frac{e^{2u}-1}{2u}$, un cálculo simple muestra, como era de esperar, que
\begin{equation}
I=\frac{8}{\pi^3}\sum_{n=0}^\infty \frac{1}{(2n+1)^3}=\frac{7\zeta(3)}{\pi^3}
\end{equation}
Otra expresión para el resultado se obtiene por derivación bajo la integral
\begin{equation}
I_f=-2i\sum_{n=0}^\infty \int_0^1x(1-x)f(x)e^{i(2n+1)\pi x}\,dx
\end{equation}
Para una función real $f$, ya que la suma debe ser real, es suficiente para mantener el imaginario contribución a la integral:
\begin{equation}
I_f=2\sum_{n=0}^\infty \int_0^1x(1-x)f(x)\sin\left( (2n+1)\pi x \right)\,dx
\end{equation}
Ahora, supongamos que una función $f_p$ es conocido que las integrales
\begin{equation}
J_p=\int_0^1x(1-x)f_p(x)\sin\left( (2n+1)\pi x \right)\,dx=\frac{A_p}{(2n+1)^{2p+1}}
\end{equation}
que da la relación
\begin{equation}
I_{f_p}=2A_p\sum_{n=0}^\infty\frac{1}{(2n+1)^{2p+1}}=2A_p\left( 1-2^{-2p-1} \right)\zeta(2p+1)
\end{equation}
Denota la función de $Q^0(x)=x(1-x)f_p(x)$ e $Q^1,Q^2(x)$ la primera y la segunda antiderivada. Dos sucesivas integraciones por partes se puede realizar:
\begin{align}
J_p&=-(2n+1)\pi \int_0^1Q^1(x)\cos\left( (2n+1)\pi x \right)\,dx\\
&=(2n+1)\pi\left[Q^2(1)+Q^2(0)\right]-(2n+1)^2\pi^2\int_0^1Q^2(x)\sin\left( (2n+1)\pi x \right)\,dx \label{eq:jp}
\end{align}
Los parámetros libres en $Q^2(x)$ puede ser elegido con el fin de que $Q^2(1)=Q^2(0)=0$. Con
\begin{equation}
Q^2(z)=\int_0^zdt\int_0^tQ^0(u)\,du+az+b
\end{equation}
uno puede optar $b=0$ e $a=-\int_0^1\,dt\int_0^tQ^0(u)\,du$. Así
\begin{equation}
Q^2(z)=\int_0^zdt\int_0^tQ^0(u)\,du-z\int_0^1\,dt\int_0^tQ^0(u)\,du
\end{equation}
Si $f_p(x)$ es un polinomio, entonces $Q^2(x)$ también. Por construcción, $x=0$ e $x=1$ están entre sus raíces. Puede ser escrita como
\begin{equation}
Q^2(x)=x(1-x)f_{p+1}(x)
\end{equation}
o
\begin{equation}
f_{p+1}(x)=\frac{\int_0^xdt\int_0^tu(1-u)f_p(u)\,du-x\int_0^1\,dt\int_0^tu(1-u)f_p(u)\,du}{x(1-x)}
\end{equation}
$J_p$ puede ser escrito como
\begin{equation}
J_p=-(2n+1)^2\pi^2\int_0^1x(1-x)f_{p+1}(x)\sin\left( (2n+1)\pi x \right)\,dx
\end{equation}
Uno obtener
\begin{equation}
\int_0^1x(1-x)f_{p+1}(x)\sin\left( (2n+1)\pi x \right)\,dx=-\frac{1}{\pi^2}\frac{A_p}{(2n+1)^{2p+3}}
\end{equation}
y así
\begin{equation}
I_{f_{p+1}}=2A_{p+1}\left( 1-2^{-2p-3} \right)\zeta(2p+3)
\end{equation}
con
\begin{equation}
A_{p+1}=-\frac{A_p}{\pi^2}
\end{equation}
A partir de $f_1(x)=1$ se obtiene
\begin{align}
f_2(x)&=\frac{1}{12}(x^2-x-1)\\
f_3(x)&=\frac{1}{360}(x^4-2x^3-2x^2+3x+3)\\
f_4(x)&=\frac{1}{20160}(x^6-3x^5-3x^4+11x^3+11x^2-17x-17)\\
...
\end{align}
lo que da
\begin{align}
\int_0^1\frac{x(1-x)}{\sin\pi x}f_2(x)\,dx&=-\frac{31}{4}\frac{\zeta(5)}{\pi^5}\\
\int_0^1\frac{x(1-x)}{\sin\pi x}f_3(x)\,dx&=\frac{127}{16}\frac{\zeta(7)}{\pi^7}\\
\int_0^1\frac{x(1-x)}{\sin\pi x}f_4(x)\,dx&=-\frac{511}{64}\frac{\zeta(9)}{\pi^9}\\
...
\end{align}
A partir de $f_1(x)=x(3-x)$, otros de la serie se puede obtener. Por ejemplo
\begin{align}
&f_2(x)=\frac{1}{60}(2x^4-10x^3+5x^2+5x+5)\\
&f_3(x)=-\frac{1}{5040}(3x^6-21x^2(x^3-x^2-x-1)-49(x+1))\\
&\int_0^1\frac{x(1-x)}{\sin\pi x}f_2(x)\,dx=-\frac{381}{4}\frac{\zeta(7)}{\pi^7}\\
&\int_0^1\frac{x(1-x)}{\sin\pi x}f_3(x)\,dx=-\frac{1533}{16}\frac{\zeta(9)}{\pi^9}
\end{align}
Otros puntos de partida pueden ser obtenidos mediante la elección de otros miembros de la lista propuesta en la pregunta anterior. Por ejemplo, a partir de $f_1(x)=x^2P_3(x)$ anterior conduce a una aparente expresión diferente para $\zeta(9)$:
\begin{align}
&f_2(x)=-\frac{1}{56}(x^6-27x^5+57x^4-13x^3-13x^2-13x-13)\\
&\int_0^1\frac{x(1-x)}{\sin\pi x}f_2(x)\,dx=-\frac{22995}{8}\frac{\zeta(9)}{\pi^9}
\end{align}
Obtenido polinomios no son de la forma $x^pP_p(x)$ como se analiza en la pregunta, sin embargo el método anterior, tal vez puede ser adaptado en este caso.
EDITAR 04/06/2017 : (perdón por la longitud de esta respuesta...)
Uno puede caracterizar más precisamente la familia de estos polinomios.
Ayuda a symmetrize las expresiones:
\begin{equation}
I_f=\int_0^1\frac{x(1-x)}{\sin \pi x}f(x)\,dx=\frac{1}{8}\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}g(y)\,dy
\end{equation}
con $g(y)=f(x)$ e $x=(1+y)/2$. En esta forma es claro que extraño contribución del polinomio $g(y)$ se desvanece. El mismo simetrización de la propuesta de la descomposición de arriba se lee:
\begin{equation}
I_f=\frac{(-1)^n}{4}\sum_{n=0}^\infty \int_{-1}^1(1-y^2)g(y)\cos\left( (2n+1) \frac{\pi y}{2} \right)\,dy
\end{equation}
Uno puede adaptar el método desarrollado anteriormente. Si $g_p(y)$ es incluso un polinomio como
\begin{equation}
\sum_{n=0}^\infty \int_{-1}^1(1-y^2)g_p(y)\cos\left( (2n+1) \frac{\pi y}{2} \right)\,dy=\frac{A_p}{(2n+1)^{2p+3}}
\end{equation}
luego, mediante la integración de dos veces por una parte, el polinomio
\begin{equation}
g_{p+1}(y)=\frac{\int_{-1}^ydt\int_{-1}^t(1-u^2)g_p(u)\,du-\frac{y+1}{2}\int_{-1}^1\,dt\int_{-1}^t(1-u^2)g_p(u)\,du}{1-y^2}
\end{equation}
es tal que
\begin{equation}
\sum_{n=0}^\infty \int_{-1}^1(1-y^2)g_{p+1}(y)\cos\left( (2n+1)\frac{\pi y}{2} \right)\,dy=\frac{A_{p+1}}{(2n+1)^{2p+5}}
\end{equation}
con $A_{p+1}=-4A_p/\pi^2$.
así
\begin{equation}
\frac{1}{8}\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}g_{p+1}(y)\,dy=2A_{p+1}\left( 1-2^{-2p-3} \right)\zeta(2p+5)
\end{equation}
Uno puede mostrar que $g(y)$ es incluso un polinomio de $y$. Para $g_0(y)=1$ uno tiene, como se esperaba
\begin{equation}
\frac{1}{8}\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}\,dy= \frac{7\zeta(3)}{\pi^3}
\end{equation}
Entonces, la recurrencia anterior produce una serie de incluso polinomios $g_p(y)$ grado $2p$ dando sucesivas integral expresiones para $\zeta(2p+3)$. Debido a la paridad de observación, se puede concluir que cualquier polinomio $Q(y)$, incluso con su coeficiente de potencia idéntica a la de $g_p(y)$, es tal, que
\begin{equation}
\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}Q(y)\,dy=\left( -1 \right)^p8\left( 2^{2p+3}-1 \right)\frac{\zeta(2p+3)}{\pi^{2p+3}}
\end{equation}
La condición lee
\begin{equation}
Q(y)+Q(-y)=2g_p(y)
\end{equation}
La primera polinomios (escrito con $Y=y^2$) son:
\begin{align}
g_0(y)&=1\\
g_1(y)&=\frac{1}{12}\left( Y-5 \right)\\
g_2(y)&=\frac{1}{360}\left( Y^2-14Y+61 \right)\\
g_3(y)&=\frac{1}{20160}\left( Y^3-27Y^2+323Y-1385 \right)\\
g_4(y)&=\frac{1}{1814400}\left( Y^4-44Y^3+1006Y^2-11804Y+50521 \right)\\
g_5(y)&=\frac{1}{239500800}\left( Y^{5}-65Y^4+2410Y^3-53954Y^2+631621Y-2702765\right)\\
g_6(y)&=\frac{1}{43589145600}\left(Y^6-90Y^5+4915Y^4-178268Y^3+3980887Y^2-46590634Y+199360981 \right)
\end{align}
En términos de la no-simétrico de la función, cualquier polinomio de la forma
\begin{equation}
f(x)=g_p\left( 2x-1 \right)+P(2x-1)
\end{equation}
donde $P(z)$ es arbitraria impar polinomio, da un resultado proporcional a $\zeta(2p+1)$ cuando se integran como en $I_f$ definido anteriormente.