We have the following key lemma. Lemma: Consider $t\in\mathbb{R}$ such that the roots of $x^2-tx+1$ are complex, but not roots of unity. Then, $P(x,y,t)$ is a polynomial in $x^2+y^2-xyt+t^2$. Proof: Let $Q(x,y)=P(x,y,t)$. It suffices to show that $Q$ is a polynomial in $x^2+y^2-xyt$. We have \[Q(x,y) = Q(x,tx-y) = Q(ty-x).\quad\quad\quad(\star)\]Since $(\star)$ preserves degree, it suffices to prove if $Q$ is homogeneous of degree $d$. Now, $(\star)$ is a polynomial identity, so it holds if we evaluate it at a complex pair. Let $\omega$ and $\bar{\omega}$ be the complex roots of $x^2-tx+1$. Let $R(y)=Q(1,y)$. We have \[R(\omega)=Q(1,\omega)=Q(1,t-\omega)=Q(1,\bar{\omega})=R(\bar{\omega}),\]as $\omega+\bar{\omega}=t$. But $\overline{R(\omega)}=R(\bar{\omega})$, so $k:=R(\omega)=R(\bar{\omega})\in\mathbb{R}$, so both $\omega$ and $\bar{\omega}$ are roots of $R(y)-k$. Thus, \[y^2-ty+1=(y-\omega)(y-\bar{\omega})\mid R(y)-k,\]so \[R(y) = k+(y^2-ty+1)\cdot A(y),\]so by homogeneity, \[Q(x,y) = kx^d + (x^2-txy+y^2)B(x,y),\]where $B$ is homogeneous of degree $d-2$. But \[Q(\omega,1)=Q(t-\omega,1)=Q(\bar{\omega},1),\]so \[k\omega^d = k\bar{\omega}^d,\]so $\omega^{2d}=1$ if $k\ne 0$, which is impossible as $\omega$ is not a root of unity. Thus, we have $k=0$, so \[Q(x,y) = (x^2-txy+y^2)B(x,y),\]so $B$ satisfies $(\star)$, and we induct down to prove $Q(x,y)=(x^2-txy+y^2)^{d/2}$. This completes the proof of the lemma. $\blacksquare$ Note that all transcendental $t$ with absolute value less than $2$ satisfy the conditions of the lemma. For there values of $t$, we have functions $\alpha_k(t)$ such that \[P(x,y,t) = \sum_{k=0}^d\alpha_k(t)(x^2+y^2-txy+t^2)^k\quad\quad\quad(\star\star).\]By matching $x^ay^b$ coefficients on the sides of the above equation, we see that it is equivalent to the vector $(\alpha_0,\ldots,\alpha_d)$ satisfying a set of linear equations with coefficients polynomials in $t$. Due to the lemma, we know the existence of at least one solution, but since all the coefficients in the equations are polynomials in $t$, by linear algebra, we may find at least one solution for $(\alpha_0,\ldots,\alpha_d)$ that are rational functions in $t$. Thus, WLOG, we may assume that the $\alpha_k(t)$s are rational functions in $t$. In particular, this implies that $(\star\star)$ actually holds for all but finitely many values of $t$. Now, $(\star\star)$ implies the existence of $A(z)$ and $B(x,y,z)$, the former a rational function in $z$, and the latter a polynomial in $x,y,z$ divided by a polynomial in $z$ such that \[P(x,y,z) = A(z) + (x^2+y^2+z^2-xyz)B(x,y,z).\]Now, by similar logic, there exists $C(y)$ and $D(x,y,z)$, the former a rational function in $y$, and the latter a polynomial in $x,y,z$ divided by a polynomial in $y$ such that \[P(x,y,z) = C(y) + (x^2+y^2+z^2-xyz)D(x,y,z).\]This implies that \[A(z)-C(y) = (x^2+y^2+z^2-xyz)(D(x,y,z)-B(x,y,z)).\]Thus, $x^2+y^2+z^2-xyz$ divides the numerator of $A(z)-C(y)$ which is a polynomial in $y$ and $z$ (in particular, $x$-degree equal to $0$), which is not possible unless it is $0$, so $A(z)=C(y)=k$ for some constant $k$, so \[x^2+y^2+z^2-xyz\mid P-k.\]It is easy to see that $\frac{P-k}{x^2+y^2+z^2-xyz}$ also satisfies the conditions of the problem, so we're done.

