ladder.py
24.9 kB · python · 611 lines
1import os2import shutil3import subprocess4import sys5import tempfile6import time7from collections import Counter, defaultdict8from concurrent.futures import ProcessPoolExecutor9from itertools import groupby10from math import comb, isqrt1112from sympy import Poly, discriminant, resultant, symbols1314T, N = symbols("t n")15WORKERS = min(8, os.cpu_count() or 1)16CHUNKS = 25617DIMS = (2, 3, 4, 5, 6)18VERBS = {}1920# BOX2122def sig_bounds(dim, cap=None):23 return [min(comb(dim, j), cap if cap else comb(dim, j)) + 1 for j in range(1, dim + 1)]2425def box_size(bounds):26 out = 127 for b in bounds:28 out *= b29 return out3031def unrank(index, bounds):32 out = [1]33 for b in bounds:34 out.append(index % b)35 index //= b36 return tuple(out)3738def multiplicity(dim, sig):39 out = 140 for j, s in enumerate(sig):41 out *= comb(comb(dim, j), s)42 return out4344def weight(sig):45 return list(sig)4647# POLYNOMIALS4849def polymul(a, b):50 out = [0] * (len(a) + len(b) - 1)51 for i, u in enumerate(a):52 for j, v in enumerate(b):53 out[i + j] += u * v54 return out5556def polypow(a, k):57 out = [1]58 for _ in range(k):59 out = polymul(out, a)60 return out6162def trim(a):63 while len(a) > 1 and a[-1] == 0:64 a.pop()65 return a6667def poly_of(coeffs, var=T):68 return Poly(list(reversed(coeffs)), var, domain="ZZ")6970def factors(coeffs):71 body = trim(list(coeffs))72 if body == [0]:73 return 0, []74 lead, parts = poly_of(body).factor_list()75 out = []76 for g, mult in parts:77 out.append((tuple(int(c) for c in reversed(g.all_coeffs())), mult))78 return int(lead), sorted(out)7980def star(g):81 d = len(g) - 182 out = [0] * (d + 1)83 for j, c in enumerate(g):84 if c:85 term = polymul([0] * j + [c], polypow([1, 1], d - j))86 for i, v in enumerate(term):87 out[i] += v88 return out8990def fill_coeffs(dim, sig):91 out = [0] * (dim + 1)92 for j, s in enumerate(sig):93 if s:94 term = polymul([0] * j + [s], polypow([1, 1], dim - j))95 for i, v in enumerate(term):96 out[i] += v97 return out9899def corners(dim, code):100 return [tuple((c >> i) & 1 for i in range(dim)) for c in range(1 << dim) if (code >> c) & 1]101102def signature_of(dim, code):103 sig = [0] * (dim + 1)104 for c in corners(dim, code):105 sig[sum(c)] += 1106 return tuple(sig)107108def text(coeffs, var="n"):109 body = []110 for i, c in reversed(list(enumerate(coeffs))):111 if not c:112 continue113 head = "" if c == 1 and i else str(c)114 body.append(head + ("" if i == 0 else var if i == 1 else "{}^{}".format(var, i)))115 return " + ".join(body) if body else "0"116117# NORM118119def norm_report(dim):120 seen = {}121 bad = []122 lifted = 0123 linear = Counter()124 stray = 0125 for code in range(1 << (1 << dim)):126 sig = signature_of(dim, code)127 if sig not in seen:128 body = trim(list(sig))129 fill = fill_coeffs(dim, sig)130 m = len(body) - 1131 lead, parts = factors(sig)132 build = polypow([1, 1], dim - m)133 for g, mult in parts:134 for _ in range(mult):135 build = polymul(build, star(g))136 build = [lead * c for c in build]137 res = Poly(resultant(poly_of(body, T).as_expr(), N - T * (N + 1), T), N, domain="ZZ")138 res = polymul([int(c) for c in reversed(res.all_coeffs())], polypow([1, 1], dim - m))139 units = Counter(g[1] + 1 for g, mult in parts for _ in range(mult) if len(g) == 2 and g[0] == 1)140 wild = sum(1 for g, mult in parts if len(g) == 2 and g[0] != 1)141 seen[sig] = (build == fill, res == fill, m < dim, units, wild)142 row = seen[sig]143 if not (row[0] and row[1]):144 bad.append((code, sig))145 lifted += row[2]146 if sig[0] == 1:147 linear += row[3]148 stray += row[4]149 return len(seen), bad, lifted, linear, stray150151def disc_pairs(dim):152 ok = 0153 off = 0154 for index in range(box_size(sig_bounds(dim))):155 sig = unrank(index, sig_bounds(dim))156 for g, mult in factors(weight(sig))[1]:157 if len(g) >= 3:158 a = int(discriminant(poly_of(list(g), T)))159 b = int(discriminant(poly_of(star(list(g)), N)))160 ok += a == b161 off += a != b162 return ok, off163164def verb_norm():165 print("NORM FORM")166 print(" P(n) = (n+1)^(D-deg W) cont(W) prod over irreducible factors g of W of g*(n), g*(n) = (n+1)^deg g g(n/(n+1))")167 print(" cont(W) is the content, which is 1 on the origin-filled box; the constant is never the leading coefficient")168 print(" the resultant form is P(n) = (n+1)^(D-deg W) Res_t(W(t), n - t(n+1))")169 for dim in (2, 3, 4):170 sigs, bad, lifted, linear, stray = norm_report(dim)171 print(" D = {}: {} designs, {} signatures, factored form exact on {}, resultant form exact on {}".format(172 dim, 1 << (1 << dim), sigs, (1 << (1 << dim)) - len(bad), (1 << (1 << dim)) - len(bad)))173 print(" designs needing the (n+1)^(D-deg W) lift, s_D = 0: {} of {}".format(lifted, 1 << (1 << dim)))174 print(" origin-filled designs: linear factors are (a n + 1) with a counted {}".format(dict(sorted(linear.items()))))175 print(" origin-filled linear factors not of the form (a n + 1): {}".format(stray))176 ok, off = disc_pairs(dim)177 print(" disc(g) = disc(g*) over the {} origin-filled signatures, {} factor slots of degree at least 2 with g(1) != 0, mismatches {}".format(178 box_size(sig_bounds(dim)), ok, off))179 sig = (0, 2, 1)180 lead, parts = factors(sig)181 print(" origin empty breaks the law: signature {} has W = {} and fill {} = {}".format(182 sig, text(list(sig), "t"), text(fill_coeffs(2, sig)),183 " ".join("({})".format(text(star(list(g)))) for g, mult in parts for _ in range(mult))))184185VERBS["norm"] = verb_norm186187# SWEEP188189def label_of(top, signs):190 if top == 1:191 return "Q"192 if top == 2:193 return "quadratic " + ("mixed" if len(signs) > 1 else ("imaginary" if -1 in signs else "real"))194 return "degree {}".format(top)195196def sweep_chunk(job):197 dim, cap, lo, hi, want = job198 bounds = sig_bounds(dim, cap)199 sigs = Counter()200 codes = Counter()201 degrees = Counter()202 keep = defaultdict(set)203 for index in range(lo, hi):204 sig = unrank(index, bounds)205 parts = factors(sig)[1]206 top = 1207 signs = set()208 for g, mult in parts:209 d = len(g) - 1210 if d < 2:211 continue212 top = max(top, d)213 degrees[d] += mult214 if d == 2:215 signs.add(-1 if g[1] * g[1] - 4 * g[2] * g[0] < 0 else 1)216 if want and 2 <= d <= 5:217 keep[d].add(g)218 tag = label_of(top, signs)219 sigs[tag] += 1220 codes[tag] += multiplicity(dim, sig)221 return sigs, codes, degrees, {d: s for d, s in keep.items()}222223def sweep(dim, cap=None, want=False):224 total = box_size(sig_bounds(dim, cap))225 edges = [total * i // CHUNKS for i in range(CHUNKS + 1)]226 jobs = [(dim, cap, edges[i], edges[i + 1], want) for i in range(CHUNKS) if edges[i] < edges[i + 1]]227 sigs, codes, degrees = Counter(), Counter(), Counter()228 keep = defaultdict(set)229 if len(jobs) == 1:230 results = [sweep_chunk(jobs[0])]231 else:232 with ProcessPoolExecutor(max_workers=WORKERS) as pool:233 results = list(pool.map(sweep_chunk, jobs))234 for a, b, c, d in results:235 sigs += a236 codes += b237 degrees += c238 for key, value in d.items():239 keep[key] |= value240 return total, sigs, codes, degrees, keep241242def order(tag):243 return ("Q", "quadratic imaginary", "quadratic real", "quadratic mixed").index(tag) if tag.startswith(("Q", "quad")) else 4 + int(tag.split()[-1])244245def verb_ladder():246 print("THE LADDER BY DEGREE")247 print(" box: s_0 = 1, 0 <= s_j <= C(D, j); a signature carries C(D, j) choose s_j oriented designs")248 for dim in DIMS:249 total, sigs, codes, degrees, _ = sweep(dim)250 print(" D = {}: {} signatures, {} oriented designs".format(dim, total, sum(codes.values())))251 for tag in sorted(sigs, key=order):252 print(" {:20s} signatures {:8d} designs {}".format(tag, sigs[tag], codes[tag]))253 print(" irreducible factor slots by degree: {}".format(dict(sorted(degrees.items()))))254255VERBS["ladder"] = verb_ladder256257# TABLES258259TABLES = {260 (3, (1, 1)): [23, 31, 44, 59, 76, 83, 87, 104, 107, 108, 116, 135, 139, 140, 152, 172, 175, 199, 200, 204, 211, 212, 216, 231, 239, 243, 244, 247, 255, 268, 283, 300, 307, 324, 327, 331, 335, 339, 351, 356, 364, 367, 379, 411, 419, 424, 431, 436, 439, 440, 451, 459, 460, 472, 484, 491, 492, 499, 503, 515, 516, 519, 524, 527, 543, 547, 563, 567, 588, 620, 628, 643, 648, 652, 655, 671, 675, 676, 679, 680, 687, 695, 696, 707, 716, 728, 731, 743, 744, 748, 751, 755, 756, 759, 771, 780, 804, 808, 812, 815],261 (3, (3, 0)): [49, 81, 148, 169, 229, 257, 316, 321, 361, 404, 469, 473, 564, 568, 621, 697, 733, 756, 761, 785, 788, 837, 892, 940, 961, 985, 993, 1016, 1076, 1101, 1129, 1229, 1257, 1300, 1304, 1345, 1369, 1373, 1384, 1396, 1425, 1436, 1489, 1492, 1509, 1524, 1556, 1573, 1593, 1620, 1708, 1765, 1772, 1825, 1849, 1901, 1929, 1937, 1940, 1944, 1957, 2021, 2024, 2057, 2089, 2101, 2177, 2213, 2228, 2233, 2241, 2292, 2296, 2300, 2349, 2429, 2505, 2557, 2589, 2597, 2636, 2673, 2677, 2700, 2708, 2713, 2777, 2804, 2808, 2836, 2857, 2917, 2920, 2941, 2981, 2993, 3021, 3028, 3124, 3132],262 (4, (0, 2)): [117, 125, 144, 189, 225, 229, 256, 257, 272, 320, 333, 392, 400, 432, 441, 512, 513, 549, 576, 576, 592, 605, 656, 657, 697, 761, 784, 788, 832, 837, 873, 892, 981, 985, 1008, 1008, 1016, 1025, 1040, 1040, 1076, 1088, 1088, 1089, 1129, 1161, 1168, 1197, 1197, 1225, 1229, 1257, 1264, 1280, 1372, 1384, 1396, 1413, 1421, 1424, 1436, 1489, 1492, 1509, 1521, 1525, 1552, 1556, 1568, 1593, 1600, 1616, 1629, 1728, 1737, 1765, 1805, 1808, 1809, 1813, 1825, 1856, 1872, 1929, 1936, 1937, 1940, 1953, 1953, 2021, 2048, 2048, 2057, 2061, 2089, 2112, 2112, 2133, 2156, 2156],263 (4, (2, 1)): [275, 283, 331, 400, 448, 475, 491, 507, 563, 643, 688, 731, 751, 775, 848, 976, 1024, 1099, 1107, 1156, 1192, 1255, 1323, 1328, 1371, 1375, 1399, 1423, 1424, 1456, 1472, 1472, 1475, 1588, 1600, 1728, 1732, 1775, 1791, 1792, 1823, 1856, 1879, 1927, 1931, 1963, 1968, 1975, 1984, 1984, 2000, 2048, 2051, 2068, 2092, 2096, 2116, 2151, 2183, 2191, 2219, 2243, 2284, 2312, 2319, 2327, 2375, 2412, 2443, 2475, 2480, 2488, 2563, 2608, 2619, 2687, 2696, 2704, 2736, 2763, 2764, 2767, 2787, 2816, 2824, 2843, 2859, 2911, 2943, 3008, 3052, 3119, 3163, 3175, 3188, 3216, 3223, 3267, 3271, 3275],264 (4, (4, 0)): [725, 1125, 1600, 1957, 2000, 2048, 2225, 2304, 2525, 2624, 2777, 3600, 3981, 4205, 4225, 4352, 4400, 4525, 4752, 4913],265 (5, (1, 2)): [1609, 1649, 1777, 2209, 2297, 2617, 2665, 2869, 3017, 3089, 3233, 3369, 3857, 3889, 4169, 4261, 4409, 4417, 4429, 4432, 4477, 4549, 4597, 4757, 4817, 4897, 5025, 5164, 5437, 5501, 5584, 5653, 5753, 5864, 5913, 6241, 6449, 6581, 6757, 6793, 7096, 7177, 7265, 7333, 7373, 7376, 7672, 7684, 7717, 7909, 8073, 8105, 8249, 8329, 8357, 8529, 8705, 8752, 8945, 8968, 9065, 9137, 9412, 9437, 9489, 9552, 9584, 9664, 9701, 9808, 9829, 10229, 10277, 10329, 10381, 10449, 10492, 10532, 10589, 10609, 10729, 10825, 10832, 10933, 11317, 11332, 11469, 11693, 11809, 11876, 11993, 12184, 12205, 12349, 12389, 12440, 12481, 12517, 12533, 12752],266 (5, (3, 1)): [4511, 4903, 5519, 5783, 7031, 7367, 7463, 8519, 8647, 9439, 9759, 10407, 11119, 11243, 11551, 12447, 13219, 13523, 13799, 13883],267 (5, (5, 0)): [14641, 24217, 36497, 38569, 65657, 70601, 81509, 81589, 89417, 101833, 106069, 117688, 122821, 124817, 126032, 135076, 138136, 138917, 144209, 147109],268}269270QUADRATIC_REACH = 200271272def squarefree(m):273 m = abs(m)274 k = 2275 while k * k <= m:276 if m % (k * k) == 0:277 return False278 k += 1279 return True280281def fundamental(d):282 if d % 4 == 1:283 return squarefree(d)284 if d % 4 == 0:285 m = d // 4286 return m % 4 in (2, 3) and squarefree(m)287 return False288289def sign_of(sig):290 return (-1) ** sig[1]291292def table(degree, sig):293 if degree == 2:294 return [d for d in range(2, QUADRATIC_REACH) if fundamental(sign_of(sig) * d)]295 return TABLES.get((degree, sig))296297# FIELDS298299def square(value):300 if value <= 0:301 return 0302 root = isqrt(value)303 return root if root * root == value else 0304305def gp_text(coeffs):306 return "+".join("{}*x^{}".format(c, i) for i, c in enumerate(coeffs) if c)307308def gp_run(lines, per):309 if not lines:310 return []311 handle, name = tempfile.mkstemp(suffix=".gp")312 with os.fdopen(handle, "w") as fh:313 fh.write("\n".join(lines) + "\nquit\n")314 done = subprocess.run(["gp", "-q", "-s", "500000000", name], capture_output=True, text=True)315 os.unlink(name)316 values = [int(v) for v in done.stdout.split()]317 if len(values) != per * len(lines):318 raise RuntimeError("gp returned {} values for {} lines".format(len(values), len(lines)))319 return [values[i:i + per] for i in range(0, len(values), per)]320321def guarded(disc, raw):322 if not disc or raw % disc or disc % 4 not in (0, 1):323 return 0324 return square(raw // disc)325326def pari_fields(keys):327 lines = ["P={};print(nfdisc(P));print(polsturm(P));print(poldisc(P))".format(gp_text(list(reversed(g))))328 for g in keys]329 out = {}330 for g, (disc, r1, raw) in zip(keys, gp_run(lines, 3)):331 d = len(g) - 1332 index = guarded(disc, raw)333 out[g] = (d, (r1, (d - r1) // 2), disc if index else None, index, raw, "value" if index else "guard")334 return out335336def pari_classes(carriers, need):337 bodies = [gp_text(list(reversed(g))) for g in carriers]338 pairs = [(i, j) for i in range(len(bodies)) for j in range(i)]339 rows = gp_run(["print(nfisisom({},{})!=0)".format(bodies[i], bodies[j]) for i, j in pairs], 1)340 same = dict(zip(pairs, [row[0] for row in rows]))341 reps = []342 for i in range(len(bodies)):343 if all(not same[(i, j)] for j in reps):344 reps.append(i)345 if len(reps) >= need:346 break347 return len(reps)348349def field_data(keys):350 keys = sorted(keys)351 if not shutil.which("gp"):352 raise SystemExit("verb fields needs PARI, and gp is not on PATH")353 out = {}354 step = 20000355 for start in range(0, len(keys), step):356 out.update(pari_fields(keys[start:start + step]))357 return out358359def candidates(delta):360 out = set()361 f = 1362 while f * f <= abs(delta):363 if delta % (f * f) == 0 and (delta // (f * f)) % 4 in (0, 1):364 out.add(delta // (f * f))365 f += 1366 return out367368CACHE = {}369370def collect():371 if CACHE:372 return CACHE["data"], CACHE["reach"]373 per = {}374 for dim in DIMS:375 per[dim] = sweep(dim, want=True)[4]376 keys = set()377 for dim in DIMS:378 for d in per[dim]:379 keys |= per[dim][d]380 data = field_data(keys)381 reach = {}382 for dim in DIMS:383 classes = defaultdict(set)384 open_slots = defaultdict(list)385 members = defaultdict(list)386 for d in per[dim]:387 for g in per[dim][d]:388 degree, sig, disc, index, raw, mode = data[g]389 if disc is None:390 open_slots[(degree, sig)].append(g)391 else:392 classes[(degree, sig)].add(disc)393 members[(degree, sig)].append(g)394 reach[dim] = (classes, open_slots, members)395 CACHE["data"], CACHE["reach"] = data, reach396 return data, reach397398def risky(key, value, open_slots, data):399 return any(value in {abs(v) for v in candidates(data[g][4])} for g in open_slots.get(key, ()))400401def run_and_gap(reached, listed, key=None, open_slots=None, data=None):402 run = []403 gap = None404 owed = []405 for value in listed:406 if value in reached:407 run.append(value)408 continue409 if data is not None and risky(key, value, open_slots, data):410 owed.append(value)411 gap = value412 break413 return run, gap, owed414415def distinct_fields(carriers, need, cap=24):416 if not carriers:417 return 0, False418 return pari_classes(carriers[:cap], need), len(carriers) > cap419420def verb_fields():421 print("FIELD DISCRIMINANTS BY DEGREE AND SIGNATURE")422 print(" every distinct irreducible factor of W of degree 2..5 over the origin-filled box at D = 2..6")423 print(" field discriminant by PARI nfdisc on the reversed monic model x^d g(1/x), signature by polsturm, guarded against poldisc")424 print(" a value is accepted only if it is 0 or 1 mod 4 and divides the polynomial discriminant with a square quotient, the quotient being the square of the index")425 print(" the polynomial discriminant is used only for that guard, for the index and for the candidate bound on an unresolved factor")426 print(" tables of smallest field discriminants: degree 2 generated here, degrees 3, 4, 5 read from the LMFDB")427 data, reach = collect()428 print(" {} distinct irreducible factors of degree 2..5".format(len(data)))429 indexes = Counter(v[3] for v in data.values() if v[2] is not None and v[3])430 print(" index [O_K : Z[1/theta]]: {}".format(dict(sorted(indexes.items())[:8])))431 modes = Counter(v[5] for v in data.values())432 bad = sum(n for m, n in modes.items() if m != "value")433 print(" nfdisc values failing the guard: {} of {}".format(bad, len(data)))434 for dim in DIMS:435 classes, open_slots, members = reach[dim]436 print(" D = {}, box {}".format(dim, sig_bounds(dim)))437 for key in sorted(set(classes) | set(open_slots)):438 degree, sig = key439 listed = table(degree, sig)440 reached = {abs(v) for v in classes.get(key, ())}441 run, gap, owed = run_and_gap(reached, listed or [], key, open_slots, data)442 print(" degree {} signature {}: {} distinct field discriminants, {} unresolved".format(443 degree, sig, len(classes.get(key, ())), len(open_slots.get(key, ()))))444 if listed:445 print(" run {}".format(" ".join(str(v) for v in run) if run else "empty"))446 print(" first gap {}{}".format(447 gap if gap is not None else "none inside the table",448 "" if not owed else " (owed: an unresolved factor could carry {}, no probe)".format(owed)))449 layer_report(dim, classes)450451VERBS["fields"] = verb_fields452453# LAYER454455def fundamental_part(d):456 f = 1457 best = d458 while f * f <= abs(d):459 if d % (f * f) == 0 and fundamental(d // (f * f)):460 best = d // (f * f)461 f += 1462 return best463464def pure_layer(dim):465 out = defaultdict(set)466 for b in range(dim + 1):467 for c in range(1, comb(dim, 2) + 1):468 d = b * b - 4 * c469 root = int(abs(d) ** 0.5)470 if d >= 0 and root * root == d:471 continue472 out["imaginary" if d < 0 else "real"].add(fundamental_part(d))473 return out474475def layer_report(dim, classes):476 reach = 4 * comb(dim, 2)477 pure = pure_layer(dim)478 told = {"imaginary": {d for d in classes.get((2, (0, 1)), set())},479 "real": {d for d in classes.get((2, (2, 0)), set())}}480 window = {d for d in range(-reach, 0) if fundamental(d)}481 print(" reach 4 C(D, 2) = {}, fundamental discriminants inside it {}".format(reach, len(window)))482 for kind in ("imaginary", "real"):483 extra = sorted(told[kind] - pure[kind])484 short = sorted(pure[kind] - told[kind])485 print(" {}: pure layer {} fields, census {} fields, census only {}, pure only {}".format(486 kind, len(pure[kind]), len(told[kind]), extra if extra else "none", short if short else "none"))487 print(" imaginary pure layer equals the window: {}".format(pure["imaginary"] == window))488489# SWAP490491def splits(coeffs):492 return all(len(g) <= 2 for g, mult in factors(coeffs)[1])493494def swap_chunk(job):495 dim, lo, hi = job496 bounds = sig_bounds(dim)497 codes = Counter()498 sigs = Counter()499 for index in range(lo, hi):500 sig = unrank(index, bounds)501 void = [comb(dim, j) - sig[j] for j in range(1, dim + 1)]502 key = (splits(list(sig)), splits(void))503 codes[key] += multiplicity(dim, sig)504 sigs[key] += 1505 return codes, sigs506507def verb_swap():508 print("THE FILL AND VOID TABLE")509 print(" P(n) is the fill, V(n) = (2n+1)^D - P(n) the void, V(0) = 0 so V = n Q(n)")510 print(" P splits means P is a product of linear factors over Q, V splits means Q is")511 print(" the empty void, W = (1+t)^D, is counted as splitting")512 for dim in (3, 4):513 total = box_size(sig_bounds(dim))514 edges = [total * i // CHUNKS for i in range(CHUNKS + 1)]515 jobs = [(dim, edges[i], edges[i + 1]) for i in range(CHUNKS) if edges[i] < edges[i + 1]]516 codes, sigs = Counter(), Counter()517 with ProcessPoolExecutor(max_workers=WORKERS) as pool:518 for a, b in pool.map(swap_chunk, jobs):519 codes += a520 sigs += b521 designs = sum(codes.values())522 print(" D = {}: {} oriented designs with the origin filled, {} signatures".format(dim, designs, total))523 for key in ((True, True), (True, False), (False, True), (False, False)):524 tag = "P{} V{}".format("+" if key[0] else "-", "+" if key[1] else "-")525 print(" {} designs {:6d} share {:.4f} signatures {}".format(526 tag, codes[key], codes[key] / designs, sigs[key]))527 full = fill_coeffs(dim, [comb(dim, j) for j in range(dim + 1)])528 solid = polypow([1, 2], dim)529 print(" fill of the full box equals (2n+1)^D: {}".format(full == solid))530531VERBS["swap"] = verb_swap532533# HUNTER534535def merged(degree):536 rows = []537 limit = None538 sigs = sorted({k[1] for k in TABLES if k[0] == degree}) if degree > 2 else [(0, 1), (2, 0)]539 for sig in sigs:540 listed = table(degree, sig)541 seen = Counter()542 for value in listed:543 seen[value] += 1544 rows.append((value, sig, seen[value]))545 limit = listed[-1] if limit is None else min(limit, listed[-1])546 return sorted(rows), limit547548def witnessed(key, value, copy, dim, reach, data, fields):549 classes, open_slots, members = reach[dim]550 reached = {abs(v) for v in classes.get(key, ())}551 if value not in reached:552 return False, value if risky(key, value, open_slots, data) else None553 if copy == 1 or not fields:554 return True, None555 carriers = [g for g in members.get(key, ()) if abs(data[g][2]) == value]556 found, capped = distinct_fields(carriers, copy)557 return found >= copy, None558559def verb_hunter():560 print("THE HUNTER BOUND")561 print(" B is the largest bound with every field of degree d and abs discriminant at most B reached by the box")562 print(" the merge is signature aware: one discriminant in two signatures is two fields, and a discriminant listed twice in one signature is two fields resolved by field isomorphism")563 print(" the merged table is valid only below the smallest per-signature table maximum, printed per degree")564 print(" the box is s_0 = 1, 0 <= s_j <= C(D, j); its height is max_j C(D, j)")565 data, reach = collect()566 for degree in (2, 3, 4, 5):567 rows, limit = merged(degree)568 print(" degree {}: {} table entries, merge valid below {}".format(degree, len(rows), limit))569 for dim in (3, 4, 5, 6):570 marks = {}571 for fields in (True, False):572 bound = 0573 miss = None574 owed = None575 for value, group in groupby(rows, key=lambda row: row[0]):576 if value > limit or miss:577 break578 for value, sig, copy in group:579 ok, short = witnessed((degree, sig), value, copy, dim, reach, data, fields)580 if not ok:581 owed = short582 miss = (value, sig, copy)583 break584 if not miss:585 bound = value586 marks[fields] = (bound, miss, owed)587 field_bound, miss, owed = marks[True]588 disc_bound = marks[False][0]589 tail = "none inside the table" if miss is None else "{} at signature {}{}".format(590 miss[0], miss[1], " listed twice" if miss[2] > 1 else "")591 print(" D = {} box {} height {}: B by field {}, by discriminant {}, first miss {}{}".format(592 dim, sig_bounds(dim), max(comb(dim, j) for j in range(dim + 1)),593 field_bound or "none", disc_bound or "none", tail,594 "" if owed is None else " (owed: an unresolved factor could carry {}, no probe)".format(owed)))595596VERBS["hunter"] = verb_hunter597598# MAIN599600def main():601 asked = sys.argv[1:] or ["all"]602 if asked == ["all"]:603 asked = ["norm", "ladder", "fields", "swap", "hunter"]604 started = time.time()605 for name in asked:606 VERBS[name]()607 print()608 print("wall time {:.1f} s on {} workers".format(time.time() - started, WORKERS))609610if __name__ == "__main__":611 main()