# quantum_plane_factorization.sage # # SageMath implementation for the quantum-plane factorization as described in # Factoring Quantum-Plane Skew Polynomials over Q(𝝎)(𝒕) # by Mark Giesbrecht # Snapshot April 24, 2026 # # Plese don't hestiate to contact me at mwg@uwaterloo.ca if you have any # thoughts, comments or concerns! # # Implementation in Sage 10.8: https://www.sagemath.org/ # # This file contains: # 1. exact skew arithmetic over K(t), right division and right gcd; # 2. exact central-bound computation by K(T)-linear algebra; # 3. central decomposition by factoring the bound in K(T)[X]; # 4. exact irreducibility certificates from the degree-bound/norm tests; # 5. a two-level modular search for irreducible-bound pieces: # T -> a, reduction at inert finite-field primes, Sage's finite-field # skew-polynomial factorization backend, CRT/rational reconstruction # inside each specialization, rational interpolation in T, and exact # right-division certification upstairs. # # The finite-field call is through Sage's finite-field skew-polynomial class: # S = lambda_field['x', sigma_bar] # fbar.factor() # Sage documents this class as implementing the fast Caruso-Le Borgne # finite-field factorization algorithm. # # Run with: # sage quantum_plane_factorization.sage # # Or load with: # load("quantum_plane_factorization.sage") # Q = QuantumPlane(2) # x, t = Q.x(), Q.t # Q.factor((x - t)*(x - (t + 1)), verbose=True) # # IMPORTANT HONESTY NOTE # ---------------------- # The attached paper describes a certifying heuristic modular stage, not a # theorem-proved complete algorithm for every generic irreducible-bound input # over Q(omega)(t). This implementation follows that design: every returned # factor is certified by exact arithmetic in the original ring, but the modular # search has budgets and may return failure for apparently hard instances. # # LICENCE AND USE # --------------- # Copyright (c) 2026 Mark Giesbrecht # # This code is released under the MIT License. # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND. # from sage.all import * from itertools import product import sys import time # Explicit finite-field skew-ring constructor. Some Sage builds do not # dispatch lam['x', sigma] to the finite-field subclass when sigma is a # nontrivial Frobenius power, which leaves fbar as SkewPolynomial_generic_dense # without .factor(). Calling the subclass directly exposes the finite-field # factorization methods. try: from sage.rings.polynomial.skew_polynomial_ring import SkewPolynomialRing_finite_field except Exception: SkewPolynomialRing_finite_field = None # Some Sage versions return a generic quotient-ring object from k.extension(...), # even when k is finite. Sage's finite-field skew factorization backend expects # a native finite-field coefficient ring whose Frobenius morphism knows its # order and fixed field. The modular context below therefore constructs the # residue extension as GF(p^(d*m)) and chooses omega/theta inside that field. # ----------------------------------------------------------------------------- # Small utility functions # ----------------------------------------------------------------------------- def _prod(items, one): out = one for z in items: out = out * z return out def _proper_divisors(n): n = ZZ(n) return [d for d in divisors(n) if d != n] def _symmetric_residue(a, M): a = ZZ(a) % ZZ(M) M = ZZ(M) if 2*a > M: a -= M return ZZ(a) def _bounded_rational_reconstruction(a, M, numerator_bound=None, denominator_bound=None): """ Reconstruct a rational q from q == a mod M. Sage usually provides rational_reconstruction; this function uses it when available and falls back to a bounded denominator search. The fallback is enough for the small examples and is useful when users run an older Sage. """ a = ZZ(a) % ZZ(M) M = ZZ(M) # Native Sage rational reconstruction, if present in this Sage version. try: return QQ(rational_reconstruction(a, M)) except Exception: pass if denominator_bound is None: denominator_bound = max(1, floor(sqrt(M) / 2)) if numerator_bound is None: numerator_bound = max(1, floor(sqrt(M) / 2)) denominator_bound = ZZ(denominator_bound) numerator_bound = ZZ(numerator_bound) for d in range(1, denominator_bound + 1): n = _symmetric_residue(a * d, M) if abs(n) <= numerator_bound and gcd(n, d) == 1: return QQ(n) / QQ(d) raise ValueError("rational reconstruction failed modulo %s for residue %s" % (M, a)) # ----------------------------------------------------------------------------- # Core skew polynomial arithmetic with left coefficients. # ----------------------------------------------------------------------------- class SkewPoly(object): """Polynomial in x with left coefficients and multiplication x*a=sigma(a)*x.""" def __init__(self, parent, coeffs=None): self.parent = parent if coeffs is None: coeffs = [] self.coeffs = [parent.F(c) for c in coeffs] self._trim() def _trim(self): while self.coeffs and self.coeffs[-1] == 0: self.coeffs.pop() return self def copy(self): return SkewPoly(self.parent, list(self.coeffs)) def degree(self): return len(self.coeffs) - 1 def is_zero(self): return len(self.coeffs) == 0 def is_one(self): return self.degree() == 0 and self.coeffs[0] == 1 def coeff(self, i): if 0 <= i < len(self.coeffs): return self.coeffs[i] return self.parent.F(0) def lc(self): if self.is_zero(): return self.parent.F(0) return self.coeffs[-1] def _coerce(self, other): if isinstance(other, SkewPoly): if other.parent is not self.parent: raise TypeError("cannot mix skew polynomials from different parents") return other return self.parent.constant(other) def __neg__(self): return SkewPoly(self.parent, [-c for c in self.coeffs]) def __add__(self, other): other = self._coerce(other) n = max(len(self.coeffs), len(other.coeffs)) return SkewPoly(self.parent, [self.coeff(i) + other.coeff(i) for i in range(n)]) def __radd__(self, other): return self.__add__(other) def __sub__(self, other): other = self._coerce(other) n = max(len(self.coeffs), len(other.coeffs)) return SkewPoly(self.parent, [self.coeff(i) - other.coeff(i) for i in range(n)]) def __rsub__(self, other): return self.parent.constant(other).__sub__(self) def left_scalar_mul(self, c): c = self.parent.F(c) return SkewPoly(self.parent, [c * a for a in self.coeffs]) def right_scalar_mul(self, c): # (a_i x^i) * c = a_i sigma^i(c) x^i. c = self.parent.F(c) return SkewPoly(self.parent, [self.coeffs[i] * self.parent.sigma(c, i) for i in range(len(self.coeffs))]) def __mul__(self, other): if not isinstance(other, SkewPoly): return self.right_scalar_mul(other) if other.parent is not self.parent: raise TypeError("cannot mix skew polynomials from different parents") if self.is_zero() or other.is_zero(): return self.parent.zero() Q = self.parent res = [Q.F(0)] * (self.degree() + other.degree() + 1) for i, ai in enumerate(self.coeffs): if ai == 0: continue for j, bj in enumerate(other.coeffs): if bj != 0: res[i + j] += ai * Q.sigma(bj, i) return SkewPoly(Q, res) def __rmul__(self, other): if isinstance(other, SkewPoly): return other.__mul__(self) return self.left_scalar_mul(other) def __pow__(self, n): n = ZZ(n) if n < 0: raise ValueError("negative powers are not implemented in this polynomial ring") result = self.parent.one() base = self while n: if n % 2: result = result * base base = base * base n = n // 2 return result def __eq__(self, other): other = self._coerce(other) return self.coeffs == other.coeffs def __ne__(self, other): return not self.__eq__(other) def monic(self): if self.is_zero(): raise ZeroDivisionError("zero skew polynomial cannot be made monic") return self.left_scalar_mul(1 / self.lc()) def left_mul_monomial(self, c, k): """Return (c*x^k) * self.""" Q = self.parent c = Q.F(c) if self.is_zero() or c == 0: return Q.zero() out = [Q.F(0)] * (self.degree() + k + 1) for i, ai in enumerate(self.coeffs): out[i + k] = c * Q.sigma(ai, k) return SkewPoly(Q, out) def right_divmod(self, divisor): """Return q,r with self = q*divisor + r and deg(r) < deg(divisor).""" Q = self.parent divisor = self._coerce(divisor) if divisor.is_zero(): raise ZeroDivisionError("right division by zero skew polynomial") f = self.copy() g = divisor q = Q.zero() dg = g.degree() while (not f.is_zero()) and f.degree() >= dg: k = f.degree() - dg # Leading coefficient of (c*x^k)*g is c*sigma^k(lc(g)). c = f.lc() / Q.sigma(g.lc(), k) term = Q.monomial(c, k) q = q + term f = f - term * g return q, f def right_rem(self, divisor): return self.right_divmod(divisor)[1] def right_quo(self, divisor): q, r = self.right_divmod(divisor) if not r.is_zero(): raise ArithmeticError("division is not exact") return q def __repr__(self): if self.is_zero(): return "0" xname = self.parent.x_name pieces = [] one = self.parent.F(1) minus_one = self.parent.F(-1) for i, c in enumerate(self.coeffs): if c == 0: continue if i == 0: term = str(c) else: xp = xname if i == 1 else "%s^%s" % (xname, i) if c == one: term = xp elif c == minus_one: term = "-" + xp else: term = "(%s)*%s" % (c, xp) pieces.append(term) return " + ".join(pieces).replace("+ -", "- ") # ----------------------------------------------------------------------------- # Exact specialized ring L_a[x; sigma_a], with theta^m=a. # ----------------------------------------------------------------------------- class SpecializedPlane(object): """Specialized skew ring L_a[x; sigma_a], theta^m=a, sigma(theta)=omega*theta.""" def __init__(self, quantum_plane, A, x_name="x"): self.Q = quantum_plane self.m = quantum_plane.m self.K = quantum_plane.K self.omega_K = quantum_plane.omega self.A = self.K(A) self.x_name = x_name RY = PolynomialRing(self.K, "theta_poly") Y = RY.gen() self.modulus = Y**self.m - self.A if not self.modulus.is_irreducible(): raise ValueError("T -> %s is not Kummer-good: y^m-a is reducible over K" % A) self.F = self.K.extension(self.modulus, names=("theta",)) self.theta = self.F.gen() self.omega = self.F(self.omega_K) def decompose_L(self, z): """Return coefficients [c_0,...,c_{m-1}] in z=sum c_j theta^j.""" z = self.F(z) raw = list(z.list()) out = [] for j in range(self.m): out.append(self.K(raw[j]) if j < len(raw) else self.K(0)) return out def element_from_theta_coeffs(self, coeffs): z = self.F(0) for j, c in enumerate(coeffs): if c != 0: z += self.F(c) * (self.theta ** j) return self.F(z) def sigma(self, z, k=1): z = self.F(z) k = ZZ(k) % self.m if k == 0 or z == 0: return z coeffs = self.decompose_L(z) out = self.F(0) for j, c in enumerate(coeffs): if c != 0: out += self.F(c) * (self.omega ** (k*j)) * (self.theta ** j) return out def zero(self): return SkewPoly(self, []) def one(self): return SkewPoly(self, [1]) def constant(self, c): return SkewPoly(self, [self.F(c)]) def monomial(self, c, k): k = ZZ(k) coeffs = [self.F(0)] * (k + 1) coeffs[k] = self.F(c) return SkewPoly(self, coeffs) def poly(self, coeffs): return SkewPoly(self, coeffs) def x(self): return self.monomial(1, 1) def verify_product(self, factors, target): p = self.one() for h in factors: p = p * h return p == target # ----------------------------------------------------------------------------- # Main quantum-plane class over K(t). # ----------------------------------------------------------------------------- class QuantumPlane(object): """The quantum-plane skew ring K(t)[x; sigma], sigma(t)=omega*t.""" def __init__(self, m, names=("t", "x", "T", "X")): self.m = ZZ(m) if self.m < 2: raise ValueError("m must be at least 2") self.t_name, self.x_name, self.T_name, self.X_name = names if self.m == 2: # Q(zeta_2)=Q and zeta_2=-1. self.K = QQ self.omega = self.K(-1) self.K_degree = ZZ(1) else: self.K = CyclotomicField(self.m, "omega") self.omega = self.K.gen() self.K_degree = ZZ(self.K.degree()) self.Rt = PolynomialRing(self.K, self.t_name) self.t_poly = self.Rt.gen() self.F = FractionField(self.Rt) self.t = self.F(self.t_poly) self.RT = PolynomialRing(self.K, self.T_name) self.T_poly = self.RT.gen() self.KT = FractionField(self.RT) self.T = self.KT(self.T_poly) self.Rz = PolynomialRing(self.KT, "z") self.z = self.Rz.gen() self.mod_z = self.z ** self.m - self.T self.RX = PolynomialRing(self.KT, self.X_name) self.X = self.RX.gen() # ------------------------------------------------------------------ # Basic constructors and automorphisms # ------------------------------------------------------------------ def zero(self): return SkewPoly(self, []) def one(self): return SkewPoly(self, [1]) def constant(self, c): return SkewPoly(self, [self.F(c)]) def monomial(self, c, k): k = ZZ(k) if k < 0: raise ValueError("negative x-degree") coeffs = [self.F(0)] * (k + 1) coeffs[k] = self.F(c) return SkewPoly(self, coeffs) def poly(self, coeffs): return SkewPoly(self, coeffs) def x(self): return self.monomial(1, 1) def X_skew(self): return self.monomial(1, self.m) def sigma(self, a, k=1): """Apply sigma^k to a in K(t). Constants in K are fixed.""" a = self.F(a) k = ZZ(k) % self.m if k == 0 or a == 0: return a scale = self.K(self.omega ** k) num = a.numerator() den = a.denominator() return self.F(num(scale * self.t) / den(scale * self.t)) # ------------------------------------------------------------------ # K(t) <-> K(T)-basis conversion, K(t)=K(T)[z]/(z^m-T) # ------------------------------------------------------------------ def _tpoly_to_zpoly(self, p): p = self.Rt(p) out = self.Rz(0) for e, c in enumerate(p.list()): if c != 0: q, r = divmod(ZZ(e), self.m) out += self.KT(c) * (self.T ** q) * (self.z ** r) return out def decompose_F_over_KT(self, a): """Return [a_0(T),...,a_{m-1}(T)] with a=sum a_j(T)*t^j.""" a = self.F(a) if a == 0: return [self.KT(0)] * self.m P = self._tpoly_to_zpoly(a.numerator()) Q = self._tpoly_to_zpoly(a.denominator()) g, s, _ = Q.xgcd(self.mod_z) if g.degree() != 0: raise ArithmeticError("denominator was not invertible modulo z^m-T") invQ = s / self.KT(g[0]) R = (P * invQ) % self.mod_z return [self.KT(R[i]) for i in range(self.m)] def KT_to_F(self, a): """Embed a(T) in K(T) into K(t) by T -> t^m.""" a = self.KT(a) num = a.numerator() den = a.denominator() return self.F(num(self.t ** self.m) / den(self.t ** self.m)) def F_to_KT_if_fixed(self, a, check=True): comps = self.decompose_F_over_KT(a) if check: for j in range(1, self.m): if comps[j] != 0: raise ValueError("element is not fixed by sigma, so it is not in K(T)") return comps[0] def _vector_over_KT(self, p, ambient_x_degree): vec = [] for i in range(ambient_x_degree): vec.extend(self.decompose_F_over_KT(p.coeff(i))) return vec # ------------------------------------------------------------------ # Central polynomials K(T)[X] <-> skew polynomials. # ------------------------------------------------------------------ def central_to_skew(self, P): P = self.RX(P) if P == 0: return self.zero() coeffs = [self.F(0)] * (self.m * P.degree() + 1) for i in range(P.degree() + 1): coeffs[self.m * i] = self.KT_to_F(P[i]) return self.poly(coeffs) def central_factor_list(self, P): P = self.RX(P) if P.degree() <= 0: return [] return [(self.RX(f).monic(), ZZ(e)) for f, e in P.factor()] def central_is_irreducible(self, P): P = self.RX(P) if P.degree() <= 0: return False facs = self.central_factor_list(P) return len(facs) == 1 and facs[0][1] == 1 and facs[0][0].degree() == P.degree() # ------------------------------------------------------------------ # Ore right Euclidean operations. # ------------------------------------------------------------------ def gcrd(self, a, b): a = a.copy() b = b.copy() while not b.is_zero(): _, r = a.right_divmod(b) a, b = b, r return a.monic() def verify_product(self, factors, target): prod_poly = self.one() for h in factors: prod_poly = prod_poly * h return prod_poly == target # ------------------------------------------------------------------ # Bound computation. # ------------------------------------------------------------------ def bound(self, f, max_X_degree=None): """ Compute the monic central bound of monic f by K(T)-linear algebra. We find the least monic P(X) in K(T)[X] such that P(x^m) right-reduces to zero modulo f. """ if f.is_zero(): raise ZeroDivisionError("zero polynomial has no monic bound") f = f.monic() n = f.degree() if n == 0: return self.RX(1) if max_X_degree is None: max_X_degree = max(1, n * self.m) Xsk = self.X_skew() cols = [] r = self.one() rows = n * self.m for d in range(max_X_degree + 1): v = self._vector_over_KT(r, n) if d > 0: A = matrix(self.KT, rows, d) for j, col in enumerate(cols): for i, val in enumerate(col): A[i, j] = val bcol = matrix(self.KT, rows, 1, v) if A.rank() == A.augment(bcol).rank(): sol = A.solve_right(vector(self.KT, v)) coeffs = [-sol[i] for i in range(d)] + [self.KT(1)] return self.RX(coeffs).monic() cols.append(v) r = (Xsk * r).right_rem(f) raise RuntimeError("failed to find a central bound up to X-degree %s" % max_X_degree) # ------------------------------------------------------------------ # Reducible-bound central decomposition. # ------------------------------------------------------------------ def strip_right_x_power(self, f): f = f.monic() eps = 0 x = self.x() g = f while (not g.is_zero()) and g.degree() > 0 and g.coeff(0) == 0: q, r = g.right_divmod(x) if not r.is_zero(): break g = q.monic() if not q.is_zero() else q eps += 1 return g, eps def central_decomposition(self, f, verbose=False): """Return rough factors with irreducible central bounds, plus stripped x's.""" if f.is_zero(): raise ZeroDivisionError("cannot factor zero") f = f.monic() g, eps = self.strip_right_x_power(f) pieces = [] if g.is_one() else self._central_decomposition_no_x(g, verbose=verbose) pieces += [self.x()] * eps return pieces def _central_decomposition_no_x(self, f, verbose=False): P = self.bound(f) facs = self.central_factor_list(P) if verbose: print(" central bound:", P) print(" central factorization:", P.factor()) if len(facs) == 1 and facs[0][1] == 1: return [f.monic()] chosen = None for pi, _ in facs: if pi.degree() > 0: chosen = pi break if chosen is None: return [f.monic()] pi_skew = self.central_to_skew(chosen) v = self.gcrd(f, pi_skew) if v.degree() <= 0 or v.degree() >= f.degree(): for pi, _ in facs[1:]: pi_skew = self.central_to_skew(pi) v = self.gcrd(f, pi_skew) if 0 < v.degree() < f.degree(): break if v.degree() <= 0 or v.degree() >= f.degree(): raise RuntimeError("central gcd did not produce a proper factor; check the input/bound") q, r = f.right_divmod(v) if not r.is_zero(): raise ArithmeticError("internal error: central gcd factor did not divide exactly") left = self.central_decomposition(q.monic(), verbose=verbose) return left + [v.monic()] # ------------------------------------------------------------------ # Optional reduced norm sufficient irreducibility test. # ------------------------------------------------------------------ def _sigma_FX(self, elem, k, FX, XX): elem = FX(elem) if elem == 0: return FX(0) num = elem.numerator() den = elem.denominator() def sig_poly(poly): out = FX(0) for e, c in enumerate(poly.list()): if c != 0: out += FX(self.sigma(c, k)) * (FX(XX) ** e) return out return sig_poly(num) / sig_poly(den) def reduced_norm(self, p): RXF = PolynomialRing(self.F, "XX") XX = RXF.gen() FX = FractionField(RXF) b = [FX(0)] * self.m for d, ad in enumerate(p.coeffs): if ad == 0: continue q, r = divmod(ZZ(d), self.m) b[r] += FX(self.sigma(ad, -r)) * (FX(XX) ** q) M = matrix(FX, self.m, self.m) for j in range(self.m): for i in range(self.m): if b[i] == 0: continue q, r = divmod(ZZ(i + j), self.m) M[r, j] += (FX(XX) ** q) * self._sigma_FX(b[i], -j, FX, XX) det = M.det() num = det.numerator() den = det.denominator() if den.degree() != 0: raise RuntimeError("norm determinant did not simplify to a polynomial") num = num / den[0] coeffs_KT = [] for _, c in enumerate(num.list()): coeffs_KT.append(self.F_to_KT_if_fixed(c, check=True)) N = self.RX(coeffs_KT) if N == 0: return N return N.monic() # ------------------------------------------------------------------ # Characteristic-zero specialization T -> A. # ------------------------------------------------------------------ def eval_KT_at(self, r, A): r = self.KT(r) A = self.K(A) den = r.denominator()(A) if den == 0: raise ZeroDivisionError("bad specialization: denominator vanishes at T=%s" % A) return self.K(r.numerator()(A) / den) def specialize_central(self, P, A): P = self.RX(P) RXK = PolynomialRing(self.K, self.X_name) coeffs = [self.eval_KT_at(P[i], A) for i in range(P.degree() + 1)] return RXK(coeffs) def _specialization_denominators_ok(self, f, P, A): try: A = self.K(A) for c in f.coeffs: for r in self.decompose_F_over_KT(c): if r.denominator()(A) == 0: return False P = self.RX(P) for i in range(P.degree() + 1): if self.KT(P[i]).denominator()(A) == 0: return False return True except Exception: return False def is_good_specialization(self, f, P, A, verbose=False): """Check the paper's S1--S4 conditions that are directly testable here.""" if not self._specialization_denominators_ok(f, P, A): if verbose: print(" T -> %s rejected: denominator vanishes" % A) return False try: Pa = self.specialize_central(P, A) if Pa.degree() != self.RX(P).degree() or Pa == 0: return False if Pa.degree() > 0 and Pa.discriminant() == 0: return False if not Pa.is_irreducible(): if verbose: print(" T -> %s rejected: pi_a reducible: %s" % (A, Pa.factor())) return False # pi_a != X RXK = Pa.parent() if Pa == RXK.gen(): return False RY = PolynomialRing(self.K, "Y_kummer") Y = RY.gen() if not (Y**self.m - self.K(A)).is_irreducible(): if verbose: print(" T -> %s rejected: y^m-a reducible over K" % A) return False return True except Exception as err: if verbose: print(" T -> %s rejected by exception: %s" % (A, err)) return False def default_specialization_candidates(self, how_many=80): vals = [] n = ZZ(2) while len(vals) < how_many: vals.append(n) n += 1 return vals def find_good_specializations(self, f, P, count, candidates=None, verbose=False): if candidates is None: candidates = self.default_specialization_candidates(max(60, 4*count + 10)) out = [] for A in candidates: if self.is_good_specialization(f, P, A, verbose=verbose): out.append(ZZ(A) if A in ZZ else A) if len(out) >= count: return out return out def specialized_plane(self, A): return SpecializedPlane(self, A, x_name=self.x_name) def specialize_skew(self, f, A, SP=None): """Specialize f in K(t)[x;sigma] to L_A[x;sigma_A].""" if SP is None: SP = self.specialized_plane(A) coeffs = [] for c in f.coeffs: comps = self.decompose_F_over_KT(c) theta_coeffs = [self.eval_KT_at(comps[j], A) for j in range(self.m)] coeffs.append(SP.element_from_theta_coeffs(theta_coeffs)) return SP.poly(coeffs) # ------------------------------------------------------------------ # K-basis conversions for CRT/reconstruction. # ------------------------------------------------------------------ def K_element_from_QQ_coords(self, coords): if self.K == QQ: return QQ(coords[0] if coords else 0) z = self.K(0) for i, q in enumerate(coords): if q != 0: z += self.K(q) * (self.omega ** i) return self.K(z) # ------------------------------------------------------------------ # Modular stage: finite-field reductions, CRT inside specializations, # interpolation back to K(t), and exact certification. # ------------------------------------------------------------------ def _default_modular_options(self): return { "interpolation_degree_bound": 0, "specializations": None, "min_specializations": None, "primes_per_specialization": 2, "prime_search_bound": 250, "min_crt_primes": 1, "height_bound": 200, "max_modular_candidates_per_prime": 60, "max_crt_states": 400, "max_exact_specialized_candidates": 40, "enumerate_linear_field_size_bound": 5000, # Debug/performance controls for the finite-field modular backend. # If set, skip finite-field factorization when |lambda| is larger. "max_lambda_cardinality": None, # Sample additional Sage finite-field skew factorizations. This is # useful because one arbitrary noncommutative factorization may not # contain the right divisor that lifts to characteristic zero. "max_uniform_factor_trials": 0, # Bounded use of fbar.factorizations(); never leave this unbounded. "max_factorization_iterator": 0, # Bounded use of fbar.right_irreducible_divisors() for degree-e # right divisors. Keep this only for small finite fields. "max_right_irreducible_divisors": 0, "verbose_modular": False, } def factor(self, f, use_modular=True, modular_options=None, coeff_bound=1, poly_degree=1, include_ratios=False, verbose=False): """ Certified factorization. The central-decomposition stage is exact. For irreducible-bound pieces, the method tries exact irreducibility tests, then the two-level modular search, then a tiny exact linear search as a last fallback. If none of these certifies a factor or irreducibility, it raises RuntimeError. """ f = f.monic() rough = self.central_decomposition(f, verbose=verbose) out = [] for g in rough: out.extend(self._factor_piece(g.monic(), use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose)) if not self.verify_product(out, f): raise ArithmeticError("internal error: returned factors do not multiply to input") return out def _factor_piece(self, g, use_modular=True, modular_options=None, coeff_bound=1, poly_degree=1, include_ratios=False, verbose=False): g = g.monic() if g.degree() <= 0: return [] if g.degree() == 1: if verbose: print(" degree-one irreducible:", g) return [g] P = self.bound(g) if not self.central_is_irreducible(P): sub = [] for h in self.central_decomposition(g, verbose=verbose): sub.extend(self._factor_piece(h, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose)) return sub # Degree-bound test, Corollary 3.2. if P.degree() == g.degree(): if verbose: print(" irreducible by degree-bound test:", g) print(" bound =", P) return [g] # Reduced-norm sufficient test. try: N = self.reduced_norm(g) if self.central_is_irreducible(N): if verbose: print(" irreducible by reduced-norm test:", g) print(" norm =", N) return [g] except Exception as err: if verbose: print(" reduced norm test skipped:", err) # Two-level modular search. This is deliberately before the tiny exact # search so the examples visibly exercise the full modular path. if use_modular: opts = self._default_modular_options() if modular_options: opts.update(modular_options) s = ZZ(P.degree()) admissible = [] for d in divisors(self.m): e = s * ZZ(d) if 0 < e < g.degree() and g.degree() % e == 0: admissible.append(e) admissible = sorted(set(admissible)) if verbose: print(" irreducible-bound piece needs modular search:", g) print(" bound =", P) print(" admissible right-factor degrees:", admissible) for e in admissible: found = self.modular_find_right_factor(g, P, e, opts, verbose=verbose) if found is not None: h, q = found if verbose: print(" modular stage certified right factor:", h) print(" exact left quotient:", q) return (self._factor_piece(q, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose) + self._factor_piece(h, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose)) # Last-resort sparse exact search for transparent higher-degree toy examples. # This is particularly useful for m=3 examples like (x^2 - t)^2, # where the true right factor is x^2 - t and the finite-field modular # factorization may return a different local right divisor. try: s_sparse = ZZ(P.degree()) sparse_degrees = [] for d_sparse in divisors(self.m): e_sparse = s_sparse * ZZ(d_sparse) if 0 < e_sparse < g.degree() and g.degree() % e_sparse == 0: sparse_degrees.append(e_sparse) for e_sparse in sorted(set(sparse_degrees)): h, q = self.find_sparse_right_factor_degree( g, e_sparse, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose) if h is not None: if verbose: print(" found certified right factor by sparse exact search:", h) return (self._factor_piece(q, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose) + self._factor_piece(h, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose)) except Exception as err: if verbose: print(" sparse exact search skipped:", err) # Last-resort tiny exact search for transparent degree-one examples. h, q = self.find_linear_right_factor(g, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios) if h is not None: if verbose: print(" found certified right factor by tiny exact search:", h) return (self._factor_piece(q, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose) + self._factor_piece(h, use_modular=use_modular, modular_options=modular_options, coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios, verbose=verbose)) raise RuntimeError( "unresolved irreducible-bound piece: %s\n" " bound = %s\n" "The implementation is certifying: no uncertified irreducibility claim is made. " "Increase modular budgets or bounds." % (g, P)) def modular_find_right_factor(self, g, P, e, opts, verbose=False): B = ZZ(opts.get("interpolation_degree_bound", 0)) min_specs = opts.get("min_specializations", None) if min_specs is None: min_specs = max(1, 2*B + 1) spec_candidates = opts.get("specializations", None) if spec_candidates is None: spec_candidates = self.default_specialization_candidates(80) good_specs = [] for A in spec_candidates: if not self.is_good_specialization(g, P, A, verbose=opts.get("verbose_modular", False)): continue if verbose or opts.get("verbose_modular", False): print(" good central specialization T ->", A) sys.stdout.flush() exact_HA = self.specialized_right_factors_via_modular(g, P, A, e, opts, verbose=verbose or opts.get("verbose_modular", False)) if exact_HA: if verbose or opts.get("verbose_modular", False): print(" exact specialized candidates certified:", len(exact_HA)) good_specs.append((self.K(A), exact_HA)) if len(good_specs) >= min_specs: # Try all small tuples from the first min_specs successful specializations. trial_specs = good_specs[:min_specs] found = self.interpolate_specialized_candidates(g, e, trial_specs, B, opts, verbose=verbose or opts.get("verbose_modular", False)) if found is not None: return found return None def specialized_right_factors_via_modular(self, g, P, A, e, opts, verbose=False): """ Algorithm 6.3-style stage for one specialization. Returns exact monic degree-e right divisors of g_A certified in L_A[x;sigma_A]. Candidate generation comes from finite-field skew factorizations and optional exhaustive linear divisor enumeration. """ SP = self.specialized_plane(A) gA = self.specialize_skew(g, A, SP=SP) contexts = self.good_modular_prime_contexts(g, P, A, count=opts.get("primes_per_specialization", 2), prime_bound=opts.get("prime_search_bound", 250), verbose=verbose) if len(contexts) == 0: return [] per_prime = [] max_lam_card = opts.get("max_lambda_cardinality", None) for ctx in contexts: if max_lam_card is not None: try: lam_card = ZZ(ctx.lam.cardinality()) if lam_card > ZZ(max_lam_card): if verbose: print(" p = %s skipped: |lambda| = %s exceeds max_lambda_cardinality = %s" % (ctx.p, lam_card, max_lam_card)) sys.stdout.flush() continue except Exception: pass fbar = ctx.reduce_skew(g) cands = ctx.modular_right_divisor_candidates(fbar, e, max_candidates=opts.get("max_modular_candidates_per_prime", 60), enumerate_linear_field_size_bound=opts.get("enumerate_linear_field_size_bound", 5000), max_uniform_factor_trials=opts.get("max_uniform_factor_trials", 0), max_factorization_iterator=opts.get("max_factorization_iterator", 0), max_right_irreducible_divisors=opts.get("max_right_irreducible_divisors", 0), verbose=verbose) data = [ctx.residue_data_from_modular_factor(hbar, e) for hbar in cands] if verbose: print(" p = %s: modular degree-%s candidates = %s" % (ctx.p, e, len(data))) sys.stdout.flush() if not data: return [] per_prime.append((ctx, data)) if not per_prime: return [] states = [] first_ctx, first_data = per_prime[0] for dct in first_data: states.append((ZZ(first_ctx.p), dct)) exact = [] seen = set() min_crt = ZZ(opts.get("min_crt_primes", 1)) max_states = ZZ(opts.get("max_crt_states", 400)) def try_states(states, used_prime_count): if used_prime_count < min_crt: return for M, dct in states[:max_states]: try: H = self.exact_specialized_candidate_from_crt(SP, e, M, dct, opts) q, r = gA.right_divmod(H) if r.is_zero(): key = repr(H) if key not in seen: seen.add(key) exact.append(H.monic()) if verbose: print(" certified specialized factor over L_A:", H.monic()) sys.stdout.flush() if len(exact) >= opts.get("max_exact_specialized_candidates", 40): return except Exception: pass try_states(states, 1) for used_idx in range(1, len(per_prime)): ctx, data_list = per_prime[used_idx] new_states = [] for M, dct in states[:max_states]: for d2 in data_list: try: merged = {} keys = set(dct.keys()).union(set(d2.keys())) for key in keys: a1 = dct.get(key, ZZ(0)) a2 = d2.get(key, ZZ(0)) merged[key] = ZZ(CRT(ZZ(a1), ZZ(a2), ZZ(M), ZZ(ctx.p))) new_states.append((ZZ(M) * ZZ(ctx.p), merged)) if len(new_states) >= max_states: break except Exception: continue if len(new_states) >= max_states: break states = new_states try_states(states, used_idx + 1) if len(exact) >= opts.get("max_exact_specialized_candidates", 40): break return exact def good_modular_prime_contexts(self, g, P, A, count=2, prime_bound=250, verbose=False): out = [] for p in prime_range(3, prime_bound + 1): if gcd(p, self.m) != 1: continue try: ctx = ModularPrimeContext(self, A, P, p) # Force coefficient reduction once to reject primes hitting denominators. _ = ctx.reduce_skew(g) out.append(ctx) if verbose: print(" good inert modular prime p = %s, k degree = %s" % (p, ctx.K_degree)) sys.stdout.flush() if len(out) >= count: return out except Exception as err: if verbose: print(" p = %s rejected: %s" % (p, err)) sys.stdout.flush() continue return out def exact_specialized_candidate_from_crt(self, SP, e, M, residue_dict, opts): height = opts.get("height_bound", 200) coeffs = [] for ell in range(e): theta_coeffs = [] for j in range(self.m): Kcoords = [] for r in range(self.K_degree): rr = residue_dict.get((ell, j, r), ZZ(0)) q = _bounded_rational_reconstruction(rr, M, numerator_bound=height, denominator_bound=height) Kcoords.append(q) theta_coeffs.append(self.K_element_from_QQ_coords(Kcoords)) coeffs.append(SP.element_from_theta_coeffs(theta_coeffs)) coeffs.append(SP.F(1)) return SP.poly(coeffs).monic() def interpolate_specialized_candidates(self, g, e, spec_data, B, opts, verbose=False): """ Algorithm 7.2-style interpolation from exact specialized candidates. spec_data is [(A, [H_A candidates]), ...]. This searches candidate tuples, interpolates K(T)-coefficient functions of degree <= B, and certifies by exact right division in K(t)[x;sigma]. """ lists = [Hs for _, Hs in spec_data] max_tuple_search = ZZ(opts.get("max_crt_states", 400)) searched = ZZ(0) for tuple_H in product(*lists): searched += 1 if searched > max_tuple_search: break try: coeff_functions = {} for ell in range(e): for j in range(self.m): values = [] for (A, _), H in zip(spec_data, tuple_H): SP = H.parent comps = SP.decompose_L(H.coeff(ell)) values.append((self.K(A), self.K(comps[j]))) rT = self.rational_interpolate_KT(values, B) coeff_functions[(ell, j)] = rT coeffs = [] for ell in range(e): c = self.F(0) for j in range(self.m): c += self.KT_to_F(coeff_functions[(ell, j)]) * (self.t ** j) coeffs.append(c) coeffs.append(self.F(1)) Hup = self.poly(coeffs).monic() q, rem = g.right_divmod(Hup) if rem.is_zero(): if verbose: print(" interpolated and certified upstairs factor:", Hup) return Hup, q.monic() except Exception as err: if verbose: print(" interpolation tuple rejected:", err) continue return None def rational_interpolate_KT(self, values, B): """Find r(T) in K(T), deg numerator/denominator <= B, matching values.""" B = ZZ(B) if all(v == 0 for _, v in values): return self.KT(0) # Homogeneous system: p(A_i) - v_i q(A_i) = 0. ncols = 2*(B + 1) M = matrix(self.K, len(values), ncols) for row, (A, v) in enumerate(values): A = self.K(A) v = self.K(v) powA = [self.K(1)] for k in range(1, B + 1): powA.append(powA[-1] * A) for k in range(B + 1): M[row, k] = powA[k] M[row, B + 1 + k] = -v * powA[k] ns = M.right_kernel().basis() if not ns: raise ValueError("no rational interpolant of degree bound %s" % B) for vec in ns: num_coeffs = [self.K(vec[k]) for k in range(B + 1)] den_coeffs = [self.K(vec[B + 1 + k]) for k in range(B + 1)] if all(c == 0 for c in den_coeffs): continue P = self.KT(0) Q = self.KT(0) for k, c in enumerate(num_coeffs): P += self.KT(c) * (self.T ** k) for k, c in enumerate(den_coeffs): Q += self.KT(c) * (self.T ** k) if Q == 0: continue r = self.KT(P / Q) ok = True for A, v in values: try: if self.eval_KT_at(r, A) != self.K(v): ok = False break except Exception: ok = False break if ok: return r raise ValueError("no valid rational interpolant found") # ------------------------------------------------------------------ # Tiny exact search fallback. # ------------------------------------------------------------------ def small_elements(self, coeff_bound=1, poly_degree=1, include_ratios=False): vals = list(range(-ZZ(coeff_bound), ZZ(coeff_bound) + 1)) polys = [] for coeffs in product(vals, repeat=ZZ(poly_degree) + 1): p = self.Rt(0) for i, c in enumerate(coeffs): if c: p += self.K(c) * (self.t_poly ** i) fp = self.F(p) if not any(fp == q for q in polys): polys.append(fp) priority = [self.F(0), self.F(1), self.F(-1), self.t, -self.t, self.t + 1, self.t - 1, -self.t + 1, -self.t - 1] elems = [] for a in priority + polys: if not any(a == b for b in elems): elems.append(a) if include_ratios: base = list(elems) for p in base: for q in base: if q != 0: a = p / q if not any(a == b for b in elems): elems.append(a) return elems def find_linear_right_factor(self, g, coeff_bound=1, poly_degree=1, include_ratios=False): x = self.x() for a in self.small_elements(coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios): h = x - a q, r = g.right_divmod(h) if r.is_zero(): return h.monic(), q.monic() return None, None def find_sparse_right_factor_degree(self, g, e, coeff_bound=1, poly_degree=1, include_ratios=False, verbose=False): """Search for a certified sparse monic right factor x^e + a(t). This is a deliberately small exact fallback for transparent examples such as the m=3 irreducible-bound test (x^2 - t)^2, bound = X^2 - T. The two-level modular stage remains the primary path. This routine is only used after modular candidate lifting fails, and every candidate is certified by exact right division in K(t)[x;sigma]. """ e = ZZ(e) if e <= 0 or e >= g.degree(): return None, None x = self.x() tried = set() for a in self.small_elements(coeff_bound=coeff_bound, poly_degree=poly_degree, include_ratios=include_ratios): h = (x ** e + a).monic() key = repr(h) if key in tried: continue tried.add(key) try: q, r = g.right_divmod(h) if r.is_zero(): if verbose: print(" sparse exact search certified degree-%s right factor: %s" % (e, h)) return h.monic(), q.monic() except Exception: continue return None, None # ----------------------------------------------------------------------------- # One finite-field reduction attached to T -> A and rational prime p. # ----------------------------------------------------------------------------- class ModularPrimeContext(object): """ Finite-field inert-prime reduction for one specialization T -> A. This context keeps two compatible finite-field models: * k = GF(p^[K:Q]) is used for the residue field of K=Q(omega), for central irreducibility and Kummer-inertness checks. * lam = GF(p^([K:Q]*m)) is a *native Sage finite field* used as coefficient field for Sage's finite-field skew-polynomial backend. Earlier versions used ``k.extension(y^m-A)`` for lam. On some Sage builds that object is a generic quotient ring; its Frobenius endomorphism then has no ``order()``/``fixed_field()`` methods, and Sage cannot construct the finite-field skew-polynomial subclass. The native GF construction avoids that backend mismatch. """ def __init__(self, Q, A, P, p): self.Q = Q self.A = Q.K(A) self.P = Q.RX(P) self.p = ZZ(p) self.Fp = GF(self.p) self.m = Q.m self.K_degree = ZZ(Q.K_degree) self.lambda_degree = ZZ(self.K_degree * self.m) if self.K_degree == 1: self.k = GF(self.p) self.omega_bar = self.k(Q.omega) else: self.k = GF(self.p ** self.K_degree, name="wbar_%s" % self.p) self.omega_bar = self._choose_omega_bar() # Check specialized central prime modulo the residue field k. Pa = Q.specialize_central(self.P, self.A) RXk = PolynomialRing(self.k, Q.X_name) Pbar = RXk([self.reduce_K(Pa[i]) for i in range(Pa.degree() + 1)]) if Pbar.degree() != Pa.degree() or Pbar == 0: raise ValueError("central polynomial loses degree modulo p") if Pbar == RXk.gen(): raise ValueError("central polynomial reduces to X") if not Pbar.is_irreducible(): raise ValueError("central polynomial reducible modulo p") self.Pbar = Pbar # Check inertness of the Kummer extension over k. RY = PolynomialRing(self.k, "theta_bar_poly") Y = RY.gen() self.Abar = self.reduce_K(self.A) self.kummer_mod = Y**self.m - self.Abar if not self.kummer_mod.is_irreducible(): raise ValueError("not inert: y^m-A reducible over residue field") # Native finite field for Sage's finite-field skew backend. self.lam = GF(self.p ** self.lambda_degree, name="lambda_bar_%s" % self.p) self._choose_native_omega_theta_and_sigma_power() self.sigma = self.lam.frobenius_endomorphism(self.sigma_power) # Sanity check: this must be a finite-field Frobenius, not the generic # prime-characteristic Frobenius. Sage documents order()/fixed_field() # on finite-field Frobenius maps; the skew ring constructor relies on it. if not hasattr(self.sigma, "order"): raise RuntimeError( "native GF construction still produced a Frobenius without " "order(); coefficient field class is %s" % type(self.lam)) if ZZ(self.sigma.order()) != self.m: raise RuntimeError("finite-field sigma has order %s, expected %s" % (self.sigma.order(), self.m)) self._prepare_basis_matrices() self._prepare_lambda_theta_omega_matrix() # Force Sage's finite-field skew-polynomial subclass. The bracket # constructor is correct on many Sage versions, but on some builds it # dispatches to the generic dense class for Frobenius powers. The # generic class has no finite-field factorization method. if SkewPolynomialRing_finite_field is None: self.S = self.lam['x', self.sigma] else: try: self.S = SkewPolynomialRing_finite_field( self.lam, self.sigma, None, names=("x",), sparse=False) except TypeError: # Older Sage constructor signature uses positional names/sparse. self.S = SkewPolynomialRing_finite_field( self.lam, self.sigma, None, ("x",), False) self.x = self.S.gen() probe = self.S(1) * self.x if not (hasattr(probe, "factor") or hasattr(probe, "factorizations") or hasattr(probe, "right_irreducible_divisor")): raise RuntimeError( "Sage did not construct the finite-field skew-polynomial " "backend; got element class %s" % type(probe)) def _ff_prime_vector(self, field, elem): """Return prime-field coordinates of elem in a finite field. Sage's small finite fields are often backed by Givaro/Zech logs. In those fields elements may not implement ``.vector()``. Use the public vector-space map first; fall back to ``.list()``/polynomial coordinates in the field generator basis. """ elem = field(elem) deg = ZZ(field.degree()) if field == self.Fp or deg == 1: return vector(self.Fp, [self.Fp(elem)]) try: V, from_V, to_V = field.vector_space(self.Fp, map=True) v = to_V(elem) return vector(self.Fp, [self.Fp(v[i]) for i in range(deg)]) except Exception: pass try: V, from_V, to_V = field.vector_space(map=True) v = to_V(elem) return vector(self.Fp, [self.Fp(v[i]) for i in range(deg)]) except Exception: pass try: coeffs = list(elem.list()) except Exception: try: coeffs = list(elem.polynomial().list()) except Exception as err: raise ValueError( "cannot get prime-field coordinates in GF(%s^%s): %s" % (self.p, field.degree(), err)) coeffs = [self.Fp(c) for c in coeffs] if len(coeffs) < deg: coeffs += [self.Fp(0)] * (deg - len(coeffs)) return vector(self.Fp, coeffs[:deg]) def _choose_omega_bar(self): q = self.p ** self.K_degree gen = self.k.multiplicative_generator() candidates = [] for r in range(1, self.m + 1): if gcd(r, self.m) == 1: candidates.append(gen ** (((q - 1) // self.m) * r)) for cand in candidates: if cand ** self.m != 1: continue if any(cand ** d == 1 for d in _proper_divisors(self.m)): continue cols = [self._ff_prime_vector(self.k, cand ** i) for i in range(self.K_degree)] B = matrix(self.Fp, self.K_degree, self.K_degree) for col, v in enumerate(cols): for row in range(self.K_degree): B[row, col] = v[row] if B.rank() == self.K_degree: return cand raise ValueError("could not choose omega_bar of full degree over F_p") def _choose_native_omega_candidates(self): """Primitive m-th roots in lam that generate the cyclotomic residue field.""" qlam = self.p ** self.lambda_degree gen = self.lam.multiplicative_generator() out = [] if (qlam - 1) % self.m != 0: return out for r in range(1, self.m + 1): if gcd(r, self.m) != 1: continue cand = gen ** (((qlam - 1) // self.m) * r) if cand ** self.m != 1: continue if any(cand ** d == 1 for d in _proper_divisors(self.m)): continue # It must lie in the degree-[K:Q] fixed field of Frobenius^K_degree. if cand ** (self.p ** self.K_degree) != cand: continue cols = [self._ff_prime_vector(self.lam, cand ** i) for i in range(self.K_degree)] B = matrix(self.Fp, self.lambda_degree, self.K_degree) for col, v in enumerate(cols): for row in range(self.lambda_degree): B[row, col] = v[row] if B.rank() == self.K_degree: out.append(cand) return out def _reduce_K_to_lam_with_omega(self, c, omega_lam): K = self.Q.K if K == QQ: return self.lam(self._reduce_QQ(c)) c = K(c) poly = c.polynomial() out = self.lam(0) for i, qi in enumerate(poly.list()): if qi != 0: out += self.lam(self._reduce_QQ(qi)) * (omega_lam ** i) return out def _choose_native_omega_theta_and_sigma_power(self): """Choose omega, theta in native lam with theta^m=A and sigma(theta)=omega*theta.""" RYlam = PolynomialRing(self.lam, "theta_native_poly") Ylam = RYlam.gen() last_reason = None for omega_lam in self._choose_native_omega_candidates(): Abar_lam = self._reduce_K_to_lam_with_omega(self.A, omega_lam) roots = (Ylam**self.m - Abar_lam).roots(multiplicities=False) if not roots: last_reason = "y^m-A has no roots in native lambda field" continue for theta in roots: if theta == 0: continue for r in range(self.m): power = self.K_degree * r if theta ** (self.p ** power) == omega_lam * theta: self.omega_lam = omega_lam self.Abar_lam = Abar_lam self.theta_bar = theta self.sigma_power = ZZ(power) if self.sigma_power == 0 and self.m > 1: continue return last_reason = "no Frobenius power sends theta to omega*theta" raise ValueError("could not choose native omega/theta in GF(p^(d*m)): %s" % last_reason) def _prepare_basis_matrices(self): # Matrix converting k elements to coordinates in 1,omega,...,omega^{d-1}. self.omega_basis_matrix = matrix(self.Fp, self.K_degree, self.K_degree) for col in range(self.K_degree): v = self._ff_prime_vector(self.k, self.omega_bar ** col) for row in range(self.K_degree): self.omega_basis_matrix[row, col] = v[row] if self.omega_basis_matrix.rank() != self.K_degree: raise ValueError("omega powers are not a basis") def _prepare_lambda_theta_omega_matrix(self): """Matrix for lam coordinates in theta^j omega^r over F_p.""" N = self.lambda_degree self.lambda_basis_pairs = [] self.lambda_basis_matrix = matrix(self.Fp, N, N) col = 0 for j in range(self.m): for r in range(self.K_degree): elem = (self.theta_bar ** j) * (self.omega_lam ** r) v = self._ff_prime_vector(self.lam, elem) for row in range(N): self.lambda_basis_matrix[row, col] = v[row] self.lambda_basis_pairs.append((j, r)) col += 1 if self.lambda_basis_matrix.rank() != N: raise ValueError("theta/omega powers are not a prime-field basis of lambda") def _find_sigma_frobenius_power(self): # Kept for compatibility with older fragments; native construction sets # self.sigma_power directly. target = self.omega_lam * self.theta_bar for r in range(self.m): power = self.K_degree * r if self.theta_bar ** (self.p ** power) == target: return ZZ(power) raise ValueError("could not realize theta -> omega*theta as a Frobenius power") def _reduce_QQ(self, q): q = QQ(q) den = ZZ(q.denominator()) if den % self.p == 0: raise ZeroDivisionError("rational denominator divisible by p") return self.Fp(ZZ(q.numerator())) / self.Fp(den) def reduce_K(self, c): K = self.Q.K if K == QQ: return self.k(self._reduce_QQ(c)) c = K(c) poly = c.polynomial() out = self.k(0) for i, qi in enumerate(poly.list()): if qi != 0: out += self.k(self._reduce_QQ(qi)) * (self.omega_bar ** i) return out def reduce_K_lam(self, c): return self._reduce_K_to_lam_with_omega(c, self.omega_lam) def reduce_KT_at_A(self, r): r = self.Q.KT(r) num = r.numerator()(self.A) den = r.denominator()(self.A) den_red = self.reduce_K_lam(den) if den_red == 0: raise ZeroDivisionError("T-denominator becomes zero modulo p") return self.reduce_K_lam(num) / den_red def reduce_F_specialized(self, c): comps = self.Q.decompose_F_over_KT(c) out = self.lam(0) for j in range(self.m): kj = self.reduce_KT_at_A(comps[j]) if kj != 0: out += self.lam(kj) * (self.theta_bar ** j) return out def reduce_skew(self, f): coeffs = [self.reduce_F_specialized(f.coeff(i)) for i in range(f.degree() + 1)] out = self.S(0) for i, c in enumerate(coeffs): if c != 0: out += self.S(c) * (self.x ** i) return out def _sage_coeff(self, fbar, i): try: return self.lam(fbar[i]) except Exception: try: coeffs = fbar.list() return self.lam(coeffs[i]) if i < len(coeffs) else self.lam(0) except Exception: return self.lam(0) def _right_divides_finite(self, fbar, hbar): try: return (fbar % hbar) == 0 except Exception: pass for name in ["is_right_divisible_by"]: if hasattr(fbar, name): try: return bool(getattr(fbar, name)(hbar)) except Exception: pass for name in ["right_quo_rem", "quo_rem"]: if hasattr(fbar, name): try: q, r = getattr(fbar, name)(hbar) return r == 0 except Exception: pass return False def _atoms_from_factorization(self, Fac): """Convert a Sage Factorization object to an ordered list of monic atoms.""" atoms = [] for fac, mult in Fac: for _ in range(ZZ(mult)): atoms.append(fac) return atoms def _finite_field_factor_atoms(self, fbar, verbose=False, uniform=False): """Return one ordered list of monic finite-field skew factors for fbar.""" last_err = None if hasattr(fbar, "factor"): try: if verbose: print(" starting Sage finite-field skew factor(uniform=%s), degree=%s, |lambda|=%s" % (uniform, fbar.degree(), self.lam.cardinality())) sys.stdout.flush() t0 = time.time() if uniform: Fac = fbar.factor(uniform=True) else: Fac = fbar.factor() atoms = self._atoms_from_factorization(Fac) if verbose: print(" finished factor() in %.3fs; atom degrees = %s" % (time.time() - t0, [h.degree() for h in atoms])) sys.stdout.flush() return atoms except Exception as err: last_err = err if verbose: print(" factor() failed/skipped:", err) sys.stdout.flush() if hasattr(fbar, "factorizations"): try: if verbose: print(" starting one step of fbar.factorizations()") sys.stdout.flush() t0 = time.time() Fac = next(fbar.factorizations()) atoms = self._atoms_from_factorization(Fac) if verbose: print(" finished factorizations() step in %.3fs; atom degrees = %s" % (time.time() - t0, [h.degree() for h in atoms])) sys.stdout.flush() return atoms except Exception as err: last_err = err if verbose: print(" Sage finite-field skew factorization unavailable/skipped:", last_err) sys.stdout.flush() return [] def modular_right_divisor_candidates(self, fbar, e, max_candidates=60, enumerate_linear_field_size_bound=5000, max_uniform_factor_trials=0, max_factorization_iterator=0, max_right_irreducible_divisors=0, verbose=False): candidates = [] seen = set() def add_candidate(h): if h.degree() != e: return try: if not self._right_divides_finite(fbar, h): return except Exception: return key = repr(h) if key not in seen: seen.add(key) candidates.append(h) if verbose: print(" kept modular right divisor candidate #%s of degree %s" % (len(candidates), e)) sys.stdout.flush() def add_from_atoms(atoms): if not atoms: return # Suffixes of this noncommutative factorization are right divisors. for start in range(len(atoms)): h = self.S(1) for fac in atoms[start:]: h = h * fac add_candidate(h) if len(candidates) >= max_candidates: return # Also add contiguous products of exactly the requested degree. for start in range(len(atoms)): h = self.S(1) for stop in range(start, len(atoms)): h = h * atoms[stop] if h.degree() == e: add_candidate(h) if h.degree() >= e or len(candidates) >= max_candidates: break if len(candidates) >= max_candidates: return # Sage finite-field skew factorization backend: default cached/random factorization. add_from_atoms(self._finite_field_factor_atoms(fbar, verbose=verbose, uniform=False)) if len(candidates) >= max_candidates: return candidates # Additional independent random/uniform factorizations. Sage documents # factor(uniform=True) for finite-field skew polynomials; this can expose # different right divisors when the first noncommutative factorization # does not lift to characteristic zero. for trial in range(ZZ(max_uniform_factor_trials)): add_from_atoms(self._finite_field_factor_atoms(fbar, verbose=verbose, uniform=True)) if len(candidates) >= max_candidates: return candidates # Bounded iterator over factorizations. This can be enormous; use only # with a small positive bound. if ZZ(max_factorization_iterator) > 0 and hasattr(fbar, "factorizations"): try: it = fbar.factorizations() for idx in range(ZZ(max_factorization_iterator)): if verbose: print(" requesting factorization iterator item", idx + 1) sys.stdout.flush() Fac = next(it) add_from_atoms(self._atoms_from_factorization(Fac)) if len(candidates) >= max_candidates: return candidates except StopIteration: pass except Exception as err: if verbose: print(" bounded factorization iterator stopped:", err) sys.stdout.flush() # Bounded enumeration of irreducible right divisors. For the m=3, # e=2, pi=X^2-T toy example, this is the most direct way to find the # degree-2 divisor x^2-theta in the small p=5 residue field. if ZZ(max_right_irreducible_divisors) > 0 and hasattr(fbar, "right_irreducible_divisors"): try: if verbose: print(" enumerating up to %s irreducible right divisors" % max_right_irreducible_divisors) sys.stdout.flush() it = fbar.right_irreducible_divisors() for idx in range(ZZ(max_right_irreducible_divisors)): h = next(it) if h.degree() == e: add_candidate(h) if len(candidates) >= max_candidates: return candidates except StopIteration: pass except Exception as err: if verbose: print(" right-divisor enumeration stopped:", err) sys.stdout.flush() # Exhaustive linear right-divisor enumeration. This is intentionally # bounded, and it makes small irreducible-bound examples robust when the # finite-field factorization chooses a different noncommutative ordering. if e == 1 and self.lam.cardinality() <= enumerate_linear_field_size_bound: for c in self.lam: h = self.x + self.S(c) add_candidate(h) if len(candidates) >= max_candidates: break return candidates def lambda_theta_omega_coords(self, z): """Return coordinates of z in theta^j omega^r over F_p.""" z = self.lam(z) v = self._ff_prime_vector(self.lam, z) sol = self.lambda_basis_matrix.solve_right(v) data = {} for idx, (j, r) in enumerate(self.lambda_basis_pairs): data[(j, r)] = ZZ(sol[idx]) % self.p return data def lambda_to_theta_coeffs(self, z): """Compatibility helper returning theta-coefficients as omega-coordinate lists.""" coords = self.lambda_theta_omega_coords(z) return [[coords[(j, r)] for r in range(self.K_degree)] for j in range(self.m)] def k_to_omega_coords(self, b): b = self.k(b) if self.K_degree == 1: return [ZZ(self.Fp(b))] v = self._ff_prime_vector(self.k, b) sol = self.omega_basis_matrix.solve_right(v) return [ZZ(sol[i]) for i in range(self.K_degree)] def residue_data_from_modular_factor(self, hbar, e): """Residues of lower coefficients in theta^j omega^r basis.""" data = {} for ell in range(e): z = self._sage_coeff(hbar, ell) coords = self.lambda_theta_omega_coords(z) for j in range(self.m): for r in range(self.K_degree): data[(ell, j, r)] = coords[(j, r)] return data # ----------------------------------------------------------------------------- # Example runner # ----------------------------------------------------------------------------- def show_factorization(Q, title, f, modular_options=None, use_modular=True): print("\n" + "=" * 78) print(title) print("m =", Q.m, " omega =", Q.omega) print("f =", f) B = Q.bound(f) print("bound(f) =", B) print("factor(bound) =", B.factor()) rough = Q.central_decomposition(f) print("rough central-decomposition pieces:") for i, g in enumerate(rough, 1): print(" g%s = %s" % (i, g)) print(" bound =", Q.bound(g)) factors = Q.factor(f, use_modular=use_modular, modular_options=modular_options, verbose=True, coeff_bound=2, poly_degree=3, include_ratios=False) print("certified factorization:") for i, h in enumerate(factors, 1): print(" h%s = %s" % (i, h)) print("verified product equals f:", Q.verify_product(factors, f)) def run_examples(): # ------------------------------------------------------------------ # m=2: K=Q, sigma(t)=-t, centre Q(T)[X]. # ------------------------------------------------------------------ Q = QuantumPlane(2) x, t = Q.x(), Q.t T_expr = t**2 # Example 1: reducible bound. This is solved entirely by the central stage. f1 = (x - t) * (x - (t + 1)) show_factorization(Q, "Example 1: reducible bound, central decomposition splits two linear factors", f1, modular_options={"interpolation_degree_bound": 0, "specializations": [2, 3, 5]}) # Example 2: irreducible central bound, but reducible skew polynomial. # X+T = (x-t)^2 for sigma(t)=-t. This exercises the modular stage. f2 = x**2 + t**2 show_factorization(Q, "Example 2: irreducible bound but reducible; modular stage recovers x-t", f2, modular_options={ "interpolation_degree_bound": 0, "specializations": [2, 3, 5, 6, 7], "primes_per_specialization": 2, "min_crt_primes": 1, "verbose_modular": True, }) # Example 3: irreducible bound, reducible, and the right factor depends on T. # Let a(t)=t*(T+1)=t*(t^2+1). Since sigma(a)=-a, # (x-a)^2 = X + a^2 = X + T*(T+1)^2, # so the bound is irreducible linear in X, while the skew polynomial is reducible. # Recovering a(t) requires interpolation in T with degree bound B=1. a = t * (T_expr + 1) f3 = (x - a) * (x - a) show_factorization(Q, "Example 3: full two-level modular lift with interpolation in T, factor x-t*(T+1)", f3, modular_options={ "interpolation_degree_bound": 1, "specializations": [2, 3, 5, 6, 7, 10, 11], "min_specializations": 3, "primes_per_specialization": 2, "min_crt_primes": 1, "height_bound": 100, "verbose_modular": True, "max_crt_states": 2000, "max_exact_specialized_candidates": 80, }) # Example 4: irreducible bound and irreducible skew polynomial by degree test. q = x**2 + Q.monomial(t, 1) + 1 show_factorization(Q, "Example 4: irreducible bound and irreducible by degree-bound test", q, modular_options={"interpolation_degree_bound": 0}) # Example 5: reducible bound with an irreducible-bound piece. f5 = q * (x - (t + 1)) show_factorization(Q, "Example 5: reducible bound; central stage splits off a linear and an irreducible piece", f5, modular_options={"interpolation_degree_bound": 0, "specializations": [2, 3, 5]}) # Optional m=3 smoke test. This shows the code can construct the finite-field # K=Q(zeta_3) modular backend. It is wrapped in try/except because available # Sage versions differ most around relative finite-field displays/coercions. print("\n" + "=" * 78) print("Optional m=3 smoke test over Q(zeta_3)") try: Q3 = QuantumPlane(3) x3, t3 = Q3.x(), Q3.t f6 = (x3 - t3) * (x3 - (t3 + 1)) print("f =", f6) print("bound(f) =", Q3.bound(f6)) print("central pieces:", Q3.central_decomposition(f6)) # Build one finite-field modular context for the first rough piece. # IMPORTANT: for m=3 and K=Q(zeta_3), rational integer values A make # y^3-A reducible at every inert rational prime p: A lies in F_p^*, and # F_p^* is contained in the cube subgroup of F_{p^2}^*. Use genuinely # K-valued specializations A=n+omega instead. g = Q3.central_decomposition(f6)[-1] P = Q3.bound(g) spec_candidates = [Q3.K(n) + Q3.omega for n in range(2, 80)] specs = Q3.find_good_specializations(g, P, 8, candidates=spec_candidates, verbose=True) if not specs: raise RuntimeError("m=3 smoke found no good non-rational T-specialization") ctxs = [] chosen_A = None for A in specs: trial = Q3.good_modular_prime_contexts(g, P, A, count=1, prime_bound=5000, verbose=True) if trial: ctxs = trial chosen_A = A break if not ctxs: raise RuntimeError("m=3 smoke found no good inert finite-field prime") print("chosen specialization T ->", chosen_A) fbar = ctxs[0].reduce_skew(g) print("finite-field skew parent:", fbar.parent()) print("finite-field skew element class:", type(fbar)) atoms = ctxs[0]._finite_field_factor_atoms(fbar, verbose=True) if not atoms: raise RuntimeError("finite-field skew factorization backend exposed no factors") print("one reduced finite-field skew factorization atom list:", atoms) assert ctxs[0].modular_right_divisor_candidates( fbar, g.degree(), max_candidates=10, enumerate_linear_field_size_bound=5000, verbose=True) print("m=3 smoke test passed") except Exception as err: print("m=3 smoke test skipped:", err) if __name__ == "__main__": run_examples()