We present the first solution from the official shortlist. Let's say a polynomial is good if $P(x,y,z) = P(x,y,xy-z)$ and is very good if it satisfies the problem condition. We will mostly only work with good polynomials. Claim: A polynomial is good if and only if it is the form $Q(x,y,z(xy-z))$ for some three-variable polynomial $Q$. Proof. One direction is clear. For the other direction, if $P$ is good, let $w = z + \frac{1}{2} xy$, i.e.\ consider the shifted polynomial \[ P_1(x,y,w) \overset{\text{def}}{=} P\left( x,y, z + \tfrac12 xy \right). \]Then $P$ is good iff $P_1$ is even in $w$, so $P_1$ is a polynomial in $x$, $y$, $w^2$. As $w^2 = -z(xy-z) + \frac14x^2y^2$, the result follows. $\blacksquare$ Claim: Suppose $P = Q(x,y,z(xy-z))$ is good. Consider all the monomials in $P$ of maximum degree. For every such monomial, the degree of $z$ is at most the degrees of $x$ and $y$. Proof. Expand $Q$ to get an expression of form \[ P = \sum_{a,b,c} \mu_{a,b,c} x^a y^b \left( z(xy-z) \right)^c. \]In each individual summand, the largest degree term is $\mu_{a,b,c} x^{a+c} y^{b+c} z^{c}$, of degree $a+b+3c$. If $m = \deg P$, then no two terms with $a+b+3c=m$ may cancel. This implies the conclusion. $\blacksquare$ Finally, assume $P$ is very good. Actually its maximum degree term must be of the form $\mu(xyz)^n$, by applying the claim three times. But then it follows \[ P(x,y,z) - \mu(x^2+y^2+z^2-xyz)^n \]satisfies the problem condition and also has lesser degree. The problem is solved by induction on $\deg P$. Remark: [Michael Ren] The problem is true in a bit of generality. Suppose $P(x_1, \dots, x_n)$ is a multivariable polynomial satisfying \[ P(x_1, \dots, x_n) = P\left(x_1, \dots, x_{n-1}, x_1x_2 \dots x_{n-1} - x_n\right) \]then one might hope that $P$ is a polynomial in $x_1 x_2 \dots x_n - x_1^2 - x_2^2 \dots - x_n^2$. This is true for all $n > 2$ and for $n = 1$, but it is not true for $n = 2$; a counterexample is given by $P(x,y) = x^6 + y^6 + (x-y)^6$.

We will use the sophisticated method of "Haha cosine identities go brrrrrrrrrrrrrrrr" We will heavily rely on the following identity :$$a+b+c=0 \implies \cos^2 a + \cos^2 b + \cos^2 c = 2 \cos a \cos b \cos c +1$$ The idea is to consider the function $Q(a,b,c)$ such that $Q(a,b,c)= P(2\cos a , 2\cos b, 2\cos c)$. We will work on the plane $a+b+c=0$. Note that the shifts $(x,y,z) \mapsto (yz-x,y,z)$ corresponds to the shift $(a,b,c) \mapsto (b-c,-b,c)$ and similarly for other cyclic variants. Also note that $Q(a,b,c)=Q(-a,-b,-c)$ Claim : $Q$ is constant over the plane $a+b+c=0$. Proof : We have $$Q(a,b,c)=Q(b-c,-b,c)=Q(c-b,-b,2b-c)=Q(a+2b,b,c-2b)$$ Hence by induction, we have for all integers $k$ the following equality : $Q(a,b,c)=Q(a+2bk,b,c-2bk)$ and other cyclic variants. Since $Q(a,b,c)=Q(a+2\pi b,b,c-2\pi b)$, we can choose $(a,b,\pi)$ to be linearly independent over $\mathbb Q$ and hence we get that $Q(x,y,z)$ is constant over $x+y+z=0$. Let this constant value be $k$. Claim : $x^2+y^2+z^2-xyz-4 \vert P(x,y,z)-k$ Proof : Suppose we have $$P(x,y,z)-k=(x^2+y^2+z^2-xyz-4)\cdot R(x,y,z) + A(x,y)\cdot z + B(x,y)$$ By considering $(x,y,z)=(2\cos a , 2\cos b, 2\cos c)$ with $a+b+c=0$ we have that $A(x,y)\cdot z + B(x,y)$ vanishes on this region. Now we will construct two different values of $z$ for a pair of $(x,y)$ such that $A(x,y)\cdot z + B(x,y)=0$. This will prove that $A(x,y)=B(x,y)=0$. Note that we can choose $a,b \in \left( \frac {\pi}3, \frac {2\pi}3\right )$. Then we have two values of $z$ namely $z=2\cos(a+b)<0$ and $z=2\cos(a-b)>0$. This means we must have $A(x,y)=B(x,y)=0$ as desired. Finally note that the polynomial $\frac {P(x,y,z)-k}{x^2+y^2+z^2-xyz-4}$ is of smaller degree than $P(x,y,z)$ and also satisfies the equation. Since the polynomial $T(x,y,z)=x^2+y^2+z^2-xyz$ itself satisfies the property, we are done by induction on the degree of the polynomial. $\blacksquare$

It is problem proposed by Russia, author - G. Chelnokov maybe together with I. Bogdanov

Claim: A polynomial is good if and only if it is the form $Q(x,y,z(xy-z))$ for some three-variable polynomial $Q$. Proof. One direction is clear. For the other direction, if $P$ is good, let $w = z + \frac{1}{2} xy$, i.e.\ consider the shifted polynomial\[ P_1(x,y,w) \overset{\text{def}}{=} P\left( x,y, z + \tfrac12 xy \right). \]Then $P$ is good iff $P_1$ is even in $w$, so $P_1$ is a polynomial in $x$, $y$, $w^2$. As $w^2 = -z(xy-z) + \frac14x^2y^2$, the result follows. $\blacksquare$ Claim: Suppose $P = Q(x,y,z(xy-z))$ is good. Consider all the monomials in $P$ of maximum degree. For every such monomial, the degree of $z$ is at most the degrees of $x$ and $y$. Proof. Expand $Q$ to get an expression of form\[ P = \sum_{a,b,c} \mu_{a,b,c} x^a y^b \left( z(xy-z) \right)^c. \]In each individual summand, the largest degree term is $\mu_{a,b,c} x^{a+c} y^{b+c} z^{c}$, of degree $a+b+3c$. If $m = \deg P$, then no two terms with $a+b+3c=m$ may cancel. This implies the conclusion. $\blacksquare$ Finally, assume $P$ is very good. Actually its maximum degree term must be of the form $\mu(xyz)^n$, by applying the claim three times. But then it follows\[ P(x,y,z) - \mu(x^2+y^2+z^2-xyz)^n \]satisfies the problem condition and also has lesser degree. The problem is solved by induction on $\deg P$. Remark: [Michael Ren] The problem is true in a bit of generality. Suppose $P(x_1, \dots, x_n)$ is a multivariable polynomial satisfying\[ P(x_1, \dots, x_n) = P\left(x_1, \dots, x_{n-1}, x_1x_2 \dots x_{n-1} - x_n\right) \]then one might hope that $P$ is a polynomial in $x_1 x_2 \dots x_n - x_1^2 - x_2^2 \dots - x_n^2$. This is true for all $n > 2$ and for $n = 1$, but it is not true for $n = 2$; a counterexample is given by $P(x,y) = x^6 + y^6 + (x-y)^6$. There are lots of sols though this official one is nice Credits to v_Enhance for making the solution as I am using his solution, sorry.

