Example: How can you guess a polynomial $p$ if you know that $p(2) = 11$? It is simple: just write 11 in binary format: 1011 and it gives the coefficients: $p(x) = x^3+x+1$. Well, of course, this polynomial is not unique, because $2x^k$ and $x^{k+1}$ give the same value at $p=2$, so for example $2x^2+x+1$, $4x+x+1$ also satisfy the condition, but their coefficients have greater absolute values!

Question 1: Assume we want to find $q(x)$ with integer coefficients, given its values at some set of primes $q(p_i)=y_i$ such that $q(x)$ has the least possible coefficients. How should we do it? Any suggestion/algorithm/software are welcome. (Least coefficients means: the least maximum of modulii of coefficients).

Question 2: Can one help to guess the polynomial $p$ such that $p(3) = 221157$, $p(5) = 31511625$ with the smallest possible integer coefficients? (Least maximum of modulii of coefficients). Does it exist?

(That example comes from the question MO404817 on count of 3x3 anticommuting matrices $AB+BA=0$ over $F_p$.)

(The degree of the polynomial seems to be 10 or 11. It seems divisible by $x^3$, and I have run a brute force search bounding absolute values of the coefficients by 3, but no polynomial satisfying these conditions is found, so I will increase the bound on the coefficients and will run this search again, but the execution time grows too quickly as the bound increases and it might be that brute force is not a good choice).

Question 3: Do conditions like $q(p_i)=y_i$ imply some bounds on coefficients? E.g., can we estimate that the coefficients are higher than some bound?

  • 1
    $\begingroup$ So no restrictions for the degree? $\endgroup$
    – Pablo H
    Sep 28 at 13:39
  • $\begingroup$ @PabloH No. But too high degree cannot be compatible with small integer coefs and fixed values p(x_i)=y_i, because if you would have x^{very high} it will be bigger than any fixed y_i $\endgroup$ Sep 28 at 13:47

You can certainly do better than brute force by considering modular constraints. If the solution is $p(x) = \sum_i a_i x^i$ then $p(x) - \sum_{j=0}^{n-1} a_j x^j$ is divisible by $x^n$ and $$\frac{p(x) - \sum_{j=0}^{n-1} a_j x^j}{x^n} = a_n \pmod x$$

Solving for $a_0$ in each of the given bases and using the Chinese remainder theorem gives an equivalence class for $a_0$; for each possible value of $a_0$ you can expand similarly for $a_1$; and traversing this tree in order of increasing cost of the coefficients gives a directed search. This works in principle for any cost function which increases when any coefficient increases in absolute value.

This Python code implements the idea and finds two polynomials with sum of absolute values of 29:

$$-2x^3 + 4x^4 + 3x^5 + 11x^6 + x^7 + 5x^8 + 3x^{10} \\ -2x^3 + 4x^4 + 3x^5 - 4x^6 + 9x^7 + 4x^8 + 3x^{10}$$

and one polynomial with maximum absolute value of 7:

$$-2x^3 + 4x^4 + 3x^5 - 4x^6 - 6x^7 - 3x^8 + 7x^9 + 2x^{10}$$

in a small fraction of a second.

Some follow-up questions in comments brought me to the realisation that if we're trying to interpolate $\{ (x_i, y_i) \}$ with the $x_i$ coprime, there is at most one polynomial with coefficients in the range $[-\lfloor \tfrac{(\operatorname{lcm} x_i) - 1}2 \rfloor, \lfloor \tfrac{\operatorname{lcm} x_i}2 \rfloor]$, because the tree collapses to a chain. This gives the following algorithm for finding such a polynomial, if it exists:

M := lcm(x_i)
while any y_i is non-zero:
    find a_0 by Chinese remainder theorem
    if a_0 > floor(M / 2):
        a_0 -= M

    output a_0
    update y_i := (y_i - a_0) / x_i

In the long term, the initial values of the $y_i$ are reduced to negligibility by the repeated division by $x_i$, so eventually each $y_i$ will be reduced to a range which is bounded by $\frac{x_i}{x_i - 1} M$. This means that for a given set of $x_i$ it's possible to compute a finite directed graph to see whether existence is guaranteed.

