I have, for a long period of time, tried to evaluate the indefinite integral $\int{\dfrac{x^m}{x^{2m}+1}dx}$ for different values of $m$, such as $m=2,3,4,6...$, and I have recently thought of generalizing this integral for all values of $m\in\mathbb{N}$. Since evaluating the integral using partial fractions for high values of $m$ by myself is very time consuming, and it happens that I can't find any direct solution on the internet, I tried to look for patterns of the integrals in WolframAlpha. After some comparisons, I found the below expression: $$\int{\dfrac{x^m}{x^{2m}+1}dx}=\sum_{n=1}^{m}(-1)^{n-1}(\dfrac{1}{2m}\sin(\dfrac{(2n-1)\pi}{2m})\ln{|x^2+2\cos(\dfrac{(2n-1)\pi}{2m})x+1|}$$$$+\dfrac{1}{m}\arctan(\csc(\dfrac{(2n-1)\pi}{2m})(x+\cos(\dfrac{(2n-1)\pi}{2m})))\cos(\dfrac{(2n-1)\pi}{2m}))$$ How could I have found the above expression without comparing the integrals?
P.S. It seems that the equivalence $x^{2m}-2x^ma^m\cos(m\theta)+a^{2m}\equiv\displaystyle\prod_{r=0}^{m-1}(x^2-2ax\cos(\theta+\dfrac{2r\pi}{m})+a^2)$ for $a\in\mathbb{R}^+, \theta\in\mathbb{R}, m\in\mathbb{N}$ can be used, yet I don't know how I can evaluate partial fractions with it.