Aryan-23 wrote: We will use the sophisticated method of "Haha cosine identities go brrrrrrrrrrrrrrrr" We will heavily rely on the following identity :$$a+b+c=0 \implies \cos^2 a + \cos^2 b + \cos^2 c = 2 \cos a \cos b \cos c +1$$ The idea is to consider the polynomial $Q(a,b,c)$ such that $Q(a,b,c)= P(2\cos a , 2\cos b, 2\cos c)$. We will work on the plane $a+b+c=0$. Note that the shifts $(x,y,z) \mapsto (yz-x,y,z)$ corresponds to the shift $(a,b,c) \mapsto (b-c,-b,c)$ and similarly for other cyclic variants. Also note that $Q(a,b,c)=Q(-a,-b,-c)$ Claim : $Q$ is constant over the plane $a+b+c=0$. Proof : We have $$Q(a,b,c)=Q(b-c,-b,c)=Q(c-b,-b,2b-c)=Q(a+2b,b,c-2b)$$ Hence by induction, we have for all integers $k$ the following equality : $Q(a,b,c)=Q(a+2bk,b,c-2bk)$ and other cyclic variants. Since $Q(a,b,c)=Q(a+2\pi b,b,c-2\pi b)$, we can choose $(a,b,\pi)$ to be linearly independent over $\mathbb Q$ and hence we get that $Q(x,y,z)$ is constant over $x+y+z=0$. Let this constant value be $k$. Claim : $x^2+y^2+z^2-xyz-4 \vert P(x,y,z)-k$ Proof : Suppose we have $$P(x,y,z)-k=(x^2+y^2+z^2-xyz-4)\cdot R(x,y,z) + A(x,y)\cdot z + B(x,y)$$ By considering $(x,y,z)=(2\cos a , 2\cos b, 2\cos c)$ with $a+b+c=0$ we have that $A(x,y)\cdot z + B(x,y)$ vanishes on this region. Now we will construct two different values of $z$ for a pair of $(x,y)$ such that $A(x,y)\cdot z + B(x,y)=0$. This will prove that $A(x,y)=B(x,y)=0$. Note that we can choose $a,b \in \left( \frac {\pi}3, \frac {2\pi}3\right )$. Then we have two values of $z$ namely $z=2\cos(a+b)<0$ and $z=2\cos(a-b)>0$. This means we must have $A(x,y)=B(x,y)=0$ as desired. Finally note that the polynomial $\frac {P(x,y,z)-k}{x^2+y^2+z^2-xyz-4}$ is of smaller degree than $P(x,y,z)$ and also satisfies the equation. Since the polynomial $T(x,y,z)=x^2+y^2+z^2-xyz$ itself satisfies the property, we are done by induction on the degree of the polynomial. $\blacksquare$ sorry , but i was thinking on this problem , and i had the same aproach as you , but i found it incorrect , as you can't say $Q$ is a polynomial , for an example if $P(x,y,z)=x+y+z$ then $Q(x,y,z)=sinx +sinz +siny$ this clearly is not a polynomial. so , can you alaborate on how your proof is complete? bc im really confused here ...

okay yeah thanks, I should have written consider the function $Q$ defined as blah blah. I dont think that affects the proof at any stage as $Q$ was really defined for notational convenience. What we are really doing is analyzing the behavior of $P$ over the region : $$\sigma = \{(2\cos a , 2\cos b , 2 \cos c )\vert a+b+c=0 \} $$ @below I'm not claiming that $x^2+y^2+z^2-xyz-4 \vert P(x,y,z)-Q(x,y,z)$ for all $(x,y,z)$. I'm claiming that $x^2+y^2+z^2-xyz-4 \vert P(x,y,z)-k$ (for all x,y,z) , where $k$ is the constant value attained by $Q(a,b,c)$ over $a+b+c=0$ . So like $k$ is fixed say $k= Q(-1, \tfrac 12 , \tfrac 12)$

well you also used $x^2+y^2+z^2-xyz-4|p(x,y,z)-Q(x,y,z)$ but its not a polynomial ...or maybe i'm again misstaking ... the proof for $Q(x,y,z)=C$ , how do you know ,its true if its not a polynomial? @above ... yeah right , i get it ... thanks