In the particular case that the $x_i$ are $\{3, 5\}$ there are three cycles, all of them loops: $(0,0) \to (0,0)$ is the terminating loop which indicates that a solution exists, but there are also loops $(2, 1) \to (2, 1)$ and $(-2, -1) \to (-2, -1)$.

  • 2
    $\begingroup$ That's a really attractive piece of Python code $\endgroup$
    – J.J. Green
    Sep 25 at 12:43
  • 2
    $\begingroup$ You are minimizing sum of absolute values; OP asks to minimize largest absolute value. $\endgroup$ Sep 25 at 13:00
  • 1
    $\begingroup$ @GerryMyerson, I misunderstood the original text. The approach should work for a wide class of cost functions including any $L^k$ norm; there may be a risk with the $L^\infty$ norm of building arbitrarily long chains of small norm but in practice it works for the example, as the updated code demonstrates. $\endgroup$ Sep 25 at 14:27
  • 2
    $\begingroup$ @AlexanderChervov, it's stronger than that: this is the only polynomial with coefficients less than 8. $\endgroup$ Sep 25 at 14:34
  • 2
    $\begingroup$ Working mod 15, every equivalence class contains only one number in the range $[-7, 7]$, so the tree is pruned to just a chain. $\endgroup$ Sep 25 at 14:40

Perhaps not the best possible polynomial but a very good polynomial can generally be computed by the LLL algorithm: Let $I$ be an integral polynomial taking the given values at the given points and let $\Lambda$ be the lattice (in $\mathbb Z^{d+1}$ for a suitable high value of the degree $d$, coordinates are of course coefficients) of integral polynomials with roots at the prescribed set of points. We search a short vector in $I+\Lambda$. Adding $I$ to a basis of $\Lambda$ and running the LLL algorithm should give a good polynomial among a basis of short vectors most of the time (It may sometimes happen that it returns only short vectors of the form $kI+P_0$ with $P_0$ in $\Lambda$ and with $k$ avoiding $\pm 1$. Changing the value of $d$ and rerunning LLL solves perhaps the issue).

Added details (after some experimentation): The problem with results giving small multiples of the required evaluations can be avoided by adjoining an additional very large coordinate to the coefficients of $I$ and by setting this coordinate to $0$ for all elements in $\Lambda$. The Maple implementation of LLL returned $$2x^{10}+7x^9-3x^8-6x^7-4x^6+3x^5+4x^4-2x^3$$ for the data of the OP. (I worked with degree $40$.) Unfortunately not a solution with coefficients in $\{-3,\ldots,3\}$ but not very far off.

