Let $C = \{ z \in \mathbb{C} : |z| = 1 \}$ and $D = \{ z \in \mathbb{C} : |z| < 1 \}$ be the unit circle and open unit disk. We will assume $a \ne \pm 1$ as their cases are trivial.
Your statement is true. We are going to prove following generalization:
For $\alpha_1, \alpha_2, \ldots, \alpha_m \in D$, define
$f(z) = \prod\limits_{k=1}^m(z - \alpha_k)$ and $g(z) = \prod\limits_{k=1}^m (1-\bar{\alpha}_k z)$.
The polynomial $f(z) - g(z)$ has all its roots belong to $C$.
Consider their ratio $h(z) = \frac{f(z)}{g(z)}$.
Since all $|\alpha_k| < 1$, $g(z)$ is never zero over $C$ and $h(z)$ is well defined there.
For $z \in C$, we have
$$|h(z)|
= \prod\limits_{k=1}^m \left|\frac{z-\alpha_k}{1-\bar{\alpha}_k z}\right|
= \prod\limits_{k=1}^m \left|\frac{z-\alpha_k}{(\bar{z} - \bar{\alpha}_k)
z}\right|
= 1$$
The ratio $h(z)$ maps $C$ to $C$.
For each factor $\frac{z-\alpha_k}{1-\bar{\alpha}_k z}$, when $z$ move long $C$ once, the factor move along $C$ also once. This implies their product $h(z)$ move along $C$ exactly $m$ times. As a result, we can find $m$ distinct $\theta_1, \ldots, \theta_m \in [ 0, 2\pi )$ such that
$$h(e^{i\theta}) = 1 \iff f(e^{i\theta}) - g(e^{i\theta}) = 0$$
Polynomial $f(z) - g(z)$ has at least $m$ distinct roots over $C$. Since degree of $f(z) - g(z)$ is $m$, counting multiplicity, it has exactly $m$ roots in $\mathbb{C}$. This means above $m$ roots on $C$ is all the roots of $f(z) - g(z)$ and all of them are simple.
On the special case $m = n + 1$ and $(\alpha_1,\alpha_2,\ldots,\alpha_m) = (a,0,\ldots,0)$ where $a \in (-1,1)$.
Polynomial $f(z) - g(z)$ reduces to
$$z^n(z - a) - (1-az) = z^{n+1} - a z^n + az - 1 = P(z)$$
and your statement follows.
IMHO, this is a good chance to introduce the concept of winding number to the students. If they are not ready for that. A standalone proof for the original statement (again $a \ne \pm 1$) goes like this.
When $a \in (-1,1)$, parameterize $C$ by $[0,2\pi) \ni \theta \mapsto z \in C$. We have
$$P(z) = z^{n+1} - az^n + az - 1 = 2ie^{i\frac{(n+1)\theta}{2}}
\left[\sin\frac{(n+1)\theta}{2} - a\sin\frac{(n-1)\theta}{2}\right]$$
Let's call what's inside the square bracket as $I(\theta)$.
When $a$ is real, $I(\theta)$ is clearly real and $\theta = 0$ is a root of it.
Let $\theta_k = \frac{(2k+1)\pi}{n+1}$ for $k = 0,\ldots,n$. When $a \in (-1,1)$, it is easy to see $I(\theta_k)$ is positive for even $k$ and
negative for odd $k$. This means $I(\theta)$ has $n$ more roots. One root
from each interval $(\theta_{k-1},\theta_k)$ for $k = 1,\ldots, n $.
As a result, $I(\theta)$ has at least $n+1$ roots over $[0,2\pi)$. This is equivalent to $P(z)$ has at least $n+1$ roots on $C$. Once again, since $P(z)$ has degree $n+1$, these are all the roots it has.
Best Answer
No need to use derivatives, if you want to continue from your partial attempt.
$(x^2-x)(x^2-x+1) = (x^2-x+\frac{1}{2})^2-\frac{1}{4} = ((x-\frac{1}{2})^2+\frac{1}{4})^2-\frac{1}{4}$.
We can obviously see that it's decreasing for $x < \frac{1}{2}$ and increasing for $x > \frac{1}{2}$.
We certainly know more, since we want to find all $x$ such that:
$((x-\frac{1}{2})^2+\frac{1}{4})^2 = c$
where $c$ is some real.
$(x-\frac{1}{2})^2 = \pm\sqrt{c}-\frac{1}{4}$.
$x = \frac{1}{2} \pm \sqrt{\pm\sqrt{c}-\frac{1}{4}}$. [The signs are independent.]
All roots come in pairs with the outer "$\pm$". So the answer just depends on how many non-real roots. If you want to determine exactly when it has $2$ or $4$ non-real roots, just split into cases based on the $\sqrt{}$, giving the cases $c < 0$, or $0 \le c < \frac{1}{16}$, or $\frac{1}{16} \le c$. You may need to handle the boundaries separately if you don't want to count double roots as separate.