The discriminant of a monic, cubic polynomial
$$p(x) := x^3 + b x^2 + c x + d$$
is the invariant
$$\Delta := -4 b^3 d+b^2 c^2+18 bc d-4 c^3-27 d^2;$$
it has the useful property that $p$ has three distinct, real roots iff $\Delta > 0$. (If $p$ has a repeated real root, then $\Delta = 0$, but in this case the nonrepeated root is always rational.)
On the other hand, by the Rational Root Theorem, if a monic polynomial $x^n + \cdots + d$ has a rational root $r$, then $r$ is an integer that divides $d$. These facts together suggest a way of generating examples:
- Pick a triple $(b, c, d)$ of integers.
- Compute $\Delta$; if $\Delta \leq 0$, $p$ does not have three real roots, so start over and pick a new triple. Otherwise, $p$ has three real roots.
- For each of the factors $s$ of $d$. Computing $p(\pm s)$. If any of these is values is zero, then $p$ has a rational root, so start over and pick a new triple. Otherwise, if none of these is zero, none of the roots of $p$ are rational, that is, $p$ satisfies the condition.
A quick Maple script shows that $2922$ ($31.2\%$) of the $21^3 = 9261$ monic, cubic polynomials with integer coefficients in $-10, \ldots, 10$ satisfy the condition, so the above procedure is efficient in the sense that in practice, one needn't try too many triples $(b, c, d)$ to produce examples.
I'm going to post the method involving trig functions. As an example, we will be using $0=x^3-3x^2+3$.
To use this method, you must know that
$$\cos(3\theta)=4\cos^3(\theta)-3\cos(\theta)$$
which is easily derived from the sum of angles formulas and the Pythagorean theorem.
We manipulate this to give us
$$\cos(3\arccos(r))=4r^3-3r\tag0$$
Anyways, we start with
$$0=ax^3+bx^2+cx+d$$
$$0=x^3-3x^2+3$$
We always start with the substitution $x=y-\frac{b}{3a}$, which for our example is $x=y+1$
$$0=(y+1)^3-3(y+1)^2+3\\=y^3-3y+1$$
In general, we get something along the lines of
$$0=a\left(y-\frac{b}{3a}\right)^3+b\left(y-\frac{b}{3a}\right)^2+c\left(y-\frac{b}{3a}\right)+d\\=ay^3+py+q$$
Where $p$ and $q$ are constants.
Make the substitution $y=uz$ and multiply both sides by $v$.
$$0=v(uz)^3-3v(uz)+3v\\=vu^3z^3-3vuz+v$$
And have $vu^3=4$ and $-3vu=-3$ so that it comes in our $(0)$ form. Solving this system of equations gives $u=2$ and $v=\frac12$.
$$0=4z^3-3z+\frac12$$
$$0=\cos(3\arccos(z))+\frac12$$
Solving for $z$ and recalling the period of cosine, we get
$$z=\cos\left(\frac{2\pi k}3+\frac13\arccos\left(\frac{-1}2\right)\right)$$
$$y=2z=2\cos\left(\frac{2\pi k}3+\frac13\arccos\left(\frac{-1}2\right)\right)$$
$$x=1+y=1+2\cos\left(\frac{2\pi k}3+\frac13\arccos\left(\frac{-1}2\right)\right)$$
$$x=1+2\cos\left(\frac{2\pi k}3+\frac{2\pi}9\right)\qquad k=1,2,3$$
This is not always possible, but under those cases, we can just switch to different trig functions and use their corresponding formulas. It's pretty easy to see which trig function you should use at the "solve for $u$ and $v$ step".
An advantage to this method is that it avoids casus irreducibilis, which is a fancy way of saying you can't factor a cubic polynomial using only real numbers, radicals, and basic arithmetic operations. This does not include trig functions, however, which can be used to obtain nice forms of the solution.
To compare, the radical solution for the above polynomial is given as
$$x_1=1-\frac12(1-i\sqrt3)\sqrt[3]{\frac12(-1+i\sqrt3)}-\frac{1+i\sqrt3}{\sqrt[3]{4(-1+i\sqrt3}}$$
$$x_2=1-\frac{1-i\sqrt3}{\sqrt[3]{4(-1+i\sqrt3)}}-\frac12(1+i\sqrt3)\sqrt[3]{\frac12(-1+i\sqrt3)}$$
$$x_3=1+\frac1{\sqrt[3]{\frac12(-1+i\sqrt3)}}+\sqrt[3]{\frac12(-1+i\sqrt3)}$$
Compared side by side, you might find one form much nicer.
Best Answer
We can go after every such polynomial all at once. The key observation is that the condition $f(x)$ divides $f(x^2)$ is equivalent to the following condition (for an algebraically closed field):
This is proven by completely factoring $f$ as $$f(x)=(x-r_1)^{a_1}(x-r_2)^{a_2}\ldots (x-r_k)^{a_k}$$ for distinct roots $r_i$ with multiplicities $a_i$. If we substitute in $x^2$ for $x$ and note $(x^2-k)=(x-\sqrt{k})(x+\sqrt{k})$ we get a complete factoring of $f^2$: $$f(x^2)=(x-\sqrt{r_1})^{a_1}(x+\sqrt{r_1})^{a_1}(x-\sqrt{r_2})^{a_2}(x+\sqrt{r_2})^{a_2}\ldots (x-\sqrt{r_k})^{a_k}(x+\sqrt{r_k})^{a_k}$$ Note that since the original list of roots was distinct, this list of square roots is also distinct, except if $r_i$ was zero. Since these polynomials are completely factored, $f(x)$ divides $f(x^2)$ if and only if every term $(x-r)^a$ in the factorization of $f(x)$ appears in the factorization of $f(x^2)$ with at least the same multiplicity. Then, noting that the statement we want is trivially true if $x=0$, we are finished.*
We can then spin this into finding every possible such $f$: first, observe that if $x$ is a root, then the sequence $x,x^2,(x^2)^2, ((x^2)^2)^2,\ldots$ must be eventually periodic because all of these must be roots of $f$. This is equivalent to asking that $x$ either be $0$ or be a root of unity.
This can be used to computationally generate every possible polynomial (over $\mathbb C$ - or any field, for that matter) of each degree satisfying the given condition.
There turn out to be a lot of polynomials of this form - though note that every root must have that the sequence of squares $\{x,x^2,x^4,\ldots\}$ has size not exceeding the degree of the polynomial, which ensures that these lists are finite for each degree. For linear terms, you get $$f(x)=x$$ $$f(x)=x-1$$ Since only $1$ and $0$ can be roots. Then, for quadratic terms, you get, letting $\gamma_{a,n}=e^{2\pi i a/n}$ be a root of unity: $$f(x)=x^2$$ $$f(x)=x(x-1)=x^2-x$$ $$f(x)=(x-1)^2=x^2-2x+1$$ $$f(x)=(x-1)(x+1)=x^2-1$$ $$f(x)=(x-\gamma_{1,3})(x-\gamma_{2,3})=x^2+x+1$$ For cubic terms, I'll just list a few of the interesting ones, because you can start combining roots from previous "generations" in numerous uninteresting ways -- note, for instance, that we could take any of the quadratic polynomials, take a square root of any of their roots, and add that as a new root, which would give, already, quite a long list! You could also, multiply any of them by $x$ or $x-1$ to get another example. If we want to look at the "primitive" polynomials which are not divisible by any polynomial in a previous generation, you get the following conjugate pair (neither of which are polynomials over $\mathbb R$): $$f(x)=(x-\gamma_{1,7})(x-\gamma_{2,7})(x-\gamma_{4,7})$$ $$f(x)=(x-\gamma_{3,7})(x-\gamma_{6,7})(x-\gamma_{5,7})$$ For fourth degree, you can expand the list of cubics, similarly. For degree $4$, you get a new real polynomial (which is a cyclotomic polynomial, not coincidentally) and two new complex ones: $$f(x)=(x-\gamma_{1,5})(x-\gamma_{2,5})(x-\gamma_{4,5})(x-\gamma_{3,5})=1+x+x^2+x^3+x^4$$ $$f(x)=(x-\gamma_{1,15})(x-\gamma_{2,15})(x-\gamma_{4,15})(x-\gamma_{8,15})$$ $$f(x)=(x-\gamma_{14,15})(x-\gamma_{13,15})(x-\gamma_{11,15})(x-\gamma_{7,15})$$
I'm fairly sure that you'll get the complete list of polynomials of degree $n$ recursively as follows:
Take the product of any two polynomials already found whose degrees sum to $n$.
Take any polynomial $f$ found in the previous generation and some $r$ such that $r^2$ is a root of $f$ of greater multiplicity than the multiplicity of $r$ (which may be $0$). Multiply $f$ by $(x-r)$.
Let $r$ be a value satisfying $r^{2^n}=r$ and such that no $n'<n$ satisfies this. Take the polynomial $(x-r)(x-r^2)(x-r^4)\ldots(x-r^{2^{n-1}})$.
though I haven't formally examined this. Note that I've only listed the final case for degrees $3$ and $4$ because the first and second cases are extremely numerous.
A stronger statement
characterizes solutions to $f | f\circ g$, proved by similar means and gives similar results for how to list such polynomials.