# 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 itertools import math from typing import Iterator, List, Optional, Tuple class KunerthSolver: """ Deterministic solver for quadratic Diophantine equations y^2 = ax^2 + bx + c. Based on Adolf Kunerth's 1877-1880 algorithms, modernized with 2025 optimizations. 1. Modular Case (a == 0): - Shift strategy c' = c + kb (Kunerth 1878) optimized with polynomial sieving. - Solves t^2 = c' via lattice reduction and local Hilbert/Jacobi filters. 2. General Case (a != 0): - Seed Discovery: Wavefront enumeration (Kunerth 1878) replaced/augmented by Gaussian reduction of discriminant forms. - Parameter Generation: Continued fraction convergents of sqrt(a) (Kunerth 1880) accelerated via Baby-Step Giant-Step and NUCOMP. - Integer Conversion: Master Formula (Hauptformel) with Stem Number divisibility constraints. """ FILTER_PRIMES = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31] QR_TABLE = { p: frozenset((i * i) % p for i in range(p)) for p in FILTER_PRIMES } def __init__(self, a: int, b: int, c: int): self.shift = 0 if a != 0: # Normalize equation via vertex shift for a!=0 k = int(math.floor(-b / (2 * a) + 0.5)) if k != 0: self.shift = k c = a * k * k + b * k + c b = 2 * a * k + b self.a = a self.b = b self.c = c # Precompute Discriminant for Modular Stem Filter # S = n^2 * D, where D = b^2 - 4ac self.D = b * b - 4 * a * c self._sieve_cache = [] for p in self.FILTER_PRIMES: try: # Solve k * b = -c (mod p) -> k = -c * b^-1 (mod p) inv_b = pow(self.b, -1, p) self._sieve_cache.append((p, inv_b, False)) except ValueError: # b is a multiple of p; congruence is static self._sieve_cache.append((p, 0, True)) @staticmethod def _jacobi(a: int, n: int) -> int: """Computes Jacobi symbol (a/n).""" if n <= 0 or n % 2 == 0: return 0 a %= n result = 1 while a != 0: while a % 2 == 0: a //= 2 n_mod_8 = n % 8 if n_mod_8 in (3, 5): result = -result a, n = n, a if a % 4 == 3 and n % 4 == 3: result = -result a %= n if n == 1: return result return 0 @staticmethod def _xgcd(a: int, b: int) -> Tuple[int, int, int]: """Extended Euclidean Algorithm.""" prev_x, x = 1, 0 prev_y, y = 0, 1 while b: q = a // b x, prev_x = prev_x - q * x, x y, prev_y = prev_y - q * y, y a, b = b, a % b return a, prev_x, prev_y def _check_hauptformel(self, seed: Tuple[int, int, int], param: Tuple[int, int]) -> Optional[Tuple[int, int]]: """Applies Kunerth's Master Formula (Formula X) with Modular Stem Filtering.""" 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 # Use abs(denom) for modulus operations abs_denom = abs(denom) # Stem Number Check (Kunerth 1880) # S = (2am + bn)^2 - 4ar^2 = K^2 - 4ar^2 # Optimization: Check if S is divisible by denom using modular arithmetic # S = n^2(b^2 - 4ac) + 4a(denom * m^2 + ...) -> simplifies to n^2*D check d_mod = self.D % abs_denom if (pow(n, 2, abs_denom) * d_mod) % abs_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 _can_be_perfect_square(self, n: int) -> bool: if n < 0: return False for p in self.FILTER_PRIMES: if (n % p) not in self.QR_TABLE[p]: return False return True def _try_seed(self, alpha: int, beta: int) -> Optional[Tuple[int, int, int]]: # 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: return None g, _, _ = self._xgcd(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: return (m, n, r) return None # Case 2: Quadratic Remainder (Discriminant condition, Kunerth 1878, Eq [3]) delta = B_p**2 - 4 * A_p * C_p if not self._can_be_perfect_square(delta): return None if delta < 0: return None s_delta = math.isqrt(delta) if s_delta**2 != delta: return None for sign in (1, -1): num = -B_p + sign * s_delta den = 2 * A_p if den == 0: continue g, _, _ = self._xgcd(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: return (m, n, r) return None def _compose_forms(self, Q1: int, P1: int, Q2: int, P2: int) -> Tuple[int, int]: """NUCOMP algorithm for composing quadratic forms (Shanks/van der Poorten).""" g, u, v = self._xgcd(Q1, Q2) limit = int(math.isqrt(math.isqrt(max(1, abs(self.a))))) m = Q2 // g if m == 0: return self._reduce_form(Q1, P1) # Safety for edge cases delta = (P2 - P1) // g x = (u * delta) % m r2, r1 = Q1 // g, x c2, c1 = 0, 1 # Euclidean descent while abs(r1) > limit and r1 != 0: q = r2 // r1 r2, r1 = r1, r2 - q * r1 c2, c1 = c1, c2 - q * c1 m_prime = Q1 // g if abs(r1) > limit: Q3 = (Q1 * Q2) // (g * g) t = (u * Q1 * (P2 - P1)) % Q3 s = (t + g * P1) % Q3 return self._reduce_form(Q3, s) # NUCOMP formulation num = r1 * r1 - self.a * c1 * c1 if num % m_prime != 0: # Fallback if non-exact (should not happen in valid group) Q3 = (Q1 * Q2) // (g * g) else: Q3 = num // m_prime abs_Q3 = abs(Q3) if abs_Q3 == 0: return (0, 0) t = (u * Q1 * delta) % abs_Q3 P3 = (t + g * P1) % abs_Q3 return self._reduce_form(Q3, P3) def _reduce_form(self, Q: int, P: int) -> Tuple[int, int]: while True: P_new = P Q_new = (self.a - P_new * P_new) // Q if abs(Q_new) >= abs(Q) and Q_new != 0: break P = int(math.floor( (-P_new + math.sqrt(self.a)) / abs(Q_new))) * abs(Q_new) - P_new Q = Q_new return Q, P def _compose_pell(self, v1: int, w1: int, v2: int, w2: int) -> Tuple[int, int]: return (v1 * v2 + self.a * w1 * w2, v1 * w2 + w1 * v2) def _lattice_reduction_gen( self, r: int, window_size: int, next_jump_dist: int, a0: Optional[int] = None ) -> Iterator[Tuple[Tuple[int, int], Tuple[int, int]]]: """Generates small-norm candidates via lattice reduction.""" if r <= 0: return if a0 is None: a0 = int(math.isqrt(r)) if a0 * a0 == r: yield (a0, 1), (0, 1) return v_curr, w_curr = 1, 0 h1, k1 = 1, 0 h2, k2 = a0, 1 P, Q = 0, 1 a_n = a0 g_v, g_w = 1, 0 g_P, g_Q = 0, 1 total_steps = 0 stop_coeff = 2 * a0 while True: final_v, final_w = self._compose_pell(g_v, g_w, h2, k2) yield (final_v, final_w), (P, Q) P = a_n * Q - P Q = (r - P * P) // Q if Q == 0: break a_n = (a0 + P) // Q h1, h2 = h2, a_n * h2 + h1 k1, k2 = k2, a_n * k2 + k1 total_steps += 1 if g_v == 1 and g_w == 0 and a_n == stop_coeff: final_v, final_w = self._compose_pell(1, 0, h2, k2) # Flush yield (final_v, final_w), (P, Q) return if total_steps >= next_jump_dist: g_v, g_w = self._compose_pell(g_v, g_w, h2, k2) g_Q, g_P = self._compose_forms(g_Q, g_P, Q, P) P, Q = g_P, g_Q h1, k1 = 1, 0 h2, k2 = 0, 1 a_n = (a0 + P) // Q h1, h2 = 1, a_n k1, k2 = 0, 1 next_jump_dist = window_size def _get_best_rational_candidates(self, limit: int) -> List[Tuple[int, int]]: raw_candidates = [] a1, a2 = self.b, self.c u1, u2 = 0, 1 while a2 != 0: raw_candidates.append((a2, u2)) q = a1 // a2 a1, a2 = a2, a1 - q * a2 u1, u2 = u2, u1 - q * u2 unique_targets = {} abs_c = abs(self.c) for n, d in raw_candidates: nd = n * d abs_nd = abs(nd) if abs_nd == abs_c: continue s = math.isqrt(abs_nd) if s * s == abs_nd: score = -1 elif n == 1: score = 1 else: score = abs_nd if nd not in unique_targets: unique_targets[nd] = (score, n, d) sorted_candidates = sorted(unique_targets.values(), key=lambda x: x[0])[:limit] result = [] # 1. Always check the integer case (c + kb) first g, d_inv, _ = self._xgcd(1, self.b) result.append((self.c, d_inv % self.b)) # 2. Then check rational approximations for _, n, d in sorted_candidates: g, d_inv, _ = self._xgcd(d, self.b) if g == 1: result.append((n * d, d_inv % self.b)) return result def _get_val_unit(self, n: int, p: int) -> Tuple[int, int]: val = 0 if n % p == 0: while n % p == 0: val += 1 n //= p return val, n % p def _compute_local_hilbert(self, p: int, val_c: int, unit_c: int, val_b: int, jac_b: int) -> int: p_eps = (p - 1) // 2 sign = 1 if (p_eps * val_c * val_b) % 2 != 0: sign = -1 if val_b % 2 != 0: if self._jacobi(unit_c, p) == -1: sign = -sign if val_c % 2 != 0: if jac_b == -1: sign = -sign return sign def _prioritized_search_gen( self, targets: List[Tuple[int, int]], k_lim: int, block_size: int) -> Iterator[Tuple[int, int, int]]: """Bytearray-based Quadratic Residue Sieve for shift parameters.""" sieve_data = [] for target_c, _ in targets: target_sieve = [] for p, inv_b, is_degenerate in self._sieve_cache: bad_k_offsets = [] t_c_mod = target_c % p if is_degenerate: if pow(t_c_mod, (p - 1) // 2, p) == p - 1: bad_k_offsets = list(range(p)) else: for r in range(p): if pow(r, (p - 1) // 2, p) == p - 1: k_bad = ((r - t_c_mod) * inv_b) % p bad_k_offsets.append(k_bad) if bad_k_offsets: target_sieve.append((p, bad_k_offsets)) sieve_data.append(target_sieve) # Iterate through blocks for k_start in range(0, k_lim + 1, block_size): k_end = min(k_start + block_size, k_lim + 1) actual_size = k_end - k_start block_results = [] for t_idx, (target_c, d_inv) in enumerate(targets): # Optimization: Use bytearray for fast C-layer slice assignment # 1 = Valid, 0 = Invalid is_valid = bytearray([1]) * actual_size sieve_rules = sieve_data[t_idx] for p, bad_offsets in sieve_rules: if not any(is_valid): break if len(bad_offsets) == p: is_valid[:] = b'\x00' * actual_size break for bad_off in bad_offsets: start_i = (bad_off - k_start) % p if start_i < actual_size: # Sieve stride assignment count = (actual_size - start_i + p - 1) // p is_valid[start_i::p] = b'\x00' * count # Collect survivors for i, valid in enumerate(is_valid): if valid: val = target_c + (k_start + i) * self.b shift = 0 if val == 0: block_results.append( (0, k_start + i, target_c, d_inv)) continue # Count trailing zeros for smoothness score temp = val while temp % 2 == 0: temp //= 2 shift += 1 # (3, 5, 7 mod 8 are never squares). if temp % 8 != 1: continue block_results.append( (val, k_start + i, target_c, d_inv)) # Sort block by smoothness (Spirit: check smooth numbers first) block_results.sort(key=lambda x: x[0]) for _, k, target_c, d_inv in block_results: yield k, target_c, d_inv def _modular_seeds(self, window_size: int, jump_dist: int, rational_limit: int, block_size: int) -> Iterator[Tuple[int, int, int]]: targets = self._get_best_rational_candidates(rational_limit) # We have no idea about b strucutre, so, jacobi will lie k_lim = abs(self.b) t_lim = abs(self.b) # Precompute Hilbert/Jacobi data for b hilbert_cache = {} for p in self.FILTER_PRIMES: val_b, unit_b = self._get_val_unit(self.b, p) jac_b = self._jacobi(unit_b, p) val_nb, unit_nb = self._get_val_unit(-self.b, p) jac_nb = self._jacobi(unit_nb, p) hilbert_cache[p] = (val_b, jac_b, val_nb, jac_nb) # Larger block size for bytearray efficiency optimized_block_size = max(block_size, 65536) # The theoretical bound is strict (2*s_prime), but our 'jump_dist' # means we might skip the exact convergent that minimizes the norm. barrier_buffer = 256 + jump_dist * 4 for k, target_c, d_inv in self._prioritized_search_gen( targets, k_lim, optimized_block_size): c_prime = target_c + k * self.b abs_c_prime = abs(c_prime) if c_prime > 0: s_prime = math.isqrt(abs_c_prime) if s_prime * s_prime == abs_c_prime: y_cand = s_prime % abs(self.b) x = k yield (x, 1, y_cand) else: s_prime = math.isqrt(abs_c_prime) # 2. The Lattice Barrier (Legendre's Bound / Coppersmith-like limit) if abs(self.b) > 2 * s_prime + barrier_buffer: continue jacobi_mod = c_prime while jacobi_mod != 0 and jacobi_mod % 2 == 0: jacobi_mod //= 2 if jacobi_mod == 0: continue allow_pos = False allow_neg = False if jacobi_mod == 1: allow_pos = True allow_neg = True else: if self._jacobi(self.b, jacobi_mod) != -1: allow_pos = True if self._jacobi(-self.b, jacobi_mod) != -1: allow_neg = True if not allow_pos and not allow_neg: continue # Check Hilbert Obstructions is_blocked = False for p in self.FILTER_PRIMES: val_c, unit_c = self._get_val_unit(c_prime, p) val_b, jac_b, val_nb, jac_nb = hilbert_cache[p] if allow_pos: if self._compute_local_hilbert(p, val_c, unit_c, val_b, jac_b) == -1: allow_pos = False if allow_neg: if self._compute_local_hilbert(p, val_c, unit_c, val_nb, jac_nb) == -1: allow_neg = False if not allow_pos and not allow_neg: is_blocked = True break if is_blocked: continue lat_gen = self._lattice_reduction_gen(c_prime, window_size, jump_dist) for param, _ in lat_gen: t, _ = param if abs(t) > t_lim: break diff = t * t - c_prime if diff > 0 and not allow_pos: continue if diff < 0 and not allow_neg: continue y_cand = (t * d_inv) % self.b if (y_cand * y_cand - self.c) % self.b == 0: x = (y_cand * y_cand - self.c) // self.b yield (x, 1, y_cand) def _wavefront_seeds( self, max_dist: Optional[int] = None) -> Iterator[Tuple[int, int, int]]: root_a = int(math.isqrt(max(0, self.a))) for dist in itertools.count(0): if max_dist is not None and dist > max_dist: return for i in range(-dist, dist + 1): alpha = root_a + i if alpha == 0: continue center = int(math.floor(self.b / (2 * alpha) + 0.5)) for beta in ([center + dist, center - dist] if dist > 0 else [center]): seed = self._try_seed(alpha, beta) if seed: yield seed for j in range(-dist + 1, dist): for sign in [1, -1]: alpha = root_a + sign * dist if alpha == 0: continue center = int(math.floor(self.b / (2 * alpha) + 0.5)) seed = self._try_seed(alpha, center + j) if seed: yield seed def _lagrange_reduce_bqf( self, A: int, B: int, C: int) -> Tuple[Tuple[int, int], Tuple[int, int, int]]: p, q = 1, 0 r, s = 0, 1 if abs(C) < abs(A): A, C = C, A B = -B p, q, r, s = q, p, s, r prev_max = None while True: if A == 0: break curr_max = max(abs(A), abs(C)) if prev_max is not None and curr_max >= prev_max: break prev_max = curr_max k = int(math.floor((B + A) / (2 * A))) if k != 0: B = B - 2 * k * A q = q - k * p s = s - k * r if abs(B) <= A and A <= C: break A, C = C, A B = -B p, q = q, p r, s = s, r return (p, r), (A, B, C) def _gaussian_reduced_seeds( self, search_radius: int) -> Iterator[Tuple[int, int, int]]: """Finds seeds via Lagrange reduction of the discriminant form.""" if self.a <= 0: return root_a = int(math.isqrt(self.a)) for base_alpha in [root_a, root_a + 1, root_a - 1] if root_a > 1 else [1, 2]: if base_alpha == 0: continue center_beta = int(math.floor(self.b / (2 * base_alpha) + 0.5)) A_form = 4 * (center_beta**2 + self.c) B_form = 8 * (base_alpha * center_beta - self.b // 2) C_form = 4 * (base_alpha**2 + self.a) if A_form == 0 or C_form == 0: continue try: (p, r), _ = self._lagrange_reduce_bqf(A_form, B_form, C_form) except (ZeroDivisionError, ValueError): continue for dx_red in range(-search_radius, search_radius + 1): for dy_red in range(-search_radius, search_radius + 1): dx = p * dx_red dy = r * dx_red + dy_red alpha = base_alpha + dx beta = center_beta + dy if alpha == 0: continue seed = self._try_seed(alpha, beta) if seed: yield seed def _modular_params(self, v_limit: int) -> Iterator[Tuple[int, int]]: # Previously hardcoded as [1, -1, 2, -2, 3, -3, 4, -4, 5, -5] for v in range(1, v_limit + 1): yield (v, 1) yield (-v, 1) def _multi_norm_lattice_gen(self, max_norm: int, w_limit: int, a0: int) -> Iterator[Tuple[int, int]]: """Direct search for small |v^2 - a*w^2| using pre-computed a0.""" if self.a <= 0: return # We know v ~= w * sqrt(a). # Using integer arithmetic: v ~= w * a0. # This prevents re-calculating isqrt for every w. for w in range(1, min(w_limit, max_norm + 1)): v_approx = w * a0 # Check neighborhood v, v+1 to find minimal norm for v in [v_approx, v_approx + 1]: if v <= 0: continue norm = v * v - self.a * w * w if abs(norm) <= max_norm and norm != 0: yield (v, w) def _pow_pell(self, v: int, w: int, n: int) -> Tuple[int, int]: res_v, res_w = 1, 0 base_v, base_w = v, w while n > 0: if n % 2 == 1: res_v, res_w = self._compose_pell(res_v, res_w, base_v, base_w) base_v, base_w = self._compose_pell(base_v, base_w, base_v, base_w) n //= 2 return res_v, res_w def _bsgs_params(self, window_size: int, jump_dist: int, max_norm: int, a0: int) -> Iterator[Tuple[int, int]]: # Baby Steps using pre-computed a0 gen = self._lattice_reduction_gen(self.a, window_size, window_size + 1, a0) jump_v, jump_w = 0, 0 jump_P, jump_Q = 0, 0 steps_taken = 0 for param, form in gen: v, w = param P, Q = form if abs(Q) <= max_norm and Q != 0: yield (v, w) jump_v, jump_w = v, w jump_P, jump_Q = P, Q steps_taken += 1 if steps_taken >= window_size: break if steps_taken < window_size: return # Giant Steps: traverse form cycle without large arithmetic curr_P, curr_Q = jump_P, jump_Q jump_count = 1 while jump_count < jump_dist: curr_Q, curr_P = self._compose_forms(curr_Q, curr_P, jump_Q, jump_P) jump_count += 1 if abs(curr_Q) <= max_norm and curr_Q != 0: # Lazy recovery of large coefficients final_v, final_w = self._pow_pell(jump_v, jump_w, jump_count) yield (final_v, final_w) def _general_params(self, window_size: int, jump_dist: int, max_norm: int, w_limit: int, a0: int) -> Iterator[Tuple[int, int]]: # For a!=0: Interleave BSGS convergents with small-norm search cf_gen = self._bsgs_params(window_size, jump_dist, max_norm, a0) lattice_gen = self._multi_norm_lattice_gen(max_norm, w_limit, a0) cf_exhausted = False lattice_exhausted = False while not (cf_exhausted and lattice_exhausted): if not cf_exhausted: try: yield next(cf_gen) except StopIteration: cf_exhausted = True if not lattice_exhausted: try: yield next(lattice_gen) except StopIteration: lattice_exhausted = True def _discriminant_matching_seeds( self, a0: int, search_radius: int) -> Iterator[Tuple[int, int, int]]: """ Solves for injection parameters (alpha, beta) by iterating alpha and calculating the required beta to make the discriminant a square. """ root_a = a0 for i in range(-search_radius, search_radius + 1): alpha = root_a + i if alpha == 0: continue A_p = self.a - alpha**2 if A_p == 0: continue delta_0 = self.b**2 - 4 * A_p * self.c # 1. Calculate analytical lower bound for S if self.a > 0: offset_limit = (self.b * alpha)**2 // self.a lower_bound_sq = delta_0 - offset_limit else: lower_bound_sq = 0 # Fallback for non-standard a<=0 start_s = math.isqrt(max(0, lower_bound_sq)) # 2. Define limited dynamic search window based on b s_range = max(2000, abs(self.b) * 10) for s_target in range(start_s, start_s + s_range): target_sq = s_target * s_target CoeffA = 4 * self.a CoeffB = -4 * self.b * alpha CoeffC = delta_0 - target_sq D_beta = CoeffB**2 - 4 * CoeffA * CoeffC if D_beta >= 0: s_beta_D = math.isqrt(D_beta) if s_beta_D * s_beta_D == D_beta: # Beta solutions: (-B +/- sqrt(D)) / 2A for sign in [1, -1]: num = -CoeffB + sign * s_beta_D den = 2 * CoeffA if den != 0 and num % den == 0: beta = num // den seed = self._try_seed(alpha, beta) if seed: yield seed def _attempt_rational_recovery( self, Y_curr: int, z: int, p: int, q: int, Num_Ap: int, target_rem: int) -> Optional[Tuple[int, int, int]]: valid_Y = None if Y_curr % (2 * self.a) == target_rem: valid_Y = Y_curr elif (-Y_curr) % (2 * self.a) == target_rem: valid_Y = -Y_curr if valid_Y is not None: B = (valid_Y + self.b * p) // (2 * self.a) Term_B = self.b * q * q - 2 * p * B Term_A = 2 * Num_Ap num = -Term_B + z den = Term_A if den != 0: g = math.gcd(num, den) m_sol, n_sol = num // g, den // g if n_sol < 0: m_sol, n_sol = -m_sol, -n_sol r_sol = p * m_sol + B * n_sol return (m_sol * q, n_sol * q, r_sol) return None def _convergent_matching_seeds( self, a0: int, limit: int) -> Iterator[Tuple[int, int, int]]: """Uses convergents p/q of sqrt(a) as rational alpha""" if a0 * a0 == self.a: return m, d, a_n = 0, 1, a0 h1, k1 = 1, 0 h2, k2 = a0, 1 D_poly = self.b**2 - 4 * self.a * self.c y_search_window = min(5000, max(200, 2 * self.a)) for _ in range(limit): p, q = h2, k2 Num_Ap = self.a * q * q - p * p N = Num_Ap * D_poly target_rem = (-self.b * p) % (2 * self.a) # Solve Y^2 - a*Z^2 = N for Y if N > 0: Y_start = math.isqrt(N) for offset in range(y_search_window): Y_curr = Y_start + offset val = Y_curr * Y_curr - N if val >= 0 and val % self.a == 0: z_sq = val // self.a z = math.isqrt(z_sq) if z * z == z_sq: res = self._attempt_rational_recovery( Y_curr, z, p, q, Num_Ap, target_rem) if res: yield res else: neg_N = -N for Y_curr in range(y_search_window): val = Y_curr * Y_curr + neg_N if val % self.a == 0: z_sq = val // self.a z = math.isqrt(z_sq) if z * z == z_sq: res = self._attempt_rational_recovery( Y_curr, z, p, q, Num_Ap, target_rem) if res: yield res # CF step m = d * a_n - m d = (self.a - m * m) // d if d == 0: break a_n = (a0 + m) // d h2, h1 = a_n * h2 + h1, h2 k2, k1 = a_n * k2 + k1, k2 def solve(self, max_iterations: Optional[int] = None, window_size: int = 50, jump_dist: int = 100, rational_limit: int = 10, modular_block_size: int = 1024, lattice_max_norm: int = 100, lattice_w_limit: int = 50, reduced_search_radius: int = 5, modular_v_max: int = 5, analytic_search_radius: int = 30) -> Optional[Tuple[int, int]]: # Early Exit Criteria (Kunerth 1880, Achterreste) if self.a == 0: if self.b % 4 == 0 and self.c % 4 in (2, 3): return None if self.b % 8 == 0 and self.c % 8 == 5: return None if self.a == 0: # Modular case: seeds already encode solutions iterations = 0 for seed in self._modular_seeds(window_size, jump_dist, rational_limit, modular_block_size): m, n, r = seed for param in self._modular_params(modular_v_max): res = self._check_hauptformel(seed, param) if res: return res iterations += 1 if max_iterations is not None and iterations >= max_iterations: break else: # General case a0 = int(math.isqrt(abs(self.a))) # Phase 1: Gaussian-reduced seeds (fast for typical cases) for seed in self._gaussian_reduced_seeds(reduced_search_radius): for param in self._general_params(window_size, jump_dist, lattice_max_norm, lattice_w_limit, a0): res = self._check_hauptformel(seed, param) if res: return (res[0] + self.shift, res[1]) # Phase 2: Rational Convergent Seeds for seed in self._convergent_matching_seeds( a0, limit=analytic_search_radius): for param in self._general_params(window_size, jump_dist, lattice_max_norm, lattice_w_limit, a0): res = self._check_hauptformel(seed, param) if res: return (res[0] + self.shift, res[1]) # Phase 3: Integer Analytic Matching for seed in self._discriminant_matching_seeds( a0, search_radius=analytic_search_radius): for param in self._general_params(window_size, jump_dist, lattice_max_norm, lattice_w_limit, a0): res = self._check_hauptformel(seed, param) if res: return (res[0] + self.shift, res[1]) # Phase 4: Wavefront seeds with distance limit for seed in self._wavefront_seeds(max_dist=max_iterations): for param in self._general_params(window_size, jump_dist, lattice_max_norm, lattice_w_limit, a0): res = self._check_hauptformel(seed, param) if res: return (res[0] + self.shift, res[1]) return None if __name__ == "__main__": def abbreviate(n: Optional[int], limit: int = 20) -> str: if n is None: return "None" s = str(n) if len(s) <= limit: return s half = (limit - 3) // 2 return f"{s[:half]}...{s[-half:]}" def verify_case(a: int, b: int, c: int, x_expected: Optional[int], y_expected: int, **kwargs): import time from datetime import timedelta parts = [] if a != 0: parts.append(f"{abbreviate(a)}x^2") if b != 0: parts.append(f"{abbreviate(b)}x") parts.append(f"{abbreviate(c)}") rhs_str = " + ".join(parts).replace("+ -", "- ") print(f"Testing: y^2 = {rhs_str} ") start = time.time() solver = KunerthSolver(a, b, c) result = solver.solve(**kwargs) end = time.time() print(f" ~ {timedelta(seconds=end-start)}") if result is None: if x_expected is None and y_expected == 0: print(f" -> OK [Impossible]") else: 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!") return if a == 0: modular_root = abs(y_calc) % abs(b) status = "OK" note = "" if modular_root == y_expected else " (Found alternative root)" print( f" -> {status}{note} [y: {abbreviate(y_calc)} , Mod: {abbreviate(modular_root)}]" ) else: status = "OK" if x_calc == x_expected else f"OK (Found different solution)" print( f" -> {status} [x: {abbreviate(x_calc)}, y: {abbreviate(y_calc)}]" ) 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) print("--- 3. Impossible Cases ---") # Theoretical Basis: "Kriterien der Unmöglichkeit", 1880, p. 346 # and Modular Exclusion, 1878b, p. 339. # y^2 = 8x + 5. Mod 8 check. verify_case(0, 8, 5, None, 0) # y^2 = 4x + 3. Mod 4 check. verify_case(0, 4, 3, None, 0) # y^2 = 100x + 2. Mod 4 check. verify_case(0, 100, 2, None, 0) # y^2 = 17x^2 + 8x - 1. Provably no solutions (requires factoring to verify). # Use max_iterations to limit search. verify_case(17, 8, -1, None, 0, max_iterations=100) print("--- 4. Large Tests ---") # Test 1: Moderate General Case (a=133), x ~ 10^12 # Checks basic wavefront and CF speed on mid-sized numbers. verify_case(133, 456, -12000000007224000001068060, 1000000000337, 11000000003767) # Test 2: Large General Case (a=991 is a classic "hard" Pell number), x ~ 10^18 # Exercises BSGS / Giant Steps infrastructure. verify_case(991, -123, -30000000000000036249000000000010846602, 1000000000000000662, 31000000000000020576) print("--- 5. Lattice Method Tests ---") # Test cases that benefit from Gaussian reduction / multi-norm lattice methods # Case where alpha is close to sqrt(a) - tests Babai's nearest plane centering verify_case(100, 19, -5, None, 1) # sqrt(100)=10, solution needs alpha near 10 # Case with perfect square a - tests edge case handling verify_case(49, 14, 1, 0, 1) # y^2 = 49x^2 + 14x + 1 = (7x+1)^2 verify_case(0, 17593, pow(123, 2, 17593), None, 123) verify_case(0, 17593, pow(12345, 2, 17593), None, 12345) verify_case(0, 1343021, pow(123, 2, 1343021), None, 123) verify_case(0, 1343021, pow(12345, 2, 1343021), None, 12345) verify_case(0, 6426786497, pow(1234, 2, 6426786497), None, 1234) verify_case(0, 6426786497, pow(12345, 2, 6426786497), None, 12345) verify_case(0, 665059573553159, pow(12345678, 2, 665059573553159), None, 12345678) verify_case(0, 665059573553159, pow(1234567890, 2, 665059573553159), None, 1234567890) verify_case(0, 249712759557946793, pow(12345678, 2, 249712759557946793), None, 12345678) verify_case(0, 249712759557946793, pow(1234567890987, 2, 249712759557946793), None, 1234567890987)