Is this too simple ... ? We first attempt to focus on one singular expression $P(x,y,z) = P(x,y,xy-z)$, to first eliminate polynomials which $\textit{couldn't even maintain one of the three equalities}$. $\textbf{Claim 1.}$ $P(x,y,z) = P(x,y,xy-z)$ implies $P$ is a polynomial in $z(xy-z),x$ and $y$. $\textbf{Proof 1.}$ We focus on the $z$ terms, and so we can let \[ P(z) = \sum z^i P_i(x,y) \]The given statement turns into $P(z) = P(xy-z)$, and further setting $z = z+\dfrac{1}{2}xy$ gets $P\left(z+\dfrac{1}{2}xy\right) = P\left(-z+\dfrac{1}{2}xy\right)$. To "translate" $P$, set $Q(z) = P\left(z + \dfrac{1}{2}xy\right)$, and express $Q$ similarly as \[Q(z) = \sum z^i Q_i(x,y) \]As $Q(z) = Q(-z)$, we can infer that $Q_i(x,y)$ is $0$ when $i$ is odd. Thus we can ignore the odd terms and say that $Q$ is a polynomial in $x,y$ and $z^2$. This brings us full circle to knowing the complete description in $P$, that $P$ satisfies the (first) equality if and only if it is a polynomial in $x,y$ and $\left(z-\dfrac{1}{2}xy\right)^2$. Since $\left(z-\dfrac{1}{2}xy\right)^2 = z^2-xyz + (\text{something in}$ $x$ $\text{and}$ $y)$, expanding again yields the desired conclusion. After applying it again for $P(x,y,z) = P(x,zx-y,z)$ and $P(x,y,z) = P(yz-x,y,z)$, we get the assertion that \[ P(x,y,z) = F_1(x,y,z(xy-z)) = F_2(y,z,x(yz-x)) = F_3 (z,x,y(zx-y)) \]Seeing that $z(xy-z) = -(x^2+y^2+z^2-xyz) + x^2+y^2$, we can use expansion $\textbf{again}$ to get that \begin{align*} F_1(x,y,z(xy-z)) &= \sum (xyz-z^2)^i X_i(x,y) = \sum (-(x^2+y^2+z^2-xyz) + x^2+y^2)^i X_i(x,y) \\ &= \sum \sum_{j \leq i} \binom{i}{j} (-(x^2+y^2+z^2-xyz))^j (x^2+y^2)^{i-j} X_i(x,y) \end{align*}to get further that \begin{align*} P(x,y,z) &\equiv G_1(x,y) \ \text{mod} \ (x^2+y^2+z^2-xyz) \\ &\equiv G_2(y,z) \ \text{mod} \ (x^2+y^2+z^2-xyz) \\ &\equiv G_3(z,x) \ \text{mod} \ (x^2+y^2+z^2-xyz) \end{align*}We then finish with the following claim: $\textbf{Claim 2.}$ $x^2+y^2+z^2-xyz \mid G_1(x,y)-G_2(y,z)$ implies $G_1(x,y) = G_2(y,z) = H_{1,2}(y)$. $\textbf{Proof 2.}$ Let the degree of $H = \dfrac{G_1(x,y)-G_2(y,z)}{x^2+y^2+z^2-xyz}$ be $d$, and now we consider the term with the highest $x-power$, with the highest$y-power$ among them. Say it is $x^iy^jz^{d-i-j}$. Simple observation yields that the term $x^{i+1}y^{j+1}z^{d-i-j+1}$ in $G_1(x,y)-G_2(y,z)$ is nonzero and uncancelled (terms of degree $d-1$ and below in $H$ will yield lesser powers, and all other terms of degree $d$ in $H$ can only match with $xyz$ to create terms of degree $d+3$ in $G_1(x,y)-G_2(y,z)$.) This directly provides us a contradiction, so $H$ must be zero. Then $G_1(x,y) = G_2(y,z)$, and this can only happen if $G_i$ only consists of terms with $y$; i.e. both of them are equal to some polynomial $H_{1,2}$ in $y$. The claim's results then yield \[ G_1(x,y) = G_2(y,z) = G_3(z,x) = H_{1,2}(y) = H_{2,3}(z) = H_{3,1}(x) \]and that gives us the final conlusion that all of them are constant polynomials. Ultimately, that means that \[ x^2+y^2+z^2-xyz | P(x,y,z)-c \]From this, we can let $P(x,y,z) = P_0(x,y,z) = (x^2+y^2+z^2-xyz)\cdot P_{-1}(x,y,z) + c$. Because $x^2+y^2+z^2-xyz = x^2+y^2+(xy-z)^2-xy(xy-z)$, we can "induct downwards", as \[ P_{-1}(x, y, z)=P_{-1}(x, y, xy-z)=P_{-1}(x, zx-y, z)=P_{-1}(yz-x, y, z) \]This cannot happen indefinitely, so there must exist a $d$ such that $P_{-d} = \text{constant}$. Retracing back finalizes the construction of the desired polynomial $F$.

Is my solution even correct ? Lemma - Let $I$ be an infinite integral domain and $P(x) \in I[x]$, then 1) If $P(c) = 0$ for $c\in I$ then $(x-c)$ divides $P(x)$. 2) If $P(t) = 0$ for infinitely many $t\in I$ then $P(x) \equiv 0$. 3) If there exists $r\in I$ such that $P(x) = P(r-x)$, then there exists $Q\in I[x]$ such that $P(x)=Q(x(x-r))$. 4) In 3), there actually exists a $Q$ such that $P(x) = Q(x(x-r)+s)$ for any fixed $s \in I$.

