In an article by Ejsmont and Lehner a method using operator trace is applied to evaluate finite sums involving the cotangent. The paper gives numerous references concerning finite trigonometric sums which seem to have a rich recent history. We use here a result obtained by Berndt and Yeap by a contour integration technique. It reads
\begin{equation}
\sum_{n=1}^{q-1}\cot^{2k}\frac{n\pi}{q}=(-1)^kq-(-1)^k2^{2k}\sum_{\substack{j_0,j_1,\ldots,j_{2k}\ge0 \\j_0+j_1+\ldots+j_{2k}=k}}q^{2j_0}\prod_{p=0}^{2k}\frac{B_{2j_p}}{(2j_p)!}\tag{BY}\label{BY}
\end{equation}
where $B_s$ are the Bernoulli numbers. It seems that no similar results with odd powers of the function were found.
To evaluate
\begin{equation}
S_{2k}=\sum _{n=1}^m \left[\cot \left(\frac{n \,\pi }{2 m+1}\right)\right]^{2k}
\end{equation}
we remark that, by symmetry,
\begin{equation}
S_{2k}=\frac{1}{2}\sum _{n=1}^{2m} \left[\cot \left(\frac{n \,\pi }{2 m+1}\right)\right]^{2k}
\end{equation}
then, by (\ref{BY}),
\begin{equation}
(-1)^kS_{2k}=m+\frac{1}{2}-2^{2k-1}\sum_{\substack{j_0,j_1,\ldots,j_{2k}\ge0 \\j_0+j_1+\ldots+j_{2k}=k}}(2m+1)^{2j_0}\prod_{p=0}^{2k}\frac{B_{2j_p}}{(2j_p)!}
\end{equation}
Under this form, it appears that the sum is polynomial in $m$.
The summation is however difficult to perform. It can be transformed by the use of the series expansion of $\coth Z$
\begin{equation}
Z\coth Z=1+\sum_{s=1}^\infty\frac{B_{2s}}{(2s)!}(2Z)^{2s}
\end{equation}
One has then
\begin{align}
S_{2k}&=(-1)^k\left( m+\frac{1}{2}\right)\left[ 1-\left[Z^{-1}\right]\left( \coth\left( (2m+1)Z \right)\coth^{2k}(Z) \right) \right]\\
&=(-1)^k\left( m+\frac{1}{2}\right)\left[ 1-\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \coth\left( (2m+1)z \right)\coth^{2k}(z)\,dz \right]\label{eq:Cauchy}
\end{align}
where $\left[z^r\right]\left( F(z) \right)$ represents the coefficient of $z^r$ in the Laurent series of $F(z)$ at $z=0$.
The latter expression for $S_{2k}$ vanishes when $m=0$, as $\operatorname{Res}\left( \coth^{2k+1}z,z=0 \right)=1$. Moreover, from the identity $\coth 2z=(\coth z+(\coth z)^{-1})/2$, it also vanishes for $m=1/2$ as proposed. Then, polynomials $Q_k(m)$ exist such that
\begin{equation}
S_{2k}=m(2m-1)Q_k(m)
\end{equation}
where
\begin{align}
Q_k(m)&=\frac{(-1)^k(2m+1)}{2m(2m-1)}\left[ 1-\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \coth\left( (2m+1)z \right)\coth^{2k}(z)\,dz \right]\\
&=\frac{(-1)^k(2m+1)}{2m(2m-1)}\left[ 1-\left[Z^{-1}\right]\left( \coth\left( (2m+1)Z \right)\coth^{2k}(Z) \right) \right]
\end{align}
Evaluated with a CAS, the latter formula gives the results which are given in the OP up to $k=5$. Higher values of $k$ were numerically tested up to $S_{20}$ for several values of $m$.
Some values of $m$ are remarkable. For instance, it is not difficult to show that $Q_k(-1)=(-1)^{k+1}/3$ and $Q_k\left( -\frac{3}{2} \right)=\frac{(-1)^{k+1}}{3}$. Moreover $Q_k\left( \frac{3}{2} \right)=\frac{1}{3}$. Indeed, for $m=3/2$,
\begin{align}
\coth\left( (2m+1)z \right)&=\coth(4z)\\
&=\frac{1}{2}\left( \coth 2z+\frac{1}{\coth 2z} \right)\\
&=\frac{1}{4}\left( \coth z+\frac{1}{\coth z} \right)+\frac{\coth z}{1+\coth^2z}
\end{align}
the integral reads
\begin{align}
\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \coth\left( 4z \right)\coth^{2k}(z)\,dz &=\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \frac14\left( \coth^{2k+1}(z)+\coth^{2k-1}(z) \right)\,dz+I_k\\
&=\frac{1}{2}+I_k
\end{align}
where
\begin{equation}
I_k=\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \frac{\coth^{2k+1}(z)}{1+\coth^2z}\,dz
\end{equation}
We have the recurrence relation
\begin{equation}
I_k=1-I_{k-1}
\end{equation}
and
\begin{align}
I_0&=\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \frac{\coth(z)}{1+\coth^2z}\,dz\\
&=\frac{1}{2}\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon}\tanh(2z)\,dz\\
&=0
\end{align}
Thus $I_{2s}=0$ and $I_{2s+1}=1$. Finally,
\begin{equation}
Q_k\left( \frac32 \right)=
\begin{cases}
\frac23\left( 1-\frac{1}{2}-0\right)=\frac{1}{3}&\text{ when } k=2s\\
-\frac23\left( 1-\frac{1}{2}-1 \right)=\frac{1}{3}&\text{ when } k=2s+1
\end{cases}
\end{equation}
Thus $Q_k\left( \frac{3}{2} \right)=\frac{1}{3}$.
This also shows, by adapting the signs, that, for $m=-5/2$, the values of the polynomials oscillate between $1/3$ and $-1/5$. Adapting to the decomposition of $\coth 3z$, should also make it possible to show (as suggested by @Blue in a comment) that $Q_k(1)=3^{-k}$, which is numerically observed.
Best Answer
It's true in general. It suffices to show that $PQ$ has odd degree. Suppose to the contrary that there exist some $P,Q\in\mathbb F_2[x]$ subject to the constraints outlined in the OP for which $PQ=1+x+\cdots+x^{2n}=\frac{x^{2n+1}-1}{x-1}$. Observe by Frobenius that $\sum a_kx^{2k}=\left(\sum a_kx^k\right)^2$ in $\mathbb F_2[x]$. We shall prove that $1+x+\cdots+x^{2n}$ is squarefree; this clearly implies the result. Indeed, \begin{align*}\gcd\left(1+x+\cdots+x^{2n},\left(1+x+\cdots+x^{2n}\right)'\right)&=\gcd\left(1+x+\cdots+x^{2n},1+x^2+\cdots+x^{2n-2}\right)\\&=\gcd\left(1+x+\cdots+x^{2n},\left(1+x+\cdots+x^{n-1}\right)^2\right)\\&=\gcd\left(\tfrac{x^{2n+1}-1}{x-1},\left(\tfrac{x^{n}-1}{x-1}\right)^2\right)=1\end{align*} Where the last equality follows from the identity $\gcd(x^m-1,x^n-1)=x^{\gcd(m,n)}-1$ $\blacksquare$