# Written in 2025 by Kirill A. Korinsky # # To the extent possible under law, the author(s) have dedicated all # copyright and related and neighboring rights to this software to the # public domain worldwide. # # This software is distributed without any warranty. import math import itertools from typing import Optional, Iterator, Tuple, List class KunerthSolver: """ A unified deterministic solver for quadratic Diophantine equations: y^2 = ax^2 + bx + c. Based on Adolf Kunerth's algorithms from 1878 and 1880. Methods: 1. Modular Case (a == 0) [Kunerth 1878]: - Strategy: Parallel 'Hare and Tortoise' search. - 'Hare' (Fast Path): Linearly checks v=1 (covers >99% of modular roots). - 'Tortoise' (Complete Path): Diagonally enumerates all pairs (t, v) to guarantee solution. 2. General Case (a != 0) [Kunerth 1880]: - Seed: Geometric Square Injection via wavefront enumeration of (alpha, beta). - Parameter p: Convergents of Continued Fractions of sqrt(a). """ def __init__(self, a: int, b: int, c: int): self.a = a self.b = b self.c = c def _check_hauptformel(self, seed: Tuple[int, int, int], param: Tuple[int, int]) -> Optional[Tuple[int, int]]: m, n, r = seed v, w = param # Denominator: v^2 - a*w^2 (Formula X denominator) denom = v*v - self.a*w*w if denom == 0: return None # Constant K = 2am + bn (Derived from Formula VIII numerator terms) K = 2 * self.a * m + self.b * n # Numerator term: w * (2rv + Kw) # Checks if the rational expression yields an integer num_term = w * (2 * r * v + K * w) if num_term % denom != 0: return None term = num_term // denom # Calculate x: nx = m + term rhs = m + term if rhs % n != 0: return None x = rhs // n # Recover y via Formula IX: y = p(x - m/n) - r/n # Simplified to integer arithmetic: y = (v(nx - m) - rw) / nw y_num = v * term - r * w y_den = n * w if y_den == 0: return None if y_num % y_den != 0: return None y = y_num // y_den # Verify if y*y == self.a*x*x + self.b*x + self.c: return x, y if (-y)**2 == self.a*x*x + self.b*x + self.c: return x, -y return None def _get_modular_seed_components(self, t: int) -> Tuple[int, int, int]: num = t*t - self.c den = self.b g = math.gcd(num, den) m, n = num // g, den // g if n < 0: m, n = -m, -n r = n * t return (m, n, r) def _hare_generator(self) -> Iterator[Tuple[int, int]]: t = 0 while True: # Check v = 1 seed = self._get_modular_seed_components(t) res = self._check_hauptformel(seed, (1, 1)) if res: yield res # Check v = -1 res = self._check_hauptformel(seed, (-1, 1)) if res: yield res t += 1 def _tortoise_generator(self) -> Iterator[Tuple[int, int]]: diagonal_sum = 2 while True: # Iterate through diagonal # v cannot be 0. We skip |v|=1 because Hare handles it. for v_abs in range(2, diagonal_sum + 1): t_abs = diagonal_sum - v_abs # Variants for t (+/-) and v (+/-) ts = [0] if t_abs == 0 else [t_abs, -t_abs] vs = [v_abs, -v_abs] for t in ts: seed = self._get_modular_seed_components(t) for v in vs: res = self._check_hauptformel(seed, (v, 1)) if res: yield res diagonal_sum += 1 def _wavefront_pairs(self) -> Iterator[Tuple[int, int]]: root_a = int(math.isqrt(max(0, self.a))) dist = 0 while True: # Generate pairs (alpha, beta) in a wavefront for i in range(-dist, dist + 1): yield (root_a + i, 0 + dist) if dist > 0: yield (root_a + i, 0 - dist) for j in range(-dist + 1, dist): yield (root_a + dist, 0 + j) yield (root_a - dist, 0 + j) dist += 1 def _wavefront_seeds(self) -> Iterator[Tuple[int, int, int]]: for alpha, beta in self._wavefront_pairs(): # Coefficients of the remainder polynomial R(x) (Kunerth 1878, Eq [2]) A_p = self.a - alpha**2 B_p = self.b - 2 * alpha * beta C_p = self.c - beta**2 # Case 1: Linear Remainder (Equation matches format y^2 = (ax+b)^2 + linear) if A_p == 0: if B_p == 0: continue g = math.gcd(C_p, B_p) m, n = -C_p // g, B_p // g if n < 0: m, n = -m, -n r = alpha * m + beta * n if n != 0 and r*r == self.a*m*m + self.b*m*n + self.c*n*n: yield m, n, r continue # Case 2: Quadratic Remainder (Discriminant condition, Kunerth 1878, Eq [3]) delta = B_p**2 - 4 * A_p * C_p if delta < 0: continue s_delta = math.isqrt(delta) if s_delta**2 != delta: continue for sign in (1, -1): num = -B_p + sign * s_delta den = 2 * A_p if den == 0: continue g = math.gcd(num, den) m, n = num // g, den // g if n < 0: m, n = -m, -n r = alpha * m + beta * n if r*r == self.a*m*m + self.b*m*n + self.c*n*n: yield m, n, r def _continued_fraction_params(self) -> Iterator[Tuple[int, int]]: if self.a <= 0: return a0 = int(math.isqrt(self.a)) if a0 * a0 == self.a: return m, d, an = 0, 1, a0 h1, h2 = 1, 0 k1, k2 = 0, 1 while True: h = an * h1 + h2 k = an * k1 + k2 yield (h, k) yield (-h, k) h2, h1 = h1, h k2, k1 = k1, k m = d * an - m d = (self.a - m * m) // d if d == 0: break an = (a0 + m) // d def solve(self) -> Optional[Tuple[int, int]]: if self.a == 0: gen_hare = self._hare_generator() gen_tortoise = self._tortoise_generator() while True: res = next(gen_hare) if res: return res res = next(gen_tortoise) if res: return res else: seed_gen = self._wavefront_seeds() cf_gen = self._continued_fraction_params() params_cache: List[Tuple[int, int]] = [] for limit in itertools.count(1): new_param = None if cf_gen: try: new_param = next(cf_gen) params_cache.append(new_param) except StopIteration: cf_gen = None if not params_cache: return None seed_gen = self._wavefront_seeds() for i, seed in enumerate(seed_gen): if i >= limit: break if i == limit - 1: for param in params_cache: if res := self._check_hauptformel(seed, param): return res elif new_param: if res := self._check_hauptformel(seed, new_param): return res def verify_case(a: int, b: int, c: int, x_expected: Optional[int], y_expected: int): eq_str = f"y^2 = {a}x^2 + {b}x + {c}" print(f"Testing: {eq_str.ljust(30)}...", end=" ") solver = KunerthSolver(a, b, c) result = solver.solve() if result is None: print(f"[FAIL] Exhausted search.") return x_calc, y_calc = result lhs = y_calc ** 2 rhs = a * (x_calc ** 2) + b * x_calc + c if lhs != rhs: print(f"[FAIL] Math Error. {y_calc}^2 != {rhs}") return matches = False if x_expected is None: # Modular if a == 0: mod = abs(b) target = abs(y_expected) % mod calc = abs(y_calc) % mod if calc == target or calc == (mod - target): matches = True else: matches = True elif x_calc == x_expected: matches = True if matches: print(f"OK [Found: x={x_calc}, y={y_calc}]") else: print(f"OK [Found: x={x_calc}, y={y_calc}. Expected: x={x_expected}]") if __name__ == "__main__": print("--- 1. Modular Cases (1878 Paper) ---") # Example 1: Page 341 verify_case(0, 20521, 1124, None, 9758) # Example 2: Page 343 verify_case(0, 21433, 16315, None, 10732) # Example 3: Page 343 verify_case(0, 239, 127, None, 137) # Example 4: Page 343 verify_case(0, 19139, 4163, None, 9179) print("--- 2. General Cases (1880 Paper) ---") # 1. Exempel: Page 358 verify_case(61, 111, -101, -210, 1633) # 2. Exempel: Page 358 verify_case(67, -290, 244, -24, 214) # 3. Exempel: Page 359 verify_case(156, 447, 127, -41, 494) # 4. Exempel: Page 359 verify_case(3245, 296, -5332, 67, 1767) # 5. Exempel: Page 360 verify_case(639, -15, -615, -8, 201) # 6. Exempel: Page 361 verify_case(126, 528, -264, -10, 84) # 7. Exempel: Page 361 verify_case(540, -78, -110, -7, 164) # 8. Exempel: Page 362 verify_case(199, 97, -20, -763, -10760) # 9. Exempel: Page 363 verify_case(124, 37, 7, -787, 8762) # 10. Exempel: Page 363 verify_case(46, -417, -324, 409, -2743) # 11. Exempel: Page 364 verify_case(351, 358, -697, 43103, 807544) # 12. Exempel: Page 365 verify_case(294, 371, 120, 8, 148) # 13. Exempel: Page 366 verify_case(637, 0, 575, 1235, 31170) # 14. Exempel: Page 370 verify_case(194, 63, 479, -83, 1154) # 15. Exempel: Page 371 verify_case(290, 199, -123, -4, 61)