cubic.py
"""cubic.py -- complex cubic fields read through binary cubic forms: the
ring of integers, the prime places, a class group checked against the
analytic class number formula, and the class of a degree-1 place.
USE. In Python, with mpmath installed, from the folder holding
cubic.py:
>>> from cubic import enumerate_fields, Field
>>> fields = enumerate_fields(300) # every one with |d_K| <= 300
>>> len(fields)
32
>>> form = next(f for f, d in fields.items() if d == -283)
>>> K = Field(form, -283)
>>> K.class_group(), K.h # checked against the formula
(True, 2)
>>> [K.place_class(P) for P in K.places(2) if P.deg == 1]
[(1,)]
WHAT IT IS FOR. A stranger holding a complex cubic field, by a
discriminant bound (enumerate_fields) or by an irreducible defining
polynomial g(x) = a x^3 - b x^2 + c x - d of negative discriminant
(the form (a, b, c, d) passed through maximal_form, which refuses
|disc| past about 6.4 * 10^13), gets its maximal order, its places
over any prime, its class group with that check, and the class of a
place of degree 1 (None when no element within place_class's search
radius places it), in pure Python. Field takes a maximal form: it
refuses one of discriminant >= 0, or one whose discriminant is not
the d_K passed beside it. triple.py reads the split primes' triples
of classes with it.
THE ROUTE.
(1) FIELDS (Hunter's theorem, Cohen, A Course in Computational
Algebraic Number Theory, Thm 6.4.2). Every cubic field is Q(theta)
with theta integral, Tr theta in {0, 1} and T2(theta) <= Tr^2/3 +
(2/3) sqrt|d_K|. Each such minimal polynomial x^3 - s1 x^2 + s2 x
- s3 is a binary cubic form, and its ring is Z[theta].
(2) THE MAXIMAL ORDER (Delone and Faddeev; Davenport and Heilbronn).
A binary cubic form f = (a, b, c, d), f(x, y) = a x^3 + b x^2 y +
c x y^2 + d y^3, is the ring with basis 1, w, v and
w v = -ad, w^2 = -ac + b w - a v, v^2 = -bd + d w - c v,
of discriminant disc(f); GL2(Z) acting by f -> det(g)^-1 f((x,y)g)
is ring isomorphism. The ring is not maximal at p exactly when f =
0 mod p or f is equivalent to a form with p^2 | a and p | b, and
then (a/p^2, b/p, c, p d) is the form of an overorder of index p.
A form maximal at every p is O_K, and p factors in O_K as f
factors mod p. With xi a root of g(x) = f(x, -1) = a x^3 - b x^2 +
c x - d, w = a xi and v = -(a xi^2 - b xi + c).
(3) ONE KEY PER FIELD. A complex cubic form has one real and two
complex roots; the quadratic factor through the complex pair is a
positive definite covariant, Gauss-reduced, and the key is the
least form (sign fixed by a > 0) over the substitutions keeping it
reduced.
(4) THE CLASS GROUP. The prime places of norm up to B, B at least the
Minkowski bound (4/pi)(2/9) sqrt|d_K|, generate Cl. Relations come
from small elements (an LLL-reduced basis under T2) whose norms are
B-smooth, and from the rational primes (p); each valuation is EXACT:
a degree-1 place unramified over p is a p-adic root of g, lifted by
Newton, and v_P(x + y w + z v) is v_p of the image; the other
places over p take what the norm leaves. The quotient of Z^B by the
relations maps ONTO Cl, and a relation reducing to zero is a unit,
whose log at the real place is carried through the reduction. So
h_found R_found = h R m with m a positive integer, and the analytic
class number formula hR = sqrt|d_K| L / (2 pi), L the Euler product
of zeta_K / zeta to 3 * 10^4, puts m = 1 when their ratio is below 2
by a margin. The truncated product carries no proved error bound,
so the check is numerical evidence for m = 1, not a proof. At
h_found = 1 the map onto Cl already gives h = 1; always_cert
asks for the ratio there too, as a control of the whole
route. The logs of the units are carried to 60 digits through a
Hermite form kept reduced, and read first off relations sharing an
exponent row, so they stay small; a regulator below 0.28 (d_K =
-23's, the least) is refused wherever the formula is read:
class_group returns False with floor_bad set.
(5) THE CLASS OF A PLACE ABOVE B: an element of its own ideal lattice
(LLL-reduced) whose norm is p times a B-smooth number places it.
"""
import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import cmath
from math import gcd, isqrt, log, pi, sqrt
import mpmath
mpmath.mp.dps = 60
MPF = mpmath.mpf
def primes_below(n):
s = bytearray([1]) * n
s[0:2] = b"\x00\x00"
for i in range(2, isqrt(n - 1) + 1):
if s[i]:
s[i * i::i] = bytearray(len(range(i * i, n, i)))
return [i for i in range(n) if s[i]]
PRIMES = primes_below(40000)
UNIT_TOL = MPF("1e-15") # a remainder below this is zero
R_FLOOR = 0.28 # d_K = -23's regulator, the least (see _certify)
# ---------------------------------------------------------------------------
# binary cubic forms
def disc(f):
a, b, c, d = f
return (b * b * c * c - 4 * a * c ** 3 - 4 * b ** 3 * d
- 27 * a * a * d * d + 18 * a * b * c * d)
def _mul(P, Q):
R = [0] * (len(P) + len(Q) - 1)
for i, u in enumerate(P):
for j, v in enumerate(Q):
R[i + j] += u * v
return R
def subst(f, al, be, ga, de):
"""det^-1 f(al X + be Y, ga X + de Y), coefficients by power of Y."""
a, b, c, d = f
x, y = [al, be], [ga, de]
x2, y2 = _mul(x, x), _mul(y, y)
out = [0, 0, 0, 0]
for k, P in ((a, _mul(x2, x)), (b, _mul(x2, y)), (c, _mul(x, y2)),
(d, _mul(y2, y))):
for i in range(4):
out[i] += k * P[i]
det = al * de - be * ga
assert det in (1, -1)
return tuple(det * o for o in out)
def maximalize(f, p):
"""The form of the p-maximal overorder of f's ring."""
while True:
a, b, c, d = f
if a % p == 0 and b % p == 0 and c % p == 0 and d % p == 0:
f = (a // p, b // p, c // p, d // p)
continue
if a % p == 0 and b % p == 0:
F = f
else:
for r in range(p):
if (a * r * r * r + b * r * r + c * r + d) % p == 0 and \
(3 * a * r * r + 2 * b * r + c) % p == 0:
F = subst(f, r, 1, 1, 0)
break
else:
return f
A, Bc, C, D = F
if A % (p * p) == 0 and Bc % p == 0:
f = (A // (p * p), Bc // p, C, D * p)
continue
return F
def squarefull_primes(n):
"""The primes p with p^2 | n; n = 0, and |n| past the trial
division's range, refused."""
n = abs(n)
if n == 0:
raise ValueError("every prime squared divides 0")
out = []
lim = round(n ** (1 / 3)) + 2
if lim > PRIMES[-1]:
raise ValueError("|n| past the trial division's range")
for p in PRIMES:
if p > lim:
break
if n % p == 0:
k = 0
while n % p == 0:
n //= p
k += 1
if k >= 2:
out.append(p)
if n > 1:
q = isqrt(n)
if q * q == n:
out.append(q)
return out
def maximal_form(f):
for p in squarefull_primes(disc(f)):
f = maximalize(f, p)
return f
def _roots_complex(f):
"""The real root and the root of positive imaginary part of f(t, 1)."""
a, b, c, d = f
rs = _cubic_roots(a, b, c, d)
real = min(rs, key=lambda z: abs(z.imag)).real
cx = max(rs, key=lambda z: z.imag)
return real, cx
def _cubic_roots(a, b, c, d):
"""All three roots of a t^3 + b t^2 + c t + d, polished by Newton."""
B, C, D = b / a, c / a, d / a
p = C - B * B / 3
q = 2 * B ** 3 / 27 - B * C / 3 + D
disc_ = q * q / 4 + p ** 3 / 27
s = cmath.sqrt(disc_)
u = (-q / 2 + s)
u = u ** (1 / 3) if u != 0 else 0
if abs(u) < 1e-300:
u = (-q / 2 - s) ** (1 / 3)
w = complex(-0.5, sqrt(3) / 2)
roots = []
for k in range(3):
uk = u * w ** k
t = uk - p / (3 * uk) if uk != 0 else 0
roots.append(t - B / 3)
out = []
for z in roots:
for _ in range(3):
fz = ((a * z + b) * z + c) * z + d
dz = (3 * a * z + 2 * b) * z + c
if dz == 0:
break
z = z - fz / dz
out.append(z)
return out
_SMALL = [(al, be, ga, de) for al in (-1, 0, 1) for be in (-1, 0, 1)
for ga in (-1, 0, 1) for de in (-1, 0, 1)
if al * de - be * ga in (1, -1)]
def canonical(f):
"""One key per GL2(Z)-class of a complex cubic form (see (3)). The
reduction runs in floats with tolerance 1e-9 A, so the key's
uniqueness rests on them; Hasse's count at f = 1 (conductor.py C3)
and the independent field counts (triple.py K1) control it."""
rho, sig = _roots_complex(f)
# J(x, y) = (x - sig y)(x - conj(sig) y) = A x^2 + B xy + C y^2
A, B, C = 1.0, -2 * sig.real, abs(sig) ** 2
for _ in range(200):
k = -int(round(B / (2 * A)))
if k:
f = subst(f, 1, k, 0, 1)
A, B, C = A, 2 * k * A + B, A * k * k + B * k + C
if C < A:
f = subst(f, 0, -1, 1, 0)
A, B, C = C, -B, A
continue
break
best = None
tol = 1e-9 * A
for (al, be, ga, de) in _SMALL:
A2 = A * al * al + B * al * ga + C * ga * ga
B2 = 2 * A * al * be + B * (al * de + be * ga) + 2 * C * ga * de
C2 = A * be * be + B * be * de + C * de * de
if abs(B2) <= A2 + tol and A2 <= C2 + tol:
g = subst(f, al, be, ga, de)
if g[0] < 0:
g = tuple(-x for x in g)
if best is None or g < best:
best = g
return best
def _irreducible_monic(s1, s2, s3):
"""x^3 - s1 x^2 + s2 x - s3 has no integer root (s3 != 0)."""
n = abs(s3)
r = 1
while r * r <= n:
if n % r == 0:
for t in (r, -r, n // r, -(n // r)):
if ((t - s1) * t + s2) * t - s3 == 0:
return False
r += 1
return True
def enumerate_fields(cap):
"""Every complex cubic field with |d_K| <= cap: {key form: d_K}."""
fields = {}
for s1 in (0, 1):
t2max = s1 * s1 / 3 + (2 / 3) * sqrt(cap) + 1e-9
s2lo = -int((t2max - s1 * s1) // 2) - 1
s2hi = int((s1 * s1 + t2max) // 2) + 1
s3max = int((t2max / 3) ** 1.5) + 1
for s2 in range(s2lo, s2hi + 1):
for s3 in range(0 if s1 == 0 else -s3max, s3max + 1):
if s3 == 0:
continue
A, Bq, Cq = -s1, s2, -s3
D = (A * A * Bq * Bq - 4 * Bq ** 3 - 4 * A ** 3 * Cq
- 27 * Cq * Cq + 18 * A * Bq * Cq)
if D >= 0:
continue
# the real root rho; the complex pair has |z|^2 = s3/rho
rs = _cubic_roots(1, -s1, s2, -s3)
rho = min(rs, key=lambda z: abs(z.imag)).real
t2 = rho * rho + 2 * s3 / rho
if t2 > t2max:
continue
if not _irreducible_monic(s1, s2, s3):
continue
# g(x) = f(x, -1) = x^3 - s1 x^2 + s2 x - s3: a, b, c, d =
# 1, s1, s2, s3
f = maximal_form((1, s1, s2, s3))
dk = disc(f)
if -dk > cap:
continue
key = canonical(f)
if key not in fields:
fields[key] = dk
return fields
# ---------------------------------------------------------------------------
# the field
def _vp(n, p):
if n == 0:
return 10 ** 9
k = 0
while n % p == 0:
n //= p
k += 1
return k
def xgcd(a, b):
x0, x1, y0, y1 = 1, 0, 0, 1
while b:
q, a, b = a // b, b, a % b
x0, x1 = x1, x0 - q * x1
y0, y1 = y1, y0 - q * y1
return a, x0, y0
def _mult(coeffs, r, p):
"""The multiplicity of r as a root mod p (coefficients from the top)."""
m = 0
cs = [x % p for x in coeffs]
while len(cs) > 1:
q, acc = [], 0
for x in cs:
acc = (acc * r + x) % p
q.append(acc)
if q[-1]:
break
m += 1
cs = q[:-1]
return m
def legendre(a, p):
a %= p
if a == 0:
return 0
return 1 if pow(a, (p - 1) // 2, p) == 1 else -1
class Place:
__slots__ = ("p", "deg", "e", "kind", "r", "idx")
def __init__(self, p, deg, e, kind, r):
self.p, self.deg, self.e, self.kind, self.r = p, deg, e, kind, r
self.idx = None
def __repr__(self):
return "P(%d,f%d,e%d,%s%s)" % (self.p, self.deg, self.e, self.kind,
self.r)
class Field:
def __init__(self, form, dk=None):
if disc(form) >= 0 or (dk is not None and dk != disc(form)):
raise ValueError("needs disc(form) = d_K < 0 (maximal_form first)")
self.f = form
self.a, self.b, self.c, self.d = form
self.dk = disc(form) if dk is None else dk
a, b, c, d = form
rs = _cubic_roots(a, -b, c, -d) # roots of g
self.xi0 = min(rs, key=lambda z: abs(z.imag)).real
self.xi1 = max(rs, key=lambda z: z.imag)
self.emb0 = (1.0, a * self.xi0,
-(a * self.xi0 ** 2 - b * self.xi0 + c))
x1 = self.xi1
self.emb1 = (1.0, a * x1, -(a * x1 * x1 - b * x1 + c))
z = MPF(self.xi0)
for _ in range(8):
z = z - (((a * z - b) * z + c) * z - d) / ((3 * a * z - 2 * b)
* z + c)
self.hp0 = (a * z, -(a * z * z - b * z + c))
self._places = {}
self.cl = None
# --- arithmetic
def norm(self, x, y, z):
a, b, c, d = self.f
m00, m01, m02 = x, -a * c * y - a * d * z, -a * d * y - b * d * z
m10, m11, m12 = y, x + b * y, d * z
m20, m21, m22 = z, -a * y, x - c * z
return (m00 * (m11 * m22 - m12 * m21)
- m01 * (m10 * m22 - m12 * m20)
+ m02 * (m10 * m21 - m11 * m20))
def mul(self, u, v):
a, b, c, d = self.f
x1, y1, z1 = u
x2, y2, z2 = v
ww, wv, vv = y1 * y2, y1 * z2 + z1 * y2, z1 * z2
return (x1 * x2 - a * c * ww - a * d * wv - b * d * vv,
x1 * y2 + y1 * x2 + b * ww + d * vv,
x1 * z2 + z1 * x2 - a * ww - c * vv)
def s0(self, u):
e = self.emb0
return u[0] + u[1] * e[1] + u[2] * e[2]
def s1(self, u):
e = self.emb1
return u[0] + u[1] * e[1] + u[2] * e[2]
def vec(self, u):
z = self.s1(u)
return (self.s0(u), sqrt(2) * z.real, sqrt(2) * z.imag)
def logabs0(self, u):
"""log |u| at the real place, at mpmath's 60 digits."""
e1, e2 = self.hp0
return mpmath.log(abs(u[0] + u[1] * e1 + u[2] * e2))
# --- places
def places(self, p):
if p in self._places:
return self._places[p]
a, b, c, d = self.f
roots = [] # (kind, r, multiplicity)
for r in range(p):
if (((a * r - b) * r + c) * r - d) % p == 0:
roots.append(("f", r, _mult([a, -b, c, -d], r, p)))
if a % p == 0:
m = 1
if b % p == 0:
m = 2
if c % p == 0:
m = 3
roots.append(("i", 0, m))
tot = sum(m for _, _, m in roots)
out = []
if tot == 3:
for kind, r, m in roots:
out.append(Place(p, 1, m, kind, r))
elif tot == 1:
kind, r, m = roots[0]
out.append(Place(p, 1, 1, kind, r))
out.append(Place(p, 2, 1, "q", 0))
elif tot == 0:
out.append(Place(p, 3, 1, "n", 0))
else:
raise ValueError("impossible splitting")
self._places[p] = out
return out
def ptype(self, p):
"""The splitting pattern at p: the sorted (degree, ramification
index) of its places, read off places()."""
return tuple(sorted((pl.deg, pl.e) for pl in self.places(p)))
def hom(self, pl, k):
"""(alpha, beta) mod p^k: the images of w, v under O -> Z/p^k at a
degree-1 place (lifted by Newton where the root is simple)."""
a, b, c, d = self.f
p = pl.p
m = p ** k
if pl.kind == "f":
r = pl.r
if pl.e == 1:
for _ in range(k.bit_length() + 1):
gv = (((a * r - b) * r + c) * r - d) % m
if gv == 0:
break
dv = ((3 * a * r - 2 * b) * r + c) % m
r = (r - gv * pow(dv, -1, m)) % m
return (a * r) % m, -(a * r * r - b * r + c) % m
# the root at infinity: eta = 1/xi, h(eta) = -d eta^3 + c eta^2
# - b eta + a, w = b - c eta + d eta^2, v = -d eta
eta = 0
if pl.e == 1:
for _ in range(k.bit_length() + 1):
hv = (((-d * eta + c) * eta - b) * eta + a) % m
if hv == 0:
break
dv = ((-3 * d * eta + 2 * c) * eta - b) % m
eta = (eta - hv * pow(dv, -1, m)) % m
return (b - c * eta + d * eta * eta) % m, (-d * eta) % m
def valuations(self, u, p, n):
"""{place: v_P(u)} over the places above p, n = v_p(N(u)) >= 1.
u is primitive."""
pls = self.places(p)
out = {}
simple = [pl for pl in pls if pl.deg == 1 and pl.e == 1]
used = 0
for pl in simple:
al, be = self.hom(pl, n + 1)
m = p ** (n + 1)
val = (u[0] + u[1] * al + u[2] * be) % m
v = _vp(val, p) if val else n + 1
out[pl] = v
used += v
rest = [pl for pl in pls if pl not in out]
left = n - used
if len(rest) == 0:
if left != 0:
raise ValueError("valuation sum off")
elif len(rest) == 1:
pl = rest[0]
if left % pl.deg:
raise ValueError("valuation not divisible")
out[pl] = left // pl.deg
else:
raise ValueError("two unresolved places")
if any(v < 0 for v in out.values()):
raise ValueError("negative valuation")
return out
# --- lattices
def lll(self, basis):
"""LLL (delta 0.99) under T2 on integer coordinate vectors."""
B = [list(v) for v in basis]
n = len(B)
def emb(v):
return self.vec(v)
def dot(x, y):
return x[0] * y[0] + x[1] * y[1] + x[2] * y[2]
k = 1
it = 0
while k < n and it < 1000:
it += 1
E = [emb(v) for v in B]
Bs, mu = [], [[0.0] * n for _ in range(n)]
for i in range(n):
w = list(E[i])
for j in range(i):
mu[i][j] = dot(E[i], Bs[j]) / dot(Bs[j], Bs[j])
w = [w[t] - mu[i][j] * Bs[j][t] for t in range(3)]
Bs.append(w)
for j in range(k - 1, -1, -1):
q = round(mu[k][j])
if q:
B[k] = [B[k][t] - q * B[j][t] for t in range(3)]
E = [emb(v) for v in B]
for jj in range(j + 1):
mu[k][jj] -= q * (mu[j][jj] if jj < j else 1.0)
bk = dot(Bs[k], Bs[k])
bk1 = dot(Bs[k - 1], Bs[k - 1])
if bk >= (0.99 - mu[k][k - 1] ** 2) * bk1:
k += 1
else:
B[k], B[k - 1] = B[k - 1], B[k]
k = max(k - 1, 1)
return [tuple(v) for v in B]
# --- class group
def minkowski(self):
return (4 / pi) * (2 / 9) * sqrt(abs(self.dk))
def build_fb(self, bound):
self.B = bound
self.fb = []
for p in PRIMES:
if p > bound:
break
for pl in self.places(p):
if pl.deg == 3:
continue
pl.idx = len(self.fb)
self.fb.append(pl)
def smooth_row(self, u, N=None):
"""The exponent row of a primitive element u over the factor base,
or None if its norm is not smooth."""
if N is None:
N = self.norm(*u)
if N == 0:
return None
m = abs(N)
row = {}
for p in PRIMES:
if p > self.B:
break
if m % p:
continue
n = 0
while m % p == 0:
m //= p
n += 1
for pl, v in self.valuations(u, p, n).items():
if v:
if pl.idx is None:
return None
row[pl.idx] = v
if m == 1:
break
if m != 1:
return None
return row
def class_group(self, R_start=3, R_max=9, always_cert=False):
"""Read Cl, checked against the analytic formula; returns True,
or False if unresolved. Units are read
first from relations sharing an exponent row (their quotient is a
unit of small exact log), so every log carried through the
reduction is reduced modulo R from the start and stays small."""
self.always_cert = always_cert
self.floor_bad = False
bound = max(int(self.minkowski()) + 1, 30)
self.build_fb(bound)
n = len(self.fb)
self.hnf = {} # pivot col -> [row list, L]
self.units = []
self.R = MPF(0)
self.nrel = 0
self.val_bad = 0
first = {}
pending = []
for p in PRIMES: # the rational primes
if p > bound:
break
row = [0] * n
any_ = False
for pl in self.places(p):
if pl.idx is not None:
row[pl.idx] = pl.e
any_ = True
if any_:
pending.append((row, mpmath.log(p)))
base = self.lll([(1, 0, 0), (0, 1, 0), (0, 0, 1)])
self.base = base
seen = set()
for R in range(R_start, R_max + 1):
for i in range(-R, R + 1):
for j in range(-R, R + 1):
for k in range(-R, R + 1):
if max(abs(i), abs(j), abs(k)) != R and \
R != R_start:
continue
if gcd(gcd(i, j), k) != 1:
continue
if (-i, -j, -k) in seen:
continue
seen.add((i, j, k))
u = tuple(i * base[0][s] + j * base[1][s]
+ k * base[2][s] for s in range(3))
N = self.norm(*u)
try:
row = self.smooth_row(u, N)
except ValueError:
self.val_bad += 1
continue
if row is None:
continue
self.nrel += 1
L = self.logabs0(u)
key = tuple(sorted(row.items()))
if key in first:
du = abs(L - first[key])
if du > UNIT_TOL:
self._add_unit(du)
else:
first[key] = L
full = [0] * n
for c_, v in row.items():
full[c_] = v
pending.append((full, L))
for full, L in pending:
self._insert(full, L)
pending = []
if len(self.hnf) == n and self._certify():
return True
return False
def _insert(self, vec, L):
"""Insert a relation into the Hermite form, kept REDUCED (every entry
above a pivot below that pivot), so entries and logs stay small; a
relation reducing to zero is a unit."""
n = len(vec)
H = self.hnf
for j in range(n):
if vec[j] == 0:
continue
if j not in H:
if vec[j] < 0:
vec = [-x for x in vec]
L = -L
H[j] = list(self._hreduce(vec, L, j))
self._reduce_above(j)
return
piv, Lp = H[j]
if vec[j] % piv[j] == 0:
q = vec[j] // piv[j]
vec = [vec[i] - q * piv[i] for i in range(n)]
L = self._redL(L - q * Lp)
continue
g, s, t = xgcd(piv[j], vec[j])
u, w = piv[j] // g, vec[j] // g
newp = [s * piv[i] + t * vec[i] for i in range(n)]
newL = s * Lp + t * L
vec = [w * piv[i] - u * vec[i] for i in range(n)]
L = self._redL(w * Lp - u * L)
if newp[j] < 0:
newp = [-x for x in newp]
newL = -newL
H[j] = list(self._hreduce(newp, newL, j))
self._reduce_above(j)
if abs(L) > UNIT_TOL:
self._add_unit(abs(L))
def _hreduce(self, row, L, j):
"""Reduce the entries of row right of column j by the pivots."""
H = self.hnf
n = len(row)
for k in range(j + 1, n):
if row[k] and k in H:
pk = H[k][0]
q = row[k] // pk[k]
if q:
row = [row[i] - q * pk[i] for i in range(n)]
L = L - q * H[k][1]
return row, self._redL(L)
def _reduce_above(self, j):
H = self.hnf
for i in list(H):
if i < j and H[i][0][j]:
H[i] = list(self._hreduce(H[i][0], H[i][1], j - 1))
def _redL(self, L):
if self.R > 0:
return L - mpmath.nint(L / self.R) * self.R
return L
def _add_unit(self, L):
"""Fold a unit's log into R, the real gcd of the units' logs. The
logs are carried to 60 digits, so a remainder under UNIT_TOL is zero;
the drift is the largest distance of a unit's log from a multiple
of R, and a sound run keeps it at the precision's floor."""
self.units.append(L)
R = self.R
if R == 0:
self.R = L
return
a, b = max(R, L), min(R, L)
while b > UNIT_TOL * a:
r = a % b
if b - r < UNIT_TOL * a:
r = MPF(0)
a, b = b, r
self.R = a
def unit_drift(self):
R = self.R
return max((float(abs(u - mpmath.nint(u / R) * R))
for u in self.units), default=0.0)
def _certify(self):
h = 1
for j, (row, _) in self.hnf.items():
h *= row[j]
self.h_found = h
if h == 1 and not self.always_cert:
self.cert = None
self._snf()
return True
if self.R == 0:
return False
if self.R < R_FLOOR:
self.floor_bad = True
return False
L = self.euler()
est = sqrt(abs(self.dk)) * L / (2 * pi)
self.cert = h * float(self.R) / est
if self.cert > 1.5 or self.cert < 0.5:
return False
self._snf()
return True
def euler(self, X=30000):
"""The Euler product of zeta_K / zeta over p < X."""
a, b, c, d = self.f
dk = self.dk
prod = 1.0
for p in PRIMES:
if p >= X:
break
if dk % p == 0 or p == 2:
pat = self.ptype(p)
if pat == ((1, 1), (1, 2)):
prod *= 1 / (1 - 1 / p)
elif pat == ((1, 1), (1, 1), (1, 1)):
prod *= 1 / (1 - 1 / p) ** 2
elif pat == ((1, 1), (2, 1)):
prod *= 1 / (1 - 1 / (p * p))
elif pat == ((3, 1),):
prod *= (1 - 1 / p) / (1 - 1 / p ** 3)
continue
ls = legendre(dk, p)
if ls == -1:
prod *= 1 / (1 - 1 / (p * p))
continue
if self._splits(p):
prod *= 1 / (1 - 1 / p) ** 2
else:
prod *= (1 - 1 / p) / (1 - 1 / p ** 3)
return prod
def _splits(self, p):
"""At p not dividing d_K with (d_K/p) = +1: total splitting."""
a, b, c, d = self.f
if a % p == 0:
return True
ia = pow(a, -1, p)
e2, e1, e0 = (-b * ia) % p, (c * ia) % p, (-d * ia) % p
def mulmod(u, v):
r = [0] * 5
for i in range(3):
if u[i]:
for j in range(3):
r[i + j] += u[i] * v[j]
for k in (4, 3):
t = r[k] % p
if t:
r[k - 1] -= t * e2
r[k - 2] -= t * e1
r[k - 3] -= t * e0
return [r[0] % p, r[1] % p, r[2] % p]
res, base, e = [1, 0, 0], [0, 1, 0], p
while e:
if e & 1:
res = mulmod(res, base)
base = mulmod(base, base)
e >>= 1
return res == [0, 1, 0]
def _snf(self):
"""Smith form of the relation matrix with the column transform V:
a row vector e is in the lattice iff (e V)_k = 0 mod d_k."""
n = len(self.fb)
A = [list(self.hnf[j][0]) for j in range(n)]
V = [[int(i == j) for j in range(n)] for i in range(n)]
for k in range(n):
while True:
piv = None
for i in range(k, n):
for j in range(k, n):
if A[i][j] and (piv is None or abs(A[i][j]) <
abs(A[piv[0]][piv[1]])):
piv = (i, j)
if piv is None:
break
i, j = piv
A[k], A[i] = A[i], A[k]
for row in A:
row[k], row[j] = row[j], row[k]
for row in V:
row[k], row[j] = row[j], row[k]
done = True
pk = A[k][k]
for i in range(k + 1, n):
q = A[i][k] // pk
if q:
A[i] = [A[i][t] - q * A[k][t] for t in range(n)]
if A[i][k]:
done = False
for j in range(k + 1, n):
q = A[k][j] // pk
if q:
for row in A:
row[j] -= q * row[k]
for row in V:
row[j] -= q * row[k]
if A[k][j]:
done = False
if not done:
continue
bad = None
for i in range(k + 1, n):
for j in range(k + 1, n):
if A[i][j] % pk:
bad = i
break
if bad is not None:
break
if bad is None:
break
A[k] = [A[k][t] + A[bad][t] for t in range(n)]
diag = [abs(A[k][k]) for k in range(n)]
self.V = V
self.inv = [(k, dk) for k, dk in enumerate(diag) if dk > 1]
self.cl = tuple(dk for _, dk in self.inv)
self.h = 1
for dk in self.cl:
self.h *= dk
def cls(self, row):
"""The class of the ideal with exponent row {col: v}."""
out = []
for k, dk in self.inv:
s = 0
for c_, v in row.items():
s += v * self.V[c_][k]
out.append(s % dk)
return tuple(out)
def place_class(self, pl, R_max=6):
"""The class of a place in the factor base, or of a degree-1
place outside it, or None if not placed."""
if pl.idx is not None:
return self.cls({pl.idx: 1})
if pl.deg != 1:
raise ValueError("past the factor base, degree 1 only")
p = pl.p
al, be = self.hom(pl, 1)
lat = self.lll([(p, 0, 0), (-al, 1, 0), (-be, 0, 1)])
for R in range(1, R_max + 1):
for i in range(-R, R + 1):
for j in range(-R, R + 1):
for k in range(-R, R + 1):
if max(abs(i), abs(j), abs(k)) != R:
continue
u = tuple(i * lat[0][s] + j * lat[1][s]
+ k * lat[2][s] for s in range(3))
if gcd(gcd(u[0], u[1]), u[2]) != 1:
continue
N = self.norm(*u)
if N == 0 or N % p or (N // p) % p == 0:
continue
row = self.smooth_row(u, N // p)
if row is None:
continue
c = self.cls(row)
return tuple((-x) % dk for x, (_, dk) in
zip(c, self.inv))
return None