I came across a formula for the partial fraction decomposition of $ \displaystyle \frac{1}{x^{2n}+a^{2n}}$.
It seems correct (at least for $n=1,2$, and $3$).
But how is it derived?
$$\frac{1}{x^{2n}+a^{2n}} = \frac{1}{na^{2n-1}} \sum_{k=0}^{n-1} \frac{a -x \cos (\frac{2k+1}{2n} \pi)}{\Big(x-a \cos (\frac{2k+1}{2n} \pi) \Big)^{2} +a^{2} \sin^{2} (\frac{2k+1}{2n} \pi) }$$
EDIT:
The denominators of the partial fractions can be rewritten as $$x^{2}-2ax\cos \left(\frac{2k+1}{2n} \pi \right) +a^{2}$$
If you factor $ \displaystyle x^{2n}+a^{2n}$ over the complex numbers and then rearrange terms and multiply pairwise you get quadratic factors like that.
So at least that makes sense.
I) First consider the well-known partial fraction decomposition formula
$$ \tag{1}\frac{1}{P(z)}~=~ \sum_{k=1}^N\frac{1}{P^{\prime}(a_k)(z-a_k)}$$
which is valid whenever
$$ \tag{2} P(z)~=~A\prod_{k=1}^N (z-a_k)$$
is an $N$th-order polynomial with $N$ distinct roots $a_1, \ldots, a_N$.
II) Next consider OP's special case $N=2n$ and
$$ \tag{3} P(z)~=~z^{2n}+a^{2n},$$
where $a>0$. The roots of the polynomial (3) are
$$ \tag{4} a_k~=~a\omega_k, \qquad \omega_k~:=~\exp\left[\frac{i\pi(2k-1)}{2n}\right], \qquad k=1, \ldots, 2n,$$
or equivalently $a_1, \ldots, a_n$, and their complex conjugate. Therefore
$$\frac{1}{P(z)}~=~ \sum_{k=1}^{2n}\frac{1}{P^{\prime}(a_k)(z-a_k)} ~=~\frac{1}{2n}\sum_{k=1}^{2n}\frac{a_k}{a_k^{2n}(z-a_k)} $$ $$ \tag{5}~=~\frac{-1}{2na^{2n-1}}\sum_{k=1}^{n}\left(\frac{\omega_k}{z-a_k} +\frac{\bar{\omega}_k}{z-\bar{a}_k} \right) ~=~\frac{1}{na^{2n-1}}\sum_{k=1}^{n}\frac{a-\frac{\omega_k+\bar{\omega}_k}{2}z}{(z-a_k)(z-\bar{a}_k)}. $$
The formula (5) readily leads to OP's partial fraction decomposition formula.