Now we use these to finish the problem. Let $t = x^2+y^2+z^2-xyz$. Taking $I = \mathbb{R}[y][z],r = yz , s = y^2+z^2$, we have $$P(x,y,z) = g(t) = g_nt^n+\cdots g_1t+g_0 \hspace{5mm} \text{where } g_i \in I$$Now we will show that all $g_i$ are in fact real numbers. Assume to the contrary . Let $k$ be the largest such that $g_k \notin \mathbb{R}$, assume that $g_k$ has a term with $y$ Now we replace $y$ by $xz - y$, note that $t$ remains unchanged and so do $g_u$ for $u>k$. Let $g_l$ transform to $h_l(x)$ for $l\leq k$, as $g_k$ had a term with $y$ we get that $h_k(x)$ is not constant poly in $x$. But we also have $$g_nt^n+\cdots g_1t+g_0 = g_nt^n+\cdots+g_{k+1}t^{k+1}+h_k(x)t^k+\cdots+h_0(x)$$$$ \implies g_kt^k+\cdots+g_0=h_k(x)t^k+\cdots+h_0(x)$$But here LHS has degree $2k$ whereas RHS has degree $(2k + $ deg$(h_k(x))) \neq 2k $ . And this is a contradiction . Hence we must have that all $g_k \in \mathbb{R}$ and we are done by taking $F = g$

First focus only on $P(x,y,z) = P(x,y,xy-z)$. I claim that this can be written as a polynomial in $xyz - z^2$ with coefficients as polynomials in $x,y$. Let $Q$ be the portion of $P$ divisible by $z$, so we have $Q(x,y,z) = Q(x,y,xy-z)$ too. Note that $Q(x,y,xy) = Q(x,y,0) = 0$. For each monomial in $Q$, say $x^ay^bz^c$, define its value to be the ordered pair $(a+c, b+c)$. Clearly, only stuff with the same value can cancel out when we put $z = xy$ ($Q$ is $0$ as a polynomial). Among things with the same value, pair up monomials (possibly by splitting them) so that each cancels the other. So we have each pair as something of the form (assume $c > m$) $c_1(x^ay^bz^c - x^{a+c-m}y^{b+c-m}z^m) = x^ay^bz^m(z^{c-m} - (xy)^{c-m})$, which has $z-xy$ as a factor. Now ignoring the $z(z-xy)$, we repeat the same thing on the remaining things so we indeed have that this is a polynomial in $x,y,xyz-z^2$. Now, the clever bit. Out of all monomials in $P$, consider the one divisible by $xyz$ and with maximal degree, call it the important term. Since we have that its a polynomial in $xyz - z^2$, the important term will be something of the form $cx^ay^b(xyz)^n$, in which clearly $z$ has minimal degree. But we can repeat this for other variables to get that all are in fact equal, so the important term is of the form $c(xyz)^n$. Now just subtract a $c(x^2 + y^2 + z^2 - xyz)^n$ from this, which clearly reduces degree, and now just induct on the degree, so we are done. $\blacksquare$

Solved with rama1728. Make the substitution $w = z - \frac{1}{2}xy.$ Then $$P(x,y,w+\frac{1}{2}xy) = P(x,y,-w+\frac{1}{2}xy).$$Treat $P$ as a polynomial in $w$ with coefficients in terms of $x$ and $y.$ Then we know it is even, so it can also be treated as a polynomial in $w^2 = z^2 - xyz + \frac{1}{4}x^2y^2.$ Thus, $P$ can be written as a polynomial in $z^2-xyz$ with coefficients in terms of $x$ and $y.$ So \[ P(x,y,z) = \sum_{i,j,k} (p_{i,j,k}) x^i y^j \left( z^2-xyz\right)^k\]where the sum is taken over nonnegative integers $i,j,k$ and $p_{i,j,k}$ is a constant depending on $i,j,k.$ The term of maximal degree in each summand is $(p_{i,j,k})x^{i+k}y^{j+k}z^k.$ So all terms of maximal degree in $P(x,y,z)$ must be of the form $(p_{i,j,k})x^{i+k}y^{j+k}z^k$ for nonnegative integers $i,j,k.$ Applying a symmetrical argument shows that all the terms of maximal degree in $P(x,y,z)$ must be of the form $C(xyz)^n$ for some nonnegative integer $n = \deg P$ and a constant $C.$ So for some constant $c,$ \[ P(x,y,z) - c(x^2+y^2+z^2-xyz)^n \]satisfies the conditions and has less degree, so we just induct down on the degree of $P$ until $F(t)$ is generated. $\blacksquare$

