How can you find an integer coefficient polynomial knowing its values only at a few points (but requiring the coefficients be small)?

Question

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?

Answer 1

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)$.

Answer 2

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:

with(IntegerRelations);
d:=50;
ii:=10^8+t*interp([3,5],[221157,31511625],t);
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:
u:=LLL([seq(li[i],i=1..nops(li))]);
vector(nops(u),i->sum('u[i,j]*3^(j-2)','j'=2..d+2)/221157);
po:=sum('u[nops(u),j]*x^(j-2)','j'=2..d+2);

)

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$.

Answer 3

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)?

Answer 4

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.

Answer 5

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.