Find $P(x,y,z)=x^n+y^n+z^n-\prod\limits_{k=0}^{n-1}(x+\omega_n^ky+\omega_n^{-k}z)$, where $\omega_n$ denotes a primitive $n$th root of unity 
Find $P(x,y,z)=x^n+y^n+z^n-\prod\limits_{k=0}^{n-1}(x+\omega_n^ky+\omega_n^{-k}z)$, where $\omega_n$ denotes a primitive $n$th root of unity.

I have manually multiplied the terms of the product and then equate the coefficients to get the polynomial but that's too cumbersome. 
Here is my method :
$(x + y + z)(x + y\omega_n + z\omega_n^{n-1})(x + y\omega_n^2 + z\omega_n^{n-2})....(x + y\omega_n^{n-1} + z\omega_n) = x^n(1 + [Y + Z])(1 + [Y\omega_n + Z\omega_n^{n-1}])(1 + [Y\omega_n^2 + Z\omega_n^{n-2}])....(1 + [Y\omega_n^{n-1} + Z\omega_n])$ where $Y=\frac {y}{x}$ and $Z=\frac {z}{x}$ 
Hence, I applied the formula:
$(1+\alpha)(1+\beta)(1+\gamma)...... = 1 + [\alpha + \beta + \gamma + ...] + [\alpha\beta + \beta\gamma + ....] + ....$
to get an expression for odd values of $n$ : 
$$P = nxyz(x^{n−3}+x^{n−5}yz+x^{n−7}y^2z^2+....)$$
My question is : How can I get a general expression for $P$ in a way better than what I have mentioned?
$$EDIT$$
I have a more general expression by now, which I think can be derived elementarily (unfortunately I still don't know how),
$$P=\frac {x^n}{t^n}(L_n(t)-t^n),$$
where $L_n(t)$ is the $n^{th}$ Lucas polynomial in $t:=\frac {ix}{\sqrt {yz}}$

We need to find a closed form expression of
$$
\prod_{k=0}^{n-1}(x+\omega^ky+\omega^{-k}z),
$$ where $\omega=e^{\frac{2\pi i}{n}}$ is the primitive $n$-th root of unity. Now let us consider the polynomial
$$
F(r) = r^n\prod_{k=0}^{n-1}(x+r\omega^ky+r^{-1}\omega^{-k}z)=\prod_{k=0}^{n-1}(xr+r^2\omega^ky+\omega^{-k}z)
$$
of degree $2n$ in $r$ regarding $x,y,z$ as fixed constants. By solving
$$
x+r\omega^ky+r^{-1}\omega^{-k}z = 0
$$ for each $0\le k\le n-1$,
we can see that $F(r)$ has simple roots
$$
r\omega^k = \frac{-x \pm \sqrt{x^2 - 4yz}}{2y}\Longrightarrow r = \omega^{-k}\frac{-x \pm \sqrt{x^2 - 4yz}}{2y}.
$$ 
If we denote $$\alpha(x,y,z) = \frac{-x + \sqrt{x^2 - 4yz}}{2y},\qquad  \beta(x,y,z)=\frac{-x - \sqrt{x^2 - 4yz}}{2y},$$ we see that the set of all roots of $F(r)$ coincides with that of $$G(r)=(r^n-\alpha^n(x,y,z))(r^n-\beta^n(x,y,z)).$$ Since neither $F$ nor $G$ has multiple roots, it follows that $F$ is a constant multiple of $G$, i.e. 
$$
F(r) = \left[\prod_{0\leq k\leq n-1} \omega^ky \right]G(r)=(-1)^{n-1}y^nG(r).
$$ by matching the leading coefficient. Plugging $r=1$ into the expression gives
\begin{align*}
\prod_{k=0}^{n-1}(x+\omega^ky+\omega^{-k}z)=&F(1)\\
=& (-1)^{n-1}y^nG(1)\\
 =&  (-1)^{n-1}y^n(1-\alpha^n-\beta^n +\alpha^n\beta^n).
\end{align*}Using $\alpha +\beta = -\tfrac x y$ and $\alpha\beta = \tfrac z y$, we get
$$
(-1)^{n-1}y^n(1-\alpha^n-\beta^n +\alpha^n\beta^n)=(-1)^{n-1}(y^n-(y\alpha)^n-(y\beta)^n+z^n ).
$$It remains to evaluate
$$
(y\alpha)^n+(y\beta)^n = \left(\frac{-x + \sqrt{x^2 - 4yz}}{2}\right)^n +  \left(\frac{-x - \sqrt{x^2 - 4yz}}{2}\right)^n.
$$ If we let $t= \frac{ix}{\sqrt{yz}}$, we have
\begin{align*}
(y\alpha)^n+(y\beta)^n =& (i\sqrt{yz})^n\left(\left(\frac{t + \sqrt{t^2 +4}}{2}\right)^n +  \left(\frac{t - \sqrt{t^2 +4}}{2}\right)^n\right) \\
=& (-1)^nx^nt^{-n}L_n(t)
\end{align*} where $L_n(t)$ is $n$-th Lucas polynomial. Finally, this gives
$$
\prod_{k=0}^{n-1}(x+\omega^ky+\omega^{-k}z) = (-1)^{n-1}y^n+ (-1)^{n-1}z^n+x^nt^{-n}L_n(t),
$$ hence
$$
P(x,y,z) = x^n(1-t^{-n}L_n(t)) + (1+(-1)^n)y^n + (1+(-1)^n)z^n.
$$