(Messy maple code for the OP example:

li:={[seq(coeff(t*ii,t^i),i=1..d+2)]};for i from 2 by 1 to d do po:=t*(-t^i+interp([3,5],[3^i,5^i],t)):li:=li union{[seq(coeff(t*po,t^i),i=1..d+2)]};od:


Added after a comment from Max Alekseev His comment is correct ( https://www.goodreads.com/quotes/123557-god-s-final-message-to-his-creation-we-apologize-for-the ): An optimal polynomial (with minimal sum of squared coefficients) is (up to sign) given by a closest lattice point in $I+\Lambda$ to the orthogonal projection of the origin onto the affine hyperplane $I+\mathbb R\otimes_{\mathbb Z}\Lambda$. The trick of adding a last large coordinate to the coordinates of $I$ (or by multiplying $I$ by a large constant $c$ as in Alekseev's answer) moves this affine hyperplane away from the hyperplane containing the codimension one lattice $\Lambda$ and ensures that the last generator of an LLL-basis is in $I+\Lambda$ (up to a sign) and close to the line $\Lambda^\perp$ perpendicular to $\Lambda$ in the ambient spanned by $I$ and $\Lambda$.

  • $\begingroup$ Thank you very much ! Can you share a mapple code please? $\endgroup$ Sep 25 at 12:16
  • $\begingroup$ It looks like you are trying to (re)invent solution to the closest vector problem. $\endgroup$ Sep 25 at 16:28
  • $\begingroup$ @Max Alekseyev Thanks for your insightful remark! $\endgroup$ Sep 26 at 10:13
  • $\begingroup$ LLL was the first thing that came to mind, but this worked example is much better than I could've done. Excellent! $\endgroup$ Sep 27 at 1:34

Given $p(x_i)=y_i$ for $i\in\{1,2,\dots,n\}$, finding a polynomial of fixed degree $d$ can be posed as finding a closest vector to $(cy_1,\dots,cy_n,\underbrace{0,\dots,0}_{d+1})^T$ in the lattice spanned by $(d+1)$ column vectors $$\begin{bmatrix} c\cdot V_{d+1}(x_1,\dots,x_n)\\ I_{d+1} \end{bmatrix}, $$ where $V_{d+1}$ is Vandermonde matrix (of size $n\times (d+1)$), $I_{d+1}$ is the identity matrix, and $c$ is a large constant. By the choice of $c$ we can guarantee that a closest vector will have form $(cy_1,\dots,cy_n,a_0,a_2,\dots,a_d)^T$ and these $a_i$ will be small and deliver us a required polynomial $p(x) = \sum_{i=0}^d a_i x^i$.

As for how to solve CVP in practice, see Package for the Closest Vector Problem (CVP)?


Here is an online SageMath implementation of @PeterTaylor's idea. It uses recursively enumerable sets and MapReduce, and is automatically parallelized by Sage (when available).

As an example it computes all polynomials $p$ of degree at most 12 with coefficients bounded by 10 by absolute value such that $p(3)=221157$ and $p(5)=31511625$.

It also confirms (when run with max_degree = +oo) that there is only one polynomial with the coefficients bounded by $7$ by absolute value.

  • $\begingroup$ Other answers show two different polynomials with coefficients bounded by 7. I did not check them. Somebody must have made a mistake (it could be me). $\endgroup$
    – user132647
    Sep 26 at 5:50
  • 2
    $\begingroup$ @user132647, I gave it in little-endian order and Roland gave it in big-endian order, but it's the same polynomial, and there's a simple proof of uniqueness in the comments on my answer. $\endgroup$ Sep 26 at 6:43
  • $\begingroup$ Oof, now I feel stupid... Thanks for pointing out my mistake! $\endgroup$
    – user132647
    Sep 26 at 6:52

In the $\ell^{\infty}$ norm for the coefficients, the problem of finding polynomial $p(x)=\sum_{k=0}^d c_kx^k$ satisfying $p(x_i)=y_i$ for $i\in\{1,2,\dots,n\}$ can also be naturally posed as an integer linear program for finding integer $c_i$ and the best bound $b$: $$\begin{cases} \textrm{minimize}\ b\\ -b \leq c_k\leq b & \text{for}\ k\in\{0,1,\dots,d\}\\ \sum_{k=0}^d c_k x_i^k = y_i & \text{for}\ i\in\{1,2,\dots,n\} \end{cases} $$

Here is a sample implementation of this approach in Sage.

  • 1
    $\begingroup$ For $d=10$, best_poly([3,5],[221157,31511625],10) gives $2*z^{10} + 7*z^9 - 3*z^8 - 6*z^7 - 4*z^6 + 3*z^5 + 4*z^4 - 2*z^3 + 2*z - 6$ (Best bound: 7), while for $d=20$ best_poly([3,5],[221157,31511625],20) gives $z^{11} + z^{10} - z^9 + z^8 - z^7 + z^6 + z^5 + z^4 + z^3 + z^2 + z + 1$ (Best bound: 1). $\endgroup$
    – Pablo H
    Sep 28 at 13:56
  • $\begingroup$ @PabloH strange $\endgroup$ Sep 28 at 14:07
  • $\begingroup$ @PabloH [ 221980., 56972656.] seems to be different values for the second polynom, than required $\endgroup$ Sep 28 at 14:17
  • 1
    $\begingroup$ @PabloH: That is an issue with the ILP solver (by default SageCell uses GPLK), which sometimes does too much relaxation and inadvertently violates some of the given constraints. This can be remediated by tightening solver parameters or switching to an exact solver (which may be time costly). $\endgroup$ Sep 28 at 17:59
  • 2
    $\begingroup$ @PabloH: I've checked with a locally installed Gurobi solver, and with degree=20 it produces the same polynomial as with degree=10. $\endgroup$ Sep 28 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.