salem.py
Checks: COLLIDE: Signed pentanomials theorem.
Cited on: COLLIDE.
To run it: save salem.py; needs flint, mpmath, numpy.
"""salem.py -- the signed pentanomials a given Salem polynomial divides,
read from the factor's side: by pairs to a degree K, and by the gap tree
at every degree.
QUESTION. signed.py's survivors to n = 100 (signed f = 1 + s1 x^q + s2 x^p
+ s3 x^m + s4 x^n, not +-reciprocal, sharing a non-cyclotomic factor with
f~) are the dilations of three coprime ones, and their shared factors are
u8 = 1 + x^3 - x^4 + x^5 + x^8 and u10 = 1 - x + x^3 - x^4 + x^5 - x^6 +
x^7 - x^9 + x^10, which at -x are the degree-8 Salem polynomial (M =
1.280638) and Lehmer's (M = 1.176281). PAIRS: does either divide any
other signed pentanomial, to degree K far past the scan's reach? TREE:
which signed pentanomials, at any degree, does P(x) or P(-x) divide, P
the minimal polynomial of one of the 47 known Salem numbers below 1.3
(M. J. Mossinghoff's list; complete for degree <= 44 by M. J.
Mossinghoff, G. Rhin and Q. Wu, Minimal Mahler measures, Experiment.
Math. 17 (2008) 451-458, which lists every irreducible non-cyclotomic
polynomial of measure below 1.3 and degree at most 44)? A Salem P is
reciprocal, so P | f puts P in G = gcd(f, f~): a hit not +-reciprocal is
a survivor, and its primitive lies on a ray.
METHOD, PAIRS (section_pairs). A u dividing f puts every root of u mod a
prime P among f's roots mod P. For each u in {u8(x), u8(-x), u10(x),
u10(-x)} take two primes P < 2^31 with a root r of u mod P. The pair
values z = s r^a + t r^b mod P, 1 <= a < b <= K, s, t = +-1, sorted once;
f vanishes at r exactly when two pairs (q, p), (m, n) with p < m have
z1 + z2 = -1 mod P. Candidates pass the second prime, then exact division
over Z (flint). Every f u divides is found; a false candidate dies at
division.
METHOD, TREE (section_tree; the gap bound). Let b be a root of u with
0 < |b| < 1, L = -log|b|, e = (e1, e2, e3, e4 = n), sigma_0 = 1 and
sigma_j = 1 + sum_(i <= j) s_i b^(e_i). If u | f then sigma_4 = f(b) = 0,
so sigma_j = -sum_(i > j) s_i b^(e_i) and |sigma_j| < (4 - j)|b|^(e_(j+1))
for j < 3, while |sigma_3| = |b|^n. No sigma_j vanishes: sigma_1 = 0
needs |b|^e1 = 1, and sigma_2 = 0 leaves s3 b^e3 + s4 b^n = 0, |b|^(n -
e3) = 1. So e_(j+1) < log((4 - j)/|sigma_j|)/L at each level and n =
log(1/|sigma_3|)/L: the f are the leaves of a finite tree, and its list
is complete at every degree. Every non-cyclotomic u dividing a signed f
has such a root: its end coefficients are +-1, so its roots' product has
modulus 1, and Kronecker forbids them all on the circle. The tree takes
the root of least modulus (for a Salem P, 1/lambda at P(x), -1/lambda at
P(-x)). A leaf is tried at the nearest integer to n and its two
neighbours, kept when |sigma_3 + s4 b^n| < TOL |sigma_3|, TOL = 1e-20,
and decided by flint division over Z. A sigma below FLOOR = 1e-25 is
decided by exact division (an exact zero has no leaf below it) or counted
UNDECIDED, and an undecided node stops the run.
THE ERROR (inner_root, tree). Precision 50 digits, eps = mpmath's
machine epsilon, 2.7e-51, above the unit roundoff. Some root b* of u
lies within r = 2 D |u(b)/u'(b)| of the computed b, u(b) and u'(b)
evaluated at 100 digits, since |u'/u| = |sum 1/(b - root)| <= D/min|b -
root|; b* lies inside the disc with b, so the bound holds at b*. Each
power is one rounded product past the last, so |b^k computed - b*^k| <=
k (r + 3 eps) for k up to KMAX, the largest power the tree forms, and a
sigma, at most four powers added with four roundings, carries ERR <= 4
KMAX (r + 3 eps) + 20 eps. Then a true zero reads below FLOOR, so exact
division decides it; above FLOOR, the loop bounds move by about ERR/(L
FLOOR), and by r/(L |b|) times themselves through L, far below the one
extra step each loop runs (to floor(bound) + 1); and a true leaf reads
below TOL |sigma_3| whenever ERR < TOL FLOOR = 1e-45 (T4, each tree).
CONTROLS. S1 the scan's cross-check, an answer found another way: at n
<= 60 the hits not +-reciprocal (the scan's population) are exactly the
elements, divisible by some u, of the orbits under reversal and x -> -x
of the three primitives and their dilations to 60 (signed.py --top 60
found no other survivor). S2 at the full K each u's hits hold a known
primitive's orbit element. T0 the list as read: 47 rows, each rebuilt
reciprocal of degree D, irreducible, one root outside the circle, equal
to the listed lambda to 25 digits. T1 for u8(+-x), u10(+-x) the tree's
hits equal the pairs' at K = max(the tree's largest n, 60), with no
undecided node. T2 a plant off the Salem family: x^3 + x^2 + 1, not
reciprocal, its inner roots complex, divides 1 + x^3 - x^4 + x^7 + x^8;
the tree holds it and agrees with the pairs. T3 the kill's reduction to
primitives reads (2, 14, 18, 22) on the rays and (1, 2, 3, 5) off them.
T4 every tree's ERR below TOL FLOOR, read as each tree closes. Q2 no
undecided node over the sweep. SIGNS: the minus signs of the three
primitives' orbit elements and of their doubles' (x -> x^k keeps the
signs, and x -> -x acts on a dilation by k as on the primitive when k
is odd and trivially when k is even, so these two cover every degree).
PREDICTIONS. PAIRS (K3): the hits to K are the known primitives' orbits
and nothing else; K3 as first frozen read +-reciprocal f as hits too,
and they entered apart, outside the scan's population, after the first
print (a reciprocal u divides itself). TREE, frozen before the engine:
Q1 the four known u give 3, 3, 6, 6 hits, all of degree <= 19; Q2 the
largest n over all 94 trees is at most 200, with no undecided node; Q3
no primitive outside the three orbits, most of the other 45 Salem P
dividing no signed f that is not +-reciprocal.
KILL. PAIRS: "new hits outside the known orbits: h" (+-reciprocal f
apart) over u8(+-x), u10(+-x) at K = 1000: h >= 1 kills "these two divide
only the known rays' primitives". TREE: "primitives outside the three
orbits: h" over the 94 trees, after T0-T3: h >= 1 kills "three rays" (a
new ray, at any degree, its factor a small Salem polynomial).
FINDINGS, PAIRS (K = 1000). S1 12 hits at n <= 60, 12 wanted. u8(x)
and u8(-x) 3 hits each, u10(x) and u10(-x) 6 each, all of degree <= 19:
the known primitives' orbit elements, and +-reciprocal f apart, (3,
4, 5, 8) for u8 (1 + x^3 - x^4 + x^5 + x^8 = u8 itself, and its twin)
and (3, 7, 11, 14), (5, 6, 7, 12) for u10, among them the all-plus 1 +
x^3 + x^7 + x^11 + x^14 = (1 + x + x^2 + x^3 + x^4) u10(x), whose
non-cyclotomic part is Lehmer's polynomial at -x, irreducible as
Filaseta and Solan's question asks. "new hits outside the known
orbits: 0": the kill misses, to n = 1000. Each primitive found at K =
1000 is the one found at K = 60.
FINDINGS, TREE. 13 of 13 checks. T0 47 rows, bad []. T1 the tree's hits
are the pairs' at K = 60 for all four (u8: 320 nodes, largest n 32;
u10: 1182 nodes, largest n 51); T2 6 tree hits, 6 pair hits, b
complex, "4 nodes below FLOOR, undecided 0", the exact-division path
read; T4 "worst 2.39e-48 against 1.0e-45, over 99 trees". SIGNS [1,
2, 3, 4]: no ray element is 0/1, at any degree. Q1 held (3,
3, 6, 6, degrees <= 19). Q2 held: "94 trees; 90 hold no f that is not
+-reciprocal; largest n over all 73; undecided nodes 0". Q3 held:
"primitives outside the three orbits: 0 []", the kill misses. Hits
print for nine P (D <= 16), the rest none: the four trees with an f
not +-reciprocal are lambda = 1.176280 (Lehmer's, D = 10), whose four
are two from each orbit, (2, 10, 11, 15) and (1, 13, 16, 19) with
their mirrors, and lambda = 1.280638 (D = 8), whose two are (1, 7, 9,
11) and its mirror. The other seven P with hits (lambda = 1.200026,
1.202616, 1.216391, 1.230391, 1.236317, 1.261230, 1.293485) divide
only +-reciprocal f. So for each of the 47, at every degree, P(+-x)
divides a signed pentanomial not +-reciprocal only on the three rays'
primitives, never a dilation. "+-reciprocal hits with two
non-cyclotomic factors: 0 []": each has P as its one non-cyclotomic
factor, and four are 0/1, (3, 7, 11, 14) at Lehmer's, (1, 9, 17, 18)
at 1.216391, (3, 5, 7, 10) at 1.230391 and (1, 7, 13, 14) at 1.261230,
each with an irreducible non-cyclotomic part.
RUN RECORD. The bare run: 13 of 13, 24.3 s, peak commit 142.9 MB;
--rehearse 13 of 13, 6.7 s. T4 fails at 46 digits ("8.85e-45 against
1.0e-45"), stopping the run at the first tree. S1 fails when the
+-reciprocal hits are counted with the rest (18 hits, 12 wanted).
RESOURCE. 2 K^2 pair values of 8 bytes, sorted with their indices: about
100 MB at K = 1000; --rehearse K = 60, and the tree over the rows of
degree <= 10.
"""
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
from math import gcd # noqa: E402
import mpmath # noqa: E402
import numpy as np # noqa: E402
from flint import fmpz_poly, nmod_poly # noqa: E402
mpmath.mp.dps = 50
FLOOR = mpmath.mpf(10) ** -25
TOL = 1e-20
CHECKS = []
ERRS = []
REHEARSE = "--rehearse" in sys.argv
K = 60 if REHEARSE else (int(sys.argv[sys.argv.index("--top") + 1])
if "--top" in sys.argv else 1000)
U8 = [1, 0, 0, 1, -1, 1, 0, 0, 1]
U10 = [1, -1, 0, 1, -1, 1, -1, 1, 0, -1, 1]
PRIMITIVES = [((1, 7, 9, 11), (1, 1, 1, -1)),
((2, 10, 11, 15), (-1, -1, 1, -1)),
((1, 13, 16, 19), (1, 1, -1, -1))]
# M. J. Mossinghoff's list of the 47 known Salem numbers below 1.3: the
# degree D, the number (truncated), and the first D/2 + 1 coefficients
# of its reciprocal minimal polynomial.
SALEM = """\
10 1.176280818259917506544070338474 1 1 0 -1 -1 -1
18 1.188368147508223588142960958629 1 -1 1 -1 0 0 -1 1 -1 1
14 1.200026523987391518902962100414 1 0 0 -1 -1 0 0 1
14 1.202616743688604261118295415948 1 0 -1 0 0 0 0 -1
10 1.216391661138265091626806311199 1 0 0 0 -1 -1
18 1.219720859040311844169606760414 1 -1 0 0 0 0 0 0 -1 1
10 1.230391434407224702790177938975 1 0 0 -1 0 -1
20 1.232613548593121003962731694807 1 -1 0 0 0 -1 1 0 0 -1 1
22 1.235664580389747308105169351531 1 0 -1 -1 0 0 0 1 1 0 -1 -1
16 1.236317931803230489899094869802 1 -1 0 0 0 0 0 0 -1
26 1.237504821217490608171021829989 1 0 -1 0 0 -1 0 0 -1 0 1 0 0 1
12 1.240726423652541392056148161575 1 -1 1 -1 0 0 -1
18 1.252775937410113900864582824053 1 0 0 0 0 0 -1 -1 -1 -1
20 1.253330650201489757028162788986 1 0 -1 0 0 -1 0 0 0 0 0
14 1.255093516763722879173091003232 1 0 -1 -1 0 1 0 -1
18 1.256221154391670233067434043309 1 -1 0 0 -1 1 0 0 0 -1
24 1.260103540354990920321649852331 1 -1 0 0 -1 1 0 -1 1 -1 0 1 -1
22 1.260284236896492963739228435283 1 -1 0 -1 1 0 0 0 -1 1 -1 1
10 1.261230961137138851946671503074 1 0 -1 0 0 -1
26 1.263038139930169261222022798085 1 -1 0 0 0 0 -1 0 0 0 0 0 0 1
14 1.267296442523068692734077407604 1 -1 0 0 0 0 -1 1
22 1.276779674019016861136497157605 1 -1 -1 1 0 0 0 0 0 -1 0 1
8 1.280638156267757596701902532710 1 0 0 -1 -1
26 1.281691371528106310055107748672 1 0 0 0 0 0 -1 -1 -1 -1 -1 -1 -1 -1
20 1.282495560639960169561207128806 1 -2 2 -2 2 -2 1 0 -1 1 -1
18 1.284616550925536736743131441485 1 0 0 0 -1 0 -1 -1 0 -1
26 1.284746821544843035729838140650 1 -2 1 1 -2 1 0 0 -1 1 0 -1 1 -1
30 1.285099363651876557117420034512 1 0 0 0 0 -1 -1 -1 -1 -1 -1 0 0 0 0 1
30 1.285121520153207532780681369174 1 -2 2 -2 1 0 -1 2 -2 1 0 -1 1 -1 1 -1
30 1.285185670752909791356387310393 1 -1 0 0 0 0 0 0 -1 0 0 0 -1 0 0 -1
26 1.285196726769853432068127270005 1 0 -1 -1 0 0 0 1 0 -1 -1 0 1 1
44 1.285199179205612167312192383918 1 -1 0 0 0 0 0 -1 0 0 0 -1 0 0 0 0 0 0 0 1 0 0 1
30 1.285235436228923770828415884707 1 0 -1 0 0 -1 -1 0 0 0 1 0 0 1 0 -1
34 1.285409064765363764030309848277 1 -1 0 0 -1 1 -1 0 1 -1 1 0 -1 1 -1 0 1 -1
18 1.286395966836277224044411092745 1 -2 2 -2 2 -2 2 -3 3 -3
26 1.286730182048201274368747282841 1 -1 0 0 -1 1 -1 0 1 -1 1 0 -1 1
24 1.291741425714500483635599066106 1 -1 0 0 0 0 -1 0 0 0 0 0 0
20 1.292039106017929461943480560567 1 0 -1 0 0 -1 0 0 -1 0 1
40 1.292418657582426546281031229140 1 0 0 -1 0 -1 0 -1 0 -1 0 -1 0 0 1 0 1 0 1 0 1
46 1.292900721780102794630870596243 1 0 0 0 -1 -1 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
10 1.293485953125454106519909883794 1 0 -1 -1 0 1
18 1.295675371944048235295741653561 1 -1 0 0 -1 1 -1 0 1 -1
34 1.296210659593309216851783179125 1 -1 0 -1 0 1 0 1 -2 0 0 1 1 -1 -1 -1 1 1
22 1.296421365194547218873224266498 1 -1 0 0 0 -1 0 0 0 0 0 1
28 1.296821373714950077456125855369 1 0 0 0 -1 -1 -1 -1 -1 0 0 0 1 1 1
36 1.298429835475111538327805425740 1 1 0 -1 -2 -2 -1 0 1 1 0 -1 -1 0 1 1 0 -1 -1
26 1.299744869472170731620096386139 1 -1 -1 0 2 0 -2 -1 2 2 -2 -2 0 3
"""
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tag = "ok " if ok else "FAIL"
print(f" [{tag}] {name}" + (f" ({detail})" if detail else ""))
def control(name, ok, detail=""):
check(name, ok, detail)
if not ok:
print(" control failed: run stopped")
raise SystemExit(1)
def at_minus(c):
return [ci * (-1) ** i for i, ci in enumerate(c)]
def penta(e, s):
c = [0] * (e[-1] + 1)
c[0] = 1
for ei, si in zip(e, s):
c[ei] = si
return fmpz_poly(c)
def recip(e, s):
"""f is +-reciprocal: its normalized reversal is itself."""
return mirror(e) == e and (s[2] * s[3], s[1] * s[3], s[0] * s[3]) == s[:3]
def mirror(e):
q, p, m, n = e
return (n - m, n - p, n - q, n)
def orbit(e, s):
rs = (s[2] * s[3], s[1] * s[3], s[0] * s[3], s[3])
out = set()
for ee, ss in ((e, s), (mirror(e), rs)):
out.add((ee, ss))
out.add((ee, tuple(si * (-1) ** ei for ei, si in zip(ee, ss))))
return out
def roots_mod(u, count):
"""count primes P < 2^31, largest first, each with a root r of u."""
out = []
P = 2 ** 31 - 1
while len(out) < count:
P -= 2
if not is_prime(P):
continue
_, fl = nmod_poly(u, P).factor()
lin = [g for g, _ in fl if g.degree() == 1]
if lin:
c = lin[0].coeffs()
out.append((P, (-int(c[0]) * pow(int(c[1]), -1, P)) % P))
return out
def is_prime(n):
if n < 2:
return False
for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
if n % a == 0:
return n == a
d, r = n - 1, 0
while d % 2 == 0:
d //= 2
r += 1
for a in (2, 3, 5, 7, 11, 13, 17):
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(r - 1):
x = x * x % n
if x == n - 1:
break
else:
return False
return True
def hits(u, top):
"""Every signed pentanomial of degree <= top that u divides."""
(P, r), (P2, r2) = roots_mod(u, 2)
pw = np.array([pow(r, k, P) for k in range(top + 1)], dtype=np.int64)
a, b = np.triu_indices(top + 1, 1)
keep = a >= 1
a, b = a[keep].astype(np.int32), b[keep].astype(np.int32)
vals, A, B, S, T = [], [], [], [], []
for s in (1, -1):
for t in (1, -1):
vals.append((s * pw[a] + t * pw[b]) % P)
A.append(a)
B.append(b)
S.append(np.full(a.size, s, dtype=np.int8))
T.append(np.full(a.size, t, dtype=np.int8))
z = np.concatenate(vals)
A, B = np.concatenate(A), np.concatenate(B)
S, T = np.concatenate(S), np.concatenate(T)
del vals
order = np.argsort(z, kind="stable")
zs = z[order]
want = (-1 - z) % P
lo = np.searchsorted(zs, want, side="left")
hi = np.searchsorted(zs, want, side="right")
found = set()
for j in np.nonzero(hi > lo)[0]:
for i in order[lo[j]:hi[j]]:
if B[i] >= A[j]:
continue
e = (int(A[i]), int(B[i]), int(A[j]), int(B[j]))
s = (int(S[i]), int(T[i]), int(S[j]), int(T[j]))
if (1 + sum(si * pow(r2, ei, P2) for ei, si in zip(e, s))) % P2:
continue
if penta(e, s) % fmpz_poly(u) == 0:
found.add((e, s))
return found
def known(top):
"""The orbits of the three primitives' dilations to n = top."""
out = set()
for e, s in PRIMITIVES:
for k in range(1, top // e[-1] + 1):
ek = tuple(k * x for x in e)
out |= orbit(ek, s)
return out
def salem_rows():
"""(D, lambda as printed, full coefficient list) for each listed row."""
out = []
for line in SALEM.strip().splitlines():
w = line.split()
half = [int(c) for c in w[2:]]
out.append((int(w[0]), w[1], half + half[-2::-1]))
return out
def inner_root(u):
"""The root b of u of least modulus with 0 < |b| < 1, to mp.dps, and
r = 2 D |u(b)/u'(b)|, evaluated at twice the precision: some root of u
lies within r of b, since |u'/u| = |sum 1/(b - root)| <= D/min|b - root|."""
rs = mpmath.polyroots(u[::-1], maxsteps=800, extraprec=400)
b = min((r for r in rs if 0 < abs(r) < 1), key=abs)
b = mpmath.findroot(lambda x: mpmath.polyval(u[::-1], x), b)
b = b.real if abs(mpmath.im(b)) < FLOOR else b
with mpmath.workdps(2 * mpmath.mp.dps):
v, dv = mpmath.polyval(u[::-1], b, derivative=True)
r = 2 * (len(u) - 1) * abs(v) / abs(dv)
return b, r
def tree(u):
"""Every signed pentanomial u divides, at every degree: the gap tree.
With b = inner_root(u), L = -log|b| and sigma_j the partial sums of
f(b), a hit has e_(j+1) < log((4 - j)/|sigma_j|)/L for j = 0, 1, 2 and
|sigma_3| = |b|^n. A node whose |sigma| falls below FLOOR is decided by
exact division (an exact zero has no leaf below it) or counted
undecided. Returns (hits, nodes, largest n reachable, undecided,
nodes below FLOOR)."""
b, r = inner_root(u)
L = -mpmath.log(abs(b))
U = fmpz_poly(u)
pw = [mpmath.mpf(1)]
def p(k):
while len(pw) <= k:
pw.append(pw[-1] * b)
return pw[k]
def end(sig, c):
return int(mpmath.floor(mpmath.log(c / abs(sig)) / L)) + 1
def zero(e, s):
"""sigma at (e, s) is exactly 0: u divides 1 + sum s_i x^e_i."""
c = [0] * (e[-1] + 1)
c[0] = 1
for ei, si in zip(e, s):
c[ei] = si
return fmpz_poly(c) % U == 0
hits, nodes, top, undecided, low = set(), 0, 0, 0, 0
for e1 in range(1, end(1, 4) + 1):
for s1 in (1, -1):
sig1 = 1 + s1 * p(e1)
for e2 in range(e1 + 1, end(sig1, 3) + 1):
for s2 in (1, -1):
sig2 = sig1 + s2 * p(e2)
nodes += 1
if abs(sig2) < FLOOR:
low += 1
undecided += not zero((e1, e2), (s1, s2))
continue
for e3 in range(e2 + 1, end(sig2, 2) + 1):
for s3 in (1, -1):
sig3 = sig2 + s3 * p(e3)
nodes += 1
if abs(sig3) < FLOOR:
low += 1
undecided += not zero((e1, e2, e3),
(s1, s2, s3))
continue
nr = mpmath.log(1 / abs(sig3)) / L
mid = int(mpmath.nint(nr))
top = max(top, mid + 1)
for n in range(max(mid - 1, e3 + 1), mid + 2):
for s4 in (1, -1):
if (abs(sig3 + s4 * p(n))
>= abs(sig3) * TOL):
continue
e = (e1, e2, e3, n)
s = (s1, s2, s3, s4)
if penta(e, s) % U == 0:
hits.add((e, s))
eps = mpmath.mp.eps
err = 4 * (len(pw) - 1) * (r + 3 * eps) + 20 * eps
ERRS.append(err)
if not err < TOL * FLOOR:
control("T4 the tree's error clears the leaf tolerance", False,
f"{mpmath.nstr(err, 3)} against {mpmath.nstr(TOL * FLOOR, 3)}")
return hits, nodes, top, undecided, low
def primitive(e, s):
g = gcd(*e)
return tuple(x // g for x in e), s
def on_rays(e, s):
"""f's primitive lies in the orbit of one of the three."""
return primitive(e, s) in set().union(*(orbit(*x) for x in PRIMITIVES))
def section_tree():
print("TREE -- the gap tree: every signed pentanomial u divides")
rows = salem_rows()
bad = []
for d, lam, c in rows:
fl = fmpz_poly(c).factor()[1]
rs = mpmath.polyroots(c[::-1], maxsteps=800, extraprec=400)
out = [r for r in rs if abs(r) > 1 + mpmath.mpf(10) ** -20]
ok = (len(c) == d + 1 and len(fl) == 1 and fl[0][1] == 1
and len(out) == 1 and abs(mpmath.im(out[0])) < FLOOR
and mpmath.nstr(mpmath.re(out[0]), 32, strip_zeros=False)
[:27] == lam[:27])
if not ok:
bad.append((d, lam))
control("T0 the list as read: 47 rows, each irreducible of degree D "
"with one root outside the circle, lambda to 25 digits",
len(rows) == 47 and not bad, f"{len(rows)} rows, bad {bad[:2]}")
known4 = {"u8(x)": U8, "u8(-x)": at_minus(U8),
"u10(x)": U10, "u10(-x)": at_minus(U10)}
agree = []
for name, u in known4.items():
h, nodes, top, und, low = tree(u)
kk = max(top, 60)
hp = hits(u, kk)
agree.append(h == hp and und == 0)
print(f" {name}: {len(h)} tree hits, {len(hp)} pair hits to "
f"K = {kk}, {nodes} nodes, largest n {top}, undecided {und}, "
f"degrees <= {max(e[-1] for e, _ in h)}")
control("T1 the tree's hits are the pair method's, for u8(+-x), "
"u10(+-x)", all(agree), f"{agree}")
plant = [1, 0, 1, 1]
h, nodes, top, und, low = tree(plant)
hp = hits(plant, max(top, 60))
want = ((3, 4, 7, 8), (1, -1, 1, 1))
control("T2 plant x^3 + x^2 + 1: the tree holds 1 + x^3 - x^4 + x^7 + "
"x^8 and agrees with the pair method",
want in h and h == hp and und == 0, f"{len(h)} tree hits, {len(hp)} pair "
f"hits to K = {max(top, 60)}, b complex: "
f"{mpmath.im(inner_root(plant)[0]) != 0}, {low} nodes below "
f"FLOOR, undecided {und}")
control("T3 the reduction reads a dilation on the rays and a stranger "
"off them", on_rays((2, 14, 18, 22), (1, 1, 1, -1))
and not on_rays((1, 2, 3, 5), (1, 1, 1, 1)))
signs = {sum(si < 0 for si in ss) for e, s in PRIMITIVES for k in (1, 2)
for _, ss in orbit(tuple(k * x for x in e), s)}
print(f" minus signs on the rays' orbit elements, dilations by 1 and 2: "
f"{sorted(signs)}")
rows = rows if not REHEARSE else [r for r in rows if r[0] <= 10]
new, tops, unds, quiet, two, plus = [], [], 0, 0, [], []
for i, (d, lam, c) in enumerate(rows, 1):
for sign, u in (("x", c), ("-x", at_minus(c))):
h, nodes, top, und, low = tree(u)
tops.append(top)
unds += und
rec = [x for x in h if recip(*x)]
for e, s in rec:
nc = sum(m for g, m in penta(e, s).factor()[1]
if not g.is_cyclotomic())
two += [(e, s)] if nc >= 2 else []
plus += [(lam[:8], e)] if s == (1, 1, 1, 1) else []
other = sorted(x for x in h if not recip(*x))
off = [x for x in other if not on_rays(*x)]
new += off
quiet += not other
if h:
print(f" D = {d:2} lambda = {lam[:8]} u({sign}): "
f"{len(h)} hits, {len(rec)} +-reciprocal, "
f"{len(other)} others {other[:4]}, off the rays "
f"{off[:4]}, {nodes} nodes, largest n {top}")
print(f" {2 * len(rows)} trees; {quiet} hold no f that is not "
f"+-reciprocal; largest n over all {max(tops)}; "
f"undecided nodes {unds}")
control("Q2 every node decided", unds == 0)
control("T4 every tree's error bound clears the leaf tolerance",
max(ERRS) < TOL * FLOOR, f"worst {mpmath.nstr(max(ERRS), 3)} "
f"against {mpmath.nstr(TOL * FLOOR, 3)}, over {len(ERRS)} trees")
print(f" +-reciprocal hits with two non-cyclotomic factors: "
f"{len(two)} {two[:4]}; the 0/1 among the +-reciprocal hits: {plus}")
print(f" primitives outside the three orbits: {len(new)} {new[:6]}")
check("KILL no primitive outside the three orbits, at any degree",
not new)
def section_pairs():
us = {"u8(x)": U8, "u8(-x)": at_minus(U8),
"u10(x)": U10, "u10(-x)": at_minus(U10)}
print(f"SALEM -- signed pentanomials divisible by u, degree <= {K}")
small = {}
for name, u in us.items():
small[name] = hits(u, 60)
union60 = {x for x in set().union(*small.values()) if not recip(*x)}
want60 = {x for x in known(60) if any(
penta(*x) % fmpz_poly(u) == 0 for u in us.values())}
control("S1 at n <= 60 the hits are the known orbits divisible by u",
union60 == want60, f"{len(union60)} hits, {len(want60)} wanted, "
f"extra {sorted(union60 - want60)[:3]}, "
f"missing {sorted(want60 - union60)[:3]}")
new = set()
for name, u in us.items():
t1 = time.time()
h = hits(u, K) if K > 60 else small[name]
prim = sorted({x for x in h if x[0] in {e for e, _ in PRIMITIVES}
or mirror(x[0]) in {e for e, _ in PRIMITIVES}})
control(f"S2 {name}: a known primitive among the hits", prim,
f"{prim[:2]}")
rec = sorted(x for x in h if recip(*x))
extra = sorted(x for x in h - known(K) if not recip(*x))
new |= set(extra)
print(f" {name}: {len(h)} hits to n = {K}, {len(rec)} +-reciprocal, "
f"{len(extra)} others outside the known orbits {extra[:6]}, "
f"{time.time() - t1:.1f} s")
for e, s in rec[:12]:
print(f" +-reciprocal e = {e} s = {s}")
print(f" new hits outside the known orbits: {len(new)}")
check(f"K3 no hit outside the known orbits to n = {K}", not new)
def main():
t0 = time.time()
section_pairs()
print()
section_tree()
print()
print(f"{sum(CHECKS)} of {len(CHECKS)} checks, {time.time() - t0:.1f} s")
return 0 if all(CHECKS) else 1
if __name__ == "__main__":
sys.exit(main())