If $f(x)\in\mathbb Q[x]$ and $f(f(x)),f(f(f(x)))\in\mathbb Z[x]$ prove that $f(x)\in\mathbb Z[x]$ 
If a polynomial $f(x)\in\mathbb Q[x]$ and $$f(f(x)),f(f(f(x)))\in\mathbb Z[x]$$
prove that $f(x)\in\mathbb Z[x]$.

This problem is from 2019 Japan Mathematical Olympiad: see problem 9 (the third day last problem) https://www.imojp.org/archive/mo2019/selection_camp/problems/2019.pdf
Here is what I tried:: Let $f(x)=\dfrac{P(x)}{d}$ where $P(x)\in\mathbb Z[x],d$ is the common denominator of the coefficients of $f(x)$. If let $P(x)=\sum_{i=0}^{n}a_{i}x^i,a_{i}\in\mathbb Z$, then we have
$$f(f(x))=\dfrac{1}{d}\sum_{i=0}^{n}\dfrac{a_{i}}{d^i}\left(\sum_{i=0}^{n}a_{i}x^i\right)^i\in\mathbb Z[x]$$
and $f(f(f(x)))$ is very ugly, so how to solve this problem?
 A: $\newcommand{\d}[1]{\mathtt{\underline{\large{#1}}}}$ This answer can be found in the Math Overflow post by David Speyer here. I give him due credit for it in the comments. I wish to bring to the notice of everybody the arguments used here, so they may attempt to use it in the future. I request feedback from people who wish for me to elaborate on further things.
The Newton Polygon of a $p$-adic polynomial
The Newton polygon of a polynomial has many uses : probably the most relatable use is to prove that the polynomial is irreducible if one can locate the roots of the polynomial. Please read the answer by Captain Lama here for more details on the Newton polygon and its uses.
Nevertheless, we must extend the notion provided in the linked post. For this purpose, I will require knowledge of the following notions.

*

*The $p$-adic numbers $\mathbb Q_p$. Details can be found on the Wikipedia page here.


*The $p$-adic integers $\mathbb Z_p$, which is the ring of integers of $\mathbb Q_p$, and to be fair is not a stretch away from $\mathbb Q_p$ in terms of the definition. See the Wolfram Mathworld post here.


*The $p$-adic valuation $\nu$ , which may have a couple of different definitions, but in the case of this answer is defined as the biggest integer $e$ such that $p^e$ divides that number. Basically, if we express the number as a traditional $p$-adic sum, then the smallest non-zero index of the sum is the valuation. There is some confusion regarding this since an alternate definition is also known : you can see here for some details. I will use the given definition.


*The notion of an algebraic closure of a field, which $\mathbb Q_p$ is. One can see here and here for details. We call the (unique up to isomorphism) algebraic closure as $\mathbb Q_p^{alg}$.


*The fact that the $p$-adic valuation extends to a valuation on $\mathbb Q_p^{alg}$. This is explained by David and many others here.
Before going into things , some motivation. The reason we wish to go into $\mathbb Q_p$ is because, by a coefficient-to-coefficient mapping ,it is easy to figure out a map from $\mathbb Z[x]$ to $\mathbb Z_p[x]$. Furthermore, note that if a polynomial $g(x) \in \mathbb Z_p[x]$ for all primes $p$ then $g(x) \in \mathbb Z[x]$, since one can see the contrapositive by taking the LCM of denominators or any prime dividing the denominator.

Definition : Let $f = f_0 + f_1x + ... + f_dx^d$  be a polynomial with coefficients in $\mathbb Q_p$. The Newton polygon of $f$ is the lower convex hull of the points $(0,\nu(f_0)), ..., (d,\nu(f_d))$. That is, it is the unique piecewise straight line path , such that each of these points either lies above the path or on the path, the points corresponding to $0$ and $d$ must lie on the path, and the path can be non-linear only at the given points. Examples can be found on the Wiki page here.

Now, the Newton polygon has slopes associated to it, which are the slopes of the lines forming it. Suppose that the polygon passes through $(j,N_j)$ for $0 \leq j \leq d$. The quantities $s_j = N_j - N_{j-1}$ are called the slopes of the polygon, and by convexity, $s_1 \leq s_2 \leq ... \leq s_d$. (Note : If one wishes to try , please calculate the slope sequence of $P(x) = 1+5x + \frac 15 x^2 + 35x^3$ with respect to the $5$-adic valuation : the Newton polygon is on the Wiki page. I will hide the answer below.)
Note that a slope sequence , along with the location of the left and rightmost points, uniquely determine the polygon.
Before moving , one notation : $\mathbb Q_p[x]$ is the set of all polynomials having coefficients in $\mathbb Q_p$, and likewise $\mathbb Z_p[x]$.
Two results concerning the Newton Polygon
What's beautiful about the Newton polygon is the fact that you can figure out the Newton polygon of a product of polynomials, from the Newton polygon of the individual polynomials. The reason why this is true is that the coefficients of the product can be easily described in terms of the coefficients of the individual polynomials, and then one just needs to do a little bit of hard work figuring out the details. The details may be found in the paper linked here.

Ostrowski's lemma (OL) : Let $f,f'$ be $p$-adic polynomials with slope sequences $s_1,...,s_n$ and $t_1,...,t_m$ respectively. Then the slope sequence of $ff'$ is the sequence $s_1,...,s_n,t_1,...,t_m$ sorted in ascending order.

The next fact links the roots of $f$, which will lie in the algebraic closure $\mathbb Q_p^{alg}$. Recall that these roots have a valuation. The next result is a genius result. There is a proof in Matthew Baker's notes, and an otherwise excellent explanation of another result, in his blog post here.

The fundamental theorem of the Newton polygon (FNTP) : Let $f$ be a $p$-adic polynomial and $\theta_1,...,\theta_d$ be the roots which lie in $\mathbb Q_p^{alg}$. Then , up to a reordering of the roots, we have $\nu(\theta_j) = -s_j$.

Truly astonishing : this links the $p$-adic values of the coefficients to the $p$-adic values of the roots of the polynomial, via the slopes of the Newton polygon. You can see why this result is so amazing : to link coefficients to roots , all you need is the Newton polygon!
All this serves to illustrate the really powerful machinery we are working with. Indeed, if I were to provide a proof sketch now of the theorem stated by the OP, then I would say this : using the Newton polygon of $f$, we can compute part of the Newton polygon of $f \circ f$. The answer then uses the link between the roots of $f$ and $f \circ f$ , and then from the FTNP transfers to the link between the coefficients of these polynomials, which is what we want to study. So a condition on the coefficients of $f$ becomes a condition on the coefficients of $f \circ f$ and so on.
The main lemma
The main lemma proved by David is the following :

Let $f \in \mathbb Q_p[x]$ not have all its coefficients in $\mathbb Z_p$ , but have the constant term $f_0 \in \mathbb Z_p$. Then $f \circ f$ cannot have all its coefficients in $\mathbb Z_p[x]$.

Proof : Note that when we write $f = f_0 +... + f_dx^d$, we get $f \circ f = f_0 + ... + f_d^{d+1}x^{d^2}$. So the leading term is $f_d^{d+1}$. Clearly if $f_d \notin \mathbb Z_p$ then $f_d^{d+1} \notin \mathbb Z_p$, and we are done. Therefore we assume that $f_d \in \mathbb Z_p$.
What can we infer about the Newton polynomial of $f$? Well, $\nu(f_0) , \nu(f_d) \geq 0$, but by assumption one of $f$'s coefficients is not an integer, so $\nu(f_j) < 0$ for some $j$. Thus, the Newton polygon of $f$ strictly dips below zero to accommodate $(j,\nu(f_j))$, and then recovers. This means that the polygon slopes are first negative and then positive : we can order them as $s_1 \leq s_2 \leq ... \leq s_k < 0 < s_{k+1} \leq ... \leq s_d$. Then $(k,N_k)$ is the lowermost point of the Newton polygon. Note that $N_k <0$  since it is the lowermost point and we know that the polygon comes below zero. We keep track of this fact.
Using FTNC, let us get the roots of $f$, labelled so that $\nu(\theta_j) = -s_j$. Make now a very key observation, which links back to the proof sketch I gave :
$$
f(x) = f_d \prod_{j} (x-\theta_j) \implies f\circ f(x) = f_d \prod_{j} (f(x)-\theta_j)
$$
and therefore $f \circ f(x)$ has been expressed as a product of polynomials, whose Newton polygons can be ascertained. We are firmly in OL territory. However, we must break into two cases, for $j$ being on the left and right hand side of $k$ on the polygon.
Let me first explain what we are trying to achieve in the next few lines : essentially, $f \circ f \notin \mathbb Z_p[x]$ means that the Newton polygon of $f \circ f$ dips below the $x$-axis somewhere, since wherever it dips, the coefficient can't be a $p$-adic integer. So we need to study height changes among the slopes of the polygons of $f(x)-\theta_j$ and transfer them to that of $f \circ f$ using OL.

*

*If $1 \leq j \leq k$, then $\nu(\theta_j) \geq 0$. The constant term of $f - \theta_j$ is $f_0 - \theta_j$ , which by the ultrametric nature of the valuation will have a non-negative valuation. The other coefficients don't change. In particular, the increasing part of the slope sequence of $f$ is retained by $f - \theta_j$, and remains $(s_{k+1},...,s_d)$. Thus, the height change of the Newton polygon of $f-\theta_j$ between  $(k,N_k)$ and $(d,N_d)$ is $s_{k+1} + ... + s_d = N_d - N_k$.


*If $k+1 \leq j \leq d$ then $\nu(\theta_j) < 0$. The constant term of $f-\theta_j$ is $\nu(f - \theta_j) = \nu(\theta_j)  = -s_j< 0$ by the ultrametric property of the valuation. The right hand end point is $(d,N_d)$. Thus, the height change of the Newton polygon of $f-\theta_j$ between $(k,N_k)$ and $(d,N_d)$ is at least $N_d + s_j$, since the polygon will dip at $(0,-s_j)$ and then go up. It turns out the estimate is sufficient.
Now, we are trying to find the dip of $f \circ f$ from the right endpoint, and if this dip is more than the right endpoint's y coordinate, we are done. By OL, we know that all the slopes of $f(x)-\theta_j$ are to be sorted into increasing order, and for finding the lowest point of the Newton polygon of $f \circ f$ we need to find the non-negative slopes and add them together to get the total drop.
What are those non-negative slopes? Well, there's those from $1 \leq j \leq k$, giving a total drop by $N_d - N_k$, and there are $k$ of those, giving a combined drop of $k(N_d - N_k)$ from these cases. For $k+1 \leq j \leq d$, we have a drop of at least $N_d + s_j$ in each case. Finally the total drop is at least $k(N_d - N_k) + \sum_{j=k+1}^d (N_d + s_j)$.
Recall that the leading term of $f \circ f$ is $f_d^{d+1}$. Using the multiplicative property of the valuation , we get that the right endpoint of the Newton polygon of $f \circ f$ is $(d^2, \nu(f_d^{d+1})) = (d^2,N_d(d+1))$.
So, we are left to ask ourselves : is $k(N_d - N_k) + \sum_{j=k+1}^d (N_d + s_j) > N_d(d+1)$?

Exercise you can give your children : Provide them the context, and let them prove that $$k(N_d - N_k) + \sum_{j=k+1}^d (N_d + s_j) - N_d(d+1) = -(k+1)N_k$$

which is a positive quantity since $N_k < 0$. Therefore, the dip is more than the right endpoint! So the Newton polygon of $f_j$ goes below the $x$ axis at some point : and exhale.
The main theorem
The main theorem can be stated and proved easily now.

Let $g$ have coefficients in $\mathbb Q_p$ and let $g \circ g , g \circ g \circ g $ have coefficients in $\mathbb Z_p$. Then $g$ has coefficients in $\mathbb Z_p$.

Proof : Note that $g(g(0)), g(g(g(0)))  \in \mathbb Z_p$, these are the constant terms of the respective polynomials. Let $f(x) = g(x + g(g(0))) - g(g(0))$. Now, note that $$
f \circ f(x) = g(g(x + g(g(0)))) - g(g(0))  = g \circ g( x + g \circ g(0)) - g \circ g(0)
$$
is a difference of polynomials in $\mathbb Z_p[x]$ hence belongs in $\mathbb Z_p[x]$. Finally, $f(0) = g \circ g \circ g(0) - g \circ g(0)$, hence $f(0) \in \mathbb Z_p$. By the contrapositive of the main lemma , we get that $f$ has coefficients in $Z_p$ (if not, then the main lemma would be contradicted).
But then $g$ has coefficients in $\mathbb Z_p$ too, from the following logic : if $f$ has coefficients in $\mathbb Z_p$,then so does $f - g(g(0))$. But then this is just a $g(g(0))$ shift of the polynomial $g(x)$, and any $\mathbb Z_p$ integer shift retains coefficient membership in $\mathbb Z_p$ from a simple binomial theorem type argument and the ultrametric nature of the valuation.
Thus, we are done. But we can go one step further. Define $h^{(r)}$ to be the $r$-fold composition of $h$, with $h^{(1)} = h$.

Let $h \in \mathbb Q_p[x]$ and $k_1,k_2$ be relatively prime such that $h^{k_1},h^{k_2} \in \mathbb Z_p[x]$ for relatively prime $k_1,k_2$. Then $h \in \mathbb Z_p[x]$.

Proof : Compositions retain membership of coefficients in $\mathbb Z_p$, so from Bezout's lemma with $k_1$ and $k_2$, we conclude that if $S = \{r : h^r \notin Z_p[x]\}$ then $S$ is bounded above. Let $r  = \max S$. But then from the previous theorem, $2r,3r \notin S$ so $r \notin S$ as well, a contradiction. Thus, $S$ must be empty i.e. $1 \notin S$, as desired.
Back to our problem
We haven't addressed the OP's problem yet! We will do that easily.

Let $f(f(x)) , f(f(f(x))) \in \mathbb Z[x]$. Then $f \in \mathbb Z[x]$.

Proof : Well, $f(f(x)) , f(f(f(x))) \in \mathbb Z_p[x]$ for all primes $p$, so $f(x) \in \mathbb Z_p[x]$ for all primes $p$, and hence $f(x) \in \mathbb Z[x]$.

Exercise you can give your children : if for co-prime $k_1,k_2$ we have $f^{(k_1)},f^{(k_2)} \in \mathbb Z[x]$ then $f \in \mathbb Z[x]$.

Well , that was worth it. There's another extension we can make as well.

Let $k_1,...,k_n$ be positive numbers such that $\gcd(k_1,...,k_n) = 1$. Then if $h^{(k_i)} \in \mathbb Z[x]$ for all $x$ then $h \in \mathbb Z[x]$.

Proof : Let $S = \{k : h^{(k)} \in \mathbb Z[x]\}$. Since each $k_i \in S$ we know that $k_i$ is closed under positive linear combinations. Thus two consecutive numbers will be contained in $S$ by a generalization of Bezout's lemma. Then we apply the child case, and we are done!
This result helps , for example, if none of the $k_i$ are pairwise coprime : for example, if $h^{(6)}, h^{(15)} , h^{(21)} \in \mathbb Z[x]$ then $h \in \mathbb Z[x]$.
Conclusions
The theory of Newton polygons is useful, since the fundamental lemma helps link the coefficients and roots of a polynomial together, and then OL helps in putting together the Newton polynomial for product polynomials. Aside from what we've done here, the theory can be used to prove a large number of irreducibility results and can give information about the nature of the Galois groups of the extension formed by the splitting field of polynomials over fields, see Baker's blog for the same.
But somewhere in Japan, I am sure they have a five-line solution waiting for me.