Let's perform a process on $P$. Take $f(y,z)x^{2k}$ to be the sum of the terms with maximal even degree in $x$. We may write \[f(y,z)x^{2k} = f(y,z)(x^2+y^2+z^2-xyz)^k - f'(x,y,z)\]where $f'$ has a maximal even degree in $x$ less than $2k$. Since $P(x,y,z)-f(y,z)x^{2k}$ also has maximal even degree in $x$ less than $2k$, $P(x,y,z)-f(y,z)(x^2+y^2+z^2-xyz)^k$ has maximal even degree in $x$ less than $2k$. We may repeat this process, eliminating the highest even degree in $x$ each time. After a finite number of steps, we obtain a polynomial with no even degrees in $x$. Effectively, this polynomial is odd in $x$. \[P(x,y,z) - \sum_{i=0}^k f_i(y,z)(x^2+y^2+z^2-xyz)^i = P_0(x,y,z),\]where $P_0(x,y,z)=-P_0(-x,y,z)$. We may now write \[P(x,y,z)=Q(x^2+y^2+z^2-xyz,y,z) + P_0(x,y,z)\]for some polynomial $Q$. For the moment, let $Q'(x,y,z)=Q(x^2+y^2+z^2-xyz,y,z)$. Note that \[Q'(x,y,z)=Q(y^2+z^2-x(yz-x))=Q'(yz-x,y,z).\]Since $P(x,y,z)=P(yz-x,y,z)$, we must have $P_0(x,y,z)=P_0(yz-x,y,z)$ as well. Hence, \[P_0(x,y,z)=-P_0(-x,y,z)=-P_0(yz+x,y,z)\]\[\implies P_0(x,y,z)=P_0(2yz+x,y,z).\]However, this means $P_0$ is periodic about $x$. Hence $P_0$ evaluates to $0$ everywhere, as it is odd. We write \[P(x,y,z) = Q(x^2+y^2+z^2-xyz,y,z),\]as desired. $\square$ We've yet to use the other two equalities. We have \[Q(x^2+y^2+z^2-xyz,y,z)=Q(x^2+y^2+z^2-xyz,y,xy-z).\]There are four quantities $a=x^2+y^2+z^2-xyz$, $b=y$, $c=z$, and $d=xy-z$ in this equation. We explore the relationship between these four quantities. In particular, if $b \neq 0$, $x=\tfrac{c+d}{b}$ and thus \[a=\left(\frac{c+d}{b}\right)^2+b^2-cd.\]By including $b=0$, we may rearrange this as $c^2-(b^2-2)dc+d^2=ab^2-b^4$. As long as this condition holds, $Q(a,b,c)=Q(a,b,d)$. We will define a sequence $\{c_n\}$ where $Q$ evaluates to a constant. Suppose we fix $c_1$ and $c_2$. Borrowing from Vieta jumping, we define the remainder of the sequence by $c_{i+1}=(b^2-2)c_i-c_{i-1}$. We find that $Q(a,b,c_i)=Q(a,b,c_{i+1})$, realizing our goal. Now, we can note that only finitely many pairs of $c_1$ and $c_2$ give periodic sequences. By picking $c_1$ and $c_2$ that gives an infinite sequence, we find that $Q$ must not depend on $z$. An analogous argument determines that $Q$ must not depend on $y$. Hence, we may write $P(x)=Q(x^2+y^2+z^2-xyz,y,z)=F(x^2+y^2+z^2-xyz)$ for some polynomial $F$.

