Skip to content

Instantly share code, notes, and snippets.

@pjt33
Last active October 4, 2021 12:37
Show Gist options
  • Save pjt33/c942040a9b8fa9bacc7f5d1bccc7ac2b to your computer and use it in GitHub Desktop.
Save pjt33/c942040a9b8fa9bacc7f5d1bccc7ac2b to your computer and use it in GitHub Desktop.
Interpolating integer polynomials with small coefficients
#!/usr/bin/pypy3
from heapq import heappush, heappop
def mod_inv(a: int, m: int) -> int:
lastremainder, remainder = a % m, m
x, lastx, y, lasty = 0, 1, 1, 0
while remainder:
lastremainder, (quotient, remainder) = remainder, divmod(lastremainder, remainder)
x, lastx = lastx - quotient*x, x
y, lasty = lasty - quotient*y, y
if lastremainder != 1:
raise ValueError
return lastx % m
def chinese_remainder(residues) -> int:
total = 0
M = 1
for m in residues:
M *= m
for m, a in residues.items():
p = M // m
total += a * mod_inv(p, m) * p
return total % M
def eval_poly(ps, x: int) -> int:
y = 0
for c in ps[-1::-1]:
y = x * y + c
return y
def poly_search(points, eval_cost, stop_on_optimal = True, max_degree = None, max_abs_coeff = None):
"""
Search for the integer polynomial with smallest cost which interpolates the points.
Assume that the x-values are coprime.
Assume that the cost evaluation is monotonic in the absolute value of each coefficient.
points: A dict from x to p(x)
eval_cost: A function to evaluate the cost of a coefficient vector, passed as a tuple
stop_on_optimal (optional): set to False to keep producing solutions after finding all of the optimal ones
max_degree (optional): the maximum degree of polynomial to consider
max_abs_coeff (optional): the maximum absolute value of coefficients to consider
"""
max_degree = float('+inf') if max_degree is None else max_degree
max_abs_coeff = float('+inf') if max_abs_coeff is None else max_abs_coeff
lcm = 1
for x in points:
lcm *= x
# (cost, len, coeffs in little-endian order, increment for the last one)
# len isn't strictly necessary, but ensures that infinite chains with fixed L^infty norm
# don't entirely suppress output.
q = [(0, 0, tuple(), None)]
push = lambda child, delta: heappush(q, (eval_cost(child), len(child), child, delta)) if abs(child[-1]) <= max_abs_coeff else None
best_cost = float('+inf')
while q:
cost, _, coeffs, delta = heappop(q)
if cost > best_cost and stop_on_optimal:
return
if len(coeffs) > max_degree + 1:
continue
if all(eval_poly(coeffs, x) == y for x, y in points.items()):
print(coeffs, "at cost", cost)
best_cost = cost
continue
if delta:
push(coeffs[:-1] + (coeffs[-1] + delta,), delta)
n = len(coeffs)
an = chinese_remainder({ x : (y - eval_poly(coeffs, x)) // x**n for x, y in points.items() })
push(coeffs + (an,), lcm)
push(coeffs + (an - lcm,), -lcm)
print("sum of absolute values")
poly_search({3: 221157, 5: 31511625}, lambda coeffs: sum(abs(c) for c in coeffs))
print()
print("maximum absolute value")
poly_search({3: 221157, 5: 31511625}, lambda coeffs: max(abs(c) for c in coeffs))
print()
print("maximum absolute value: extra solutions at max abs coeff 8")
poly_search({3: 221157, 5: 31511625}, lambda coeffs: max(abs(c) for c in coeffs), False, max_degree = 24, max_abs_coeff = 8)
print()
print("maximum absolute value: extra solutions at degree 10")
poly_search({3: 221157, 5: 31511625}, lambda coeffs: max(abs(c) for c in coeffs), False, max_degree = 10, max_abs_coeff = 13)
@chervov
Copy link

chervov commented Oct 3, 2021

Thanks again for very nice idea and code !

Here are some issues occured:

1

For [2:75776, 3: 24659883] I am getting error:
an = chinese_remainder({ x : (y - eval_poly(coeffs, x)) // x**n for x, y in points.items() })
OverflowError: int too large to convert to float

2

For dict_res = {2:1184, 3:33827 }
It quickly finds FIRST polynom, but after waiting couple minutes I stopped it, because it was afraid of falling in infinite loop,
can you please check it ?
p^9+2p^8+3p^5+3p^4+2p^3-p^2+p+2

@pjt33
Copy link
Author

pjt33 commented Oct 3, 2021

1

For [2:75776, 3: 24659883] I am getting error: an = chinese_remainder({ x : (y - eval_poly(coeffs, x)) // x**n for x, y in points.items() }) OverflowError: int too large to convert to float

There shouldn't be any floats involved at all. Are you using Python 2? If so, try it with Python 3. Otherwise try changing x**n to pow(x, n). Note that with the maximum absolute value cost function this will loop infinitely unless you supply the max_degree parameter, because it gets in a loop with (0, 0, 0, 0, 0, 0, 2, 1, -1, 2, -3, 2, -2, -2, 0, -1, -2, -2, -2, ...) and (0, 0, 0, 0, 0, 0, 2, 1, -1, 2, 3, -3, -1, -2, 0, -1, -2, -2, -2, ...).

2

For dict_res = {2:1184, 3:33827 } It quickly finds FIRST polynom, but after waiting couple minutes I stopped it, because it was afraid of falling in infinite loop, can you please check it ? p^9+2p^8+3p^5+3p^4+2p^3-p^2+p+2

Yes, in this case it gets in a loop with (2, 1, -1, 2, -3, 2, -2, -2, 0, -1, -2, -2, -2, ...) and (2, 1, -1, 2, 3, -3, -1, -2, 0, -1, -2, -2, -2, ...).

@chervov
Copy link

chervov commented Oct 4, 2021

Thanks for reply ! No I am using Python 3.
May be problem that in my code 75776 , 24659883 were floats (actually integers) ? So example should be [2:75776.0, 3: 24659883.0] .
How that can be cured ? Any hint ?

@pjt33
Copy link
Author

pjt33 commented Oct 4, 2021

May be problem that in my code 75776 , 24659883 were floats (actually integers) ?

Could be. There's a particular gotcha in Python that / always gives a float. If you want integer division you have to use //.

You can coerce a float to an integer as x = int(x).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment