Writing a program to find polynomial $f(x)$ from $f(1)$ and $f(f(1))$
The conditions $f(1) = p$ and $f(p) = q$ imply the following two equations:
$$
\begin{align}
&amp;\sum_{i=0}^n a_i = p, \\
&amp;\sum_{i=0}^n a_i p^i = q.
\end{align}
$$
When $p &lt; 0$ or $p$ is not an integer, there are no solutions. When $p = 0$, the first equation implies that the polynomial must be $0$ and so $q = 0$.
When $p = 1$, the equations imply that $p = q$, and the first equation implies that the solutions are the polynomials $1,x,\ldots,x^n$. From now on, we assume that $p \geq 2$.

Intuitively, the polynomial minimizing $\sum_i a_i$ under the second constraint above is obtained by writing $q$ in base $p$ notation (first $n$ digits) and putting the rest as a coefficient of $a_n$. Denote this solution by $b_0,\ldots,b_n$. If $\sum_i b_i &gt; p$ then there is no solution. If $\sum_i b_i = p$ then we have found the unique solution. From now on, suppose that $\sum_i b_i &lt; p$; this implies that the base $p$ expansion of $q$ has length at most $n+1$.

The given equations imply more constraints. First of all, $q \geq p$. Second, since $p \equiv 1 \pmod{p-1}$, we must have $q \equiv p \pmod{p-1}$. Starting with the polynomial $g(x) = \sum_{i=0}^n b_i x^i$, we can increase $g(1)$ by applying the update $(b'_{i+1},b'_i) \gets (b_{i+1} - \alpha,b_i + p\alpha)$ (when $b_{i+1} \geq \alpha$); this increases $g(1)$ by $\alpha(p-1)$. We can apply these updates systematically, moving the mass all the way to the polynomial $q$. Since $\sum_i b_i &lt; p \leq q$ and $q \equiv p \pmod{p-1}$, somewhere along the way we will find a (unique) polynomial satisfying both requirements above.

How do we perform this process algorithmically? We go over the indices from $b_n$ to $b_1$ in order. At each point, we subtract from $b_i$ the maximal amount $\alpha$ which keeps the invariant $g(1) \leq p$. After processing $b_1$, we should reach the desired polynomial. (We are basically looking at the base $p-1$ expansion of $p - g(1)$.)