Let $x^ay^bz^c$ be a monomial with greatest degree in $P$ (hence has non-zero coefficient). Claim 1: $a = b = c$. Proof: Suppose FTSOC that $a \le b \le c$ with $a < c$ (other variations follow by symmetrical arguments). Let $f(\alpha, \beta, \gamma)$ denote $[x^\alpha y^\beta z^\gamma]P(x, y, z)$, with the abbreviation $f(t) = f(a-t, b-t, c+t)$. Note that if $a - t < 0$, $f(t) = 0$. Furthermore, by maximality of $a+b+c$, for $t < 0$ we also have $f(t) = 0$, otherwise $\deg(x^{a-t}y^{b-t}z^{c+t}) = a+b+c-t > a+b+c$. \[x^\alpha y^\beta (xy-z)^\gamma = x^\alpha y^\beta \sum_{r = 0}^\gamma {\gamma \choose r} (-1)^{\gamma-r} (xy)^r z^{\gamma - r} = \sum_{r = 0}^\gamma {\gamma \choose r} (-1)^{\gamma-r} x^{\alpha+r} y^{\beta+r} z^{\gamma - r}.\]Hence, if we have monomial $A(x, y, z)$ and any monomial $B(x, y, z)$ in $A(x, y, xy-z)$, then $\deg_y(A) - \deg_x(A) = \deg_y(B) - \deg_x(B)$. Furthermore, $\deg_x(A) + \deg_z(A) = \deg_x(B) + \deg_z(B)$; these properties are invariant. We focus on $P(x, y, z) = P(x, y, xy-z)$. Consider $Q(x, y, z)$, the subset of monomials in $P(x, y, z)$ with $\deg_y - \deg_x = b-a$ and $\deg_x + \deg_z = a+c$. Thus we have $Q(x, y, z) = Q(x, y, xy-z)$. However, by the invariant properties above, a monomial lies in $Q(x, y, xy-z)$ iff its 'pre-image' before $z \to xy-z$ lies in $Q(x, y, z)$. Hence, we get: \begin{align*} \sum_{t = 0}^a f(t) x^{a-t}y^{b-t}z^{c+t} &= Q(x, y, z) = Q(x, y, xy-z) \\ &= \sum_{t = 0}^a f(t) \left(\sum_{r = 0}^{c+t} {c+t \choose r} (-1)^{c+t-r} x^{a-t+r}y^{b-t+r}z^{c+t-r}\right) \\ &= \sum_{p = -c}^a x^{a-p}y^{b-p}z^{c+p} \left(\sum_{t = p}^a f(t) {c+t \choose t-p} (-1)^{c+p} \right). \end{align*}Now by equating coefficients, for $p = -c, -c+1, ..., -1$, we necessarily have: \[\sum_{t=p}^a f(t) {c+t \choose t-p} (-1)^{c+p} = 0 \Rightarrow \sum_{t = 0}^a f(t) {t+c \choose t-p} = 0.\]We have $a+1$ variables $f(0), f(1), ..., f(a)$, and $c \ge a+1$ equations. We focus on $p = -c, ..., -c+a$, and define $S_p = \sum_{t=0}^a f(t) {t+c \choose t-p} = \sum_{t=0}^a f(t) {t+c \choose c+p}$. Suppose $S_{-c}, S_{-c+1}, ..., S_{-c+a}$ aren't linearly independent, i.e. there exist $\mu_{-c}, ..., \mu_{-c+a}$ not all zero such that: \[\sum_{p = -c}^{-c+a} \mu_pS_p = 0.\]This occurs iff for every $t = 0, 1, ..., a$, we have: \[\sum_{p = -c}^{-c+a} \mu_p {t+c \choose c+p} = 0 \iff \sum_{q = 0}^a \mu_{q-c} {t+c \choose q} = 0.\]But notice that ${t+c \choose q} = \frac{(t+c)(t+c-1)...(t+c-q+1)}{q!}$ is a degree $q$ polynomial in $t$, hence if $k-c$ is the largest in $\{-c, ..., -c+a\}$ whereby $\mu_{k-c} \neq 0$, then above we have a degree $k$ polynomial in $t$. If $k = 0$, then $\mu_{-c} {t+c \choose 0} = \mu_{-c} = 0$ contradiction. Hence, $0 < k \le a$, yet there are $a+1$ distinct roots $0, 1, ..., a$, contradiction. Thus, $S_{-c}, ..., S_{-c+a}$ uniquely define one solution, that being $f(t) = 0$ for all $t = 0, 1,..., a$. Yet this means $f(0) = 0$, so the coefficient of $x^ay^bz^c$ is 0, contradiction. $\square$ We finish by induction on the degree of $P$. The base case is $\deg(P) = 0$ whereby $F$ is a constant function. For the inductive step, by Claim 1 we find that the unique monomial with maximum degree is $Kx^ny^nz^n$. Let $P'(x, y, z) = P(x, y, z) - (-1)^nK(x^2 + y^2 + z^2 - xyz)^n$. The unique term of maximal degree $3n$ in $(x^2 + y^2 + z^2 - xyz)^n$ is $x^ny^nz^n$ as well, hence $\deg(P') < \deg(P)$. Furthermore, $P(x, y, z)$ and $(x^2 + y^2 + z^2 - xyz)^n$ both satisfy the problem conditions, so $P'(x, y, z)$ must too. By the inductive hypothesis, $P'(x, y, z) = F'(x^2 + y^2 + z^2 - xyz)$, so $P(x, y, z) = F'(x^2 + y^2 + z^2 - xyz) + (-1)^nK(x^2 + y^2 + z^2 - xyz)^n$ is a polynomial in $x^2 + y^2 + z^2 - xyz$. $\blacksquare$

WLOG the constant term of $P$ is $0$, otherwise we could shift. Also, it now suffices to show that $A=x^2+y^2+z^2-xyz$ is a divisor of $P$, because if it is, we can simply divide out by it and repeatedly shift by a constant and scale by $A$. Call a polynomial "totally radical" if it satisfies the property in the problem condition. Now consider the polynomial $Q$ such that $Q(x,y,z)=P(x,y,z)$, but $Q$'s domain is restricted to $(x,y,z)$ such that $x^2+y^2+z^2-xyz=0$. Then $Q(x,y,z) + AR(x,y,z)$ is totally radical on its domain. Now assume we fix $z>2$. Let $r$ be the larger root of the quadratic $r^2-rz+1=0$. Note that such an $r$ must always exist for $z>2$ by looking at the discriminant. Since the roots are reciprocals of each other, we must have $r>1$. Now pick a starting triple $(a_1,a_2,z)$ in the domain of $Q$ such that $\frac{a_1}{a_2}\neq r$ and $\frac{a_1}{a_2}\neq \frac{1}{r}$. Let $K=Q(a_1,a_2,z)$ Now, in general, let $a_n=za_{n-1}-a_{n-2}$. Since $Q$ is totally radical, $K = Q(a_3,a_2,z)$. Again applying the totally radical property we get $K=Q(a_3,a_4,z)$. We can continue this to get that in general both $Q(a_{2n+1},a_{2n},z)$ and $Q(a_{2n-1},a_{2n},z)$ are equal to $K$. However, it is well known that the linear recurrence above has the solution $a_n=c_1r^n+c_2r^{-n}$ since the roots are reciprocals. This means that $Q(c_1r^{2n+1}+c_2r^{-2n-1},c_1r^{2n}+c_2r^{-2n},z)=K$. Consider the polynomial $S_z(x,y)$ that consists of the set of terms of $Q$ that have the highest total degree in $x$ and $y$ and replacing the $z$'s with the value of $z$ we fixed. If we pick $n$ large in the previous equation then we get that either the highest total degree is $0$ or $(x-ry)$ is a factor of $S_z(x,y)$. Similarly, we get that $(rx-y)$ is a factor by considering the times when $Q(a_{2n-1},a_{2n},z)$ . This means that $(x-ry)(rx-y)=rx^2-(r^2+1)xy+ry^2=r(x^2+y^2-xyz)$ is a factor of $R_z(x,y)$. Since we could do this for every $z>2$, this implies that $(x^2+y^2-xyz)$ is a factor of the sum of the highest degree terms in $Q$. So we could subtract off a multiple of $A$ and be left with a polynomial $Q'$ with a lower maximal degree in $x$ and $y$. Since we restricted the domain of $Q$, $Q'$ is also totally radical. In this way we can repeat the argument until we get a totally radical $Q^*(x,y,z)$ which can be written as a polynomial in $z$. I now claim $Q^*$ must be constant. To do this we can simply vary $(x,y)$ along the solution curve to $x^2+y^2+z^2-xyz=0$. Implicit differentiation means that if we nudge $x$ by $h$ small then $y$ must be nudged by $-\frac{2x-yz}{2y-yz}h$. It is easy to check that under these nudges that $z$ can be chosen so that the derivative of $xz-y$ is nonzero. This means that under the transformation of $(x,y,z)\rightarrow(x,y,xz-y)$, we can find an interval such that $Q^*(x,y,z)$ is constant. This means $Q^*$ is constant. Since we assumed earlier the constant term of $Q$ was $0$, then $Q^*=0$. Since we only subtracted off multiples of $A$ to get to $Q^*$, $P$ must have been divisible by $A$.

Very Nice Problem! Rough Sketch: 1. use the substitution $(2\cos(a),2\cos(b),2\cos(c))$ motivated by identity $a+b+c=0 \implies \cos^2 (a) + \cos^2 (b) + \cos^2 (c) = 2 \cos (a) \cos (b) \cos (c) +1$ 2. Let $R(a,b,c)=P(2\cos(a),2\cos(b),2\cos(c))$ show that it is constant over $a+b+c=0$ 3. Show that $a^2+b^2+c^2-abc-4 | P(a,b,c)-t$ where $t$ is the constant value by division algorithm.

First we prove that $P$ is a polynomial in $x,y, z(xy - z)$ (and you can similarly show it's also a polynomial in $x,y(xz - y), z$ and $x(yz - x), y, z$). It's clear that any polynomial satisfying these three conditions satisfies the problem conditions. Let $Q(x,y,z) = P\left(x,y,z + \frac{xy}{2} \right)$, we have $Q(x,y,z) = Q(x,y,-z) $. Viewing this as a polynomial in $z$, we see that it's even, so all its terms have an even exponent of $z$. Hence $Q$ can be written as a polynomial in $x, y, z^2$, so $P$ can be written as a polynomial in $x, y, \left( z - \frac{xy}{2} \right)^2$, simplifying this gives $x, y, z(z - xy) + \frac{x^2y^2}{4}$, but we can ignore the $\frac{x^2y^2}{4}$ as it is already a polynomial in $x,y$. Therefore, $P(x,y,z)$ can be written as a real polynomial in $x, y, z(z-xy)$. Now, consider some largest degree term of $P$. Let $a, b, c$ be the exponents of $x,y,z$ in this term. For some $x_1, x_2, x_3$, we have $x^a y^b z^c$ is a term of $ x^{x_1} y^{x_2} (z^2 - xyz)^{x_3}$. However, the largest degree term of this is $x^{x_1} y^{x_2} (xyz)^{x_3}$, which is $x^{x_1 + x_3} y^{x_2 + x_3} z^{x_3}$. Thus, $c = x_3$, $a = x_1 + x_3, b = x_2 + x_3$ and so, $a + b \ge 2c$. Similarly, we can use the fact that $P$ is a polynomial in $x,y(xz - y), z$ and in $x(yz - x) , y, z$, to receive that $a + c \ge 2b$ and $b+c \ge 2a$, respectively. Therefore, $2(a+b+c) = (a+b) + (a+c) + (b+c) \ge 2c + 2b + 2a = 2(a+b+c)$, so equality must hold at all points. Hence $a + b = 2c, b + c = 2a, c + a = 2b$. Hence subtracting the first two equations gives that $a - c = 2c - 2a$, so $a = c$. But then $c + a = 2a =2b$, so $a = b = c$. Hence the largest degree term of $P$ is $(xyz)^c$ for some nonnegative integer $c$ (and thus must be unique). Suppose some $P$ satisfied the problem conditions and could not be written as a polynomial in $x^2 + y^2 + z^2 - xyz$. Since any constant polynomial $P$ can be written in this form, choose such $P$ of minimal degree. Let $a$ be the leading coefficient of the largest $(xyz)^c$ term in $P$. Consider the polynomial $Q = P \pm a \cdot (x^2 + y^2 + z^2 - xyz)^c$, where the $\pm$ is $+$ if $c$ is odd, and $-$ if $c$ is even. Notice that $x^2 + y^2 + z^2 - xyz$ can be written as a polynomial in $x,y,z(xy - z)$, as well as $x, y(xz - y) , z$, and $x(yz - x), y, z$, meaning that $Q$ satisfies the problem conditions. Now, $P = Q \pm a \cdot (x^2 + y^2 + z^2 - xyz )^c $ and $Q$ can be written as a polynomial in $x^2 + y^2 + z^2 - xyz$, so $P$ can also be written as a polynomial in $x^2 +y^2 + z^2 - xyz$, contradiction. Hence all polynomials satisfying the problem conditions can be written in the form $F(x^2 + y^2 + z^2 - xyz)$ for some polynomial $F$.

Here's a deceptively simple solution? Note that we can express \[ P(x, y, z) = F(x^2 + y^2 + z^2 - xyz) + P(x, y) + R(y, z) + Q(z, x). \]by repeatedly dividing out any monomial which is divisible by $xyz$ (which also decreases the monovariant of the number of maximal degree terms divisible by $xyz$). Since $F(x^2 + y^2 + z^2 - xyz)$ always satisfies the conditions, it remains to show that $T(x, y, z) \coloneq P(x, y) + R(y, z) + Q(z, x)$ is constant, which also satisfies the given conditions. Note that this implies that \[ P(x, y) + R(y, z) = P(x, xz - y) + R(xz - y, z). \] Claim: This holds for any $xz = w$, ie $P(x, y) + R(y, z) = P(x, w - y) + R(w - y, z)$ holds for any $w$. Proof. Pretend we replace $xz$ with $w$ in the RHS. The map from $x^az^b$ to $w^a x^b z^c$ where one of $b, c$ is zero is injective on our polynomial, so this remains true by considering coefficients in the expansion. This means that our polynomial $T(x, y, z)$ is fixed under shifts of $y$ and the same for $x, z$ and is thus constant.