Interpolating integer polynomials with small coefficients

MO404760.py

#!/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)

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