"""flatten.py -- the least height of an integer polynomial forced to vanish to a given order at 1: the exact chart, where products of x^g - 1 attain it and where they do not, and the inverse dial read as the shortest integer-tap filter that kills polynomial signals. QUESTION. A vector c of M integers is the polynomial P(x) = sum c_i x^i of degree < M. Its FLATTENING is the multiplicity of the root 1 and its HEIGHT is max |c_i|. h(M, J) is the least height of a nonzero P of degree < M flattened to depth J. How small is it, what attains it, and read the other way, how short can a vector of taps bounded by k be if it must annihilate every polynomial signal of degree < J? THE OBJECTS. THE LATTICE. The vectors flattened to J are exactly (x - 1)^J Z[x] cut at degree < M, a lattice of rank r = M - J with basis x^i (x - 1)^J, i < r, so h(M, J) is a shortest vector in the sup norm. A CELL is a pair (M, J) with J < M. THE PURE PRODUCTS. prod (x^{g_i} - 1) over a multiset of positive PARTS g_i: it has degree sum g_i and flattening equal to the number of parts. One is ADMISSIBLE at a cell when its degree is below M and its part count at least J; ph(M, J) is the least height of an admissible one, so h <= ph, and a cell FAILS when h < ph. THE INVERSE DIAL. d_k(J) is the least degree of a nonzero multiple of (x - 1)^J of height at most k, so h(M, J) <= k iff M > d_k(J). A tap vector c annihilates every signal n -> n^j, j < J, under sum_i c_i s(n - i) exactly when sum_i c_i i^j = 0 for j < J, which is exactly flattening to depth J; so d_k(J) + 1 is the least number of taps of such a filter with integer taps bounded by k (property, by the moment identity below). THE HAND ATTACK, on paper before the engine. (1) THE MOMENT IDENTITY. P(x) = sum_j m_j (x - 1)^j with m_j = sum_i C(i, j) c_i, and sum_i C(i, j) c_i = 0 for all j < J is the same linear system as sum_i i^j c_i = 0 for j < J (the binomials C(i, j), j < J, and the powers i^j span the same polynomials in i of degree < J). So flattening is annihilation. (2) THE FREE RADIUS. A vector of sup norm at most H has Euclidean norm at most H sqrt(M), so the Euclidean ball of radius (q - 1) sqrt(M) holds every lattice vector of height below q. An enumeration of that ball, shrunk each time a vector of height below q is met, returns h exactly in exact arithmetic; its starting q is the least height among the reduced basis's own vectors, so the pure family under test supplies nothing to the number it is compared with. (3) RANK 1. At M = J + 1 the lattice is generated by (x - 1)^J alone, h = C(J, floor(J/2)), and the only pure product with J parts and degree J is (x - 1)^J, so ph = h there (property). (4) THE FINAL COEFFICIENTS. Multiplying P by a part x^g - 1 sends the coefficient at i < g to -c_i. So in a product built in nondecreasing part order, once the last part is g every coefficient below index g is final up to sign, and a walk may discard a partial product with such a coefficient above k. That makes the least degree of a J-part pure product of height <= k an exact branch-and-bound over partitions. (5) THE CYCLOTOMIC FLOOR. Let P have height <= k, degree <= T and (x - 1)^J | P, and let zeta be a primitive p^u-th root of unity, p prime. If Phi_{p^u} does not divide P, the norm of P(zeta) is a nonzero integer divisible by the norm of (zeta - 1)^J, which is p^J, while |P(zeta)| <= k(T + 1) at every conjugate; so J log p > phi(p^u) log(k(T + 1)) forces Phi_{p^u} | P. The forced factors are coprime to (x - 1)^J, so J + sum phi(p^u) <= T, and the least T passing that test is a floor on d_k(J) (property). (6) THE BOX PRINCIPLE. Over 0/1 vectors the moment m_j ranges over 0..C(M, j + 1); if 2^M exceeds prod_{j= 21. P3 RANK 1. h = ph = C(J, floor(J/2)) at every rank-1 cell. P4 THE DIAL. For k = 1..8 and J = 1..14 where the older table has a value, d_k(J) equals the least degree of a pure product of height <= k with exactly J parts, and the table reads J k=1 k=2 k=3 k=4 k=5 k=6 k=7 k=8 1 1 1 1 1 1 1 1 1 2 3 2 2 2 2 2 2 2 3 6 4 3 3 3 3 3 3 4 11 7 5 5 5 4 4 4 5 15 11 8 8 6 6 6 6 6 22 16 11 10 9 9 9 9 7 30 22 17 14 14 11 11 11 8 41 26 20 19 18 16 15 15 9 48 34 27 22 22 21 20 19 10 61 45 33 29 27 23 23 23 11 69 51 43 34 30 30 25 25 12 . . . 44 38 36 32 31 13 . . . . 46 46 37 37 14 . . . . . . 49 47 P5 THE FLOOR sits at or below every dial value read. P6 THE BOX. Every cell the box principle reaches has h = 1. KILLS, as prints. K1 A corner cell printing h < ph. K2 A dial cell where the lattice prints a vector of height <= k and degree below the walk's least pure degree: the pure law broken on the inverse dial (a finding, not a bug; inside the older table, whose cells were decided, it fails the run). KB (bugs, never findings) h > ph; a witness not cleared by J synthetic divisions (counted in the chart and break stages; the dial's K2 witness is printed, not counted; the chart's length and height tests hold by construction and cannot fire); a floor above a dial value, not tested at a K2 cell; a box cell with h > 1. CONTROLS, run before any verdict is read. C1 (POSITIVE, BRUTE FORCE) At every cell with M <= 6, h equals the least height found by exhausting the box [-3, 3]^M (or reports h > 3), and at M <= 9 over [-1, 1]^M the height-1 cells agree. C3 (POSITIVE, THE DIAL'S ROUTE) At every dial cell the decision route, run one width up, finds a vector of height <= k. C2 (POSITIVE, THE CLOSED FORM) The rank-1 cells match C(J, J//2) from the lattice route, which never sees the closed form; at rank 1 the lattice is its one basis vector, so C2 reads the reduction, not the enumeration, and prints only a count of mismatches. C4 (POSITIVE, THE SPLIT) The cyclotomic split finds Phi_6 beside 1 + 2x, Phi_14 and Phi_22 beside B = 2 + 3x + 2x^2, and Phi_22 alone, though deg Phi_n is not monotone in n: a loop stopping at the first Phi_n longer than what is left stops at Phi_5, Phi_11, Phi_17 and Phi_13 and finds none of the four. THE DESIGN. Exact integers throughout: an integral LLL (Cohen's fraction-free form, delta = 99/100) returns a reduced basis with its Gram-Schmidt data as integers; the enumeration runs Schnorr-Euchner over those data converted to floats with a relative slack of 1e-9 on the radius, and every vector it returns is recomputed exactly from integer coefficients, so a float error can only cost nodes, never admit a wrong vector; missing one would need an error above the slack. The pure products are one depth-first walk over the multisets of degree sum <= 39, keeping the least height per (part count, degree). Stages by argv: chart (C1, C2, C4, P1, P2, P3, P6), dial (P4, P5, the older table's cells) and extend (each column past the older table to J = 16, under a node cap per decision and a wall per column) and break (the one cell where the dial's pure law fails, decided on its own); all but extend by default. RESOURCE NOTE. Estimated: the chart a minute or less, a few MB; the dial dominated by its deepest decisions (k = 1 at J = 11 is a rank-58 enumeration), capped per decision by a node budget and printed as undecided rather than guessed; run under a memory guard at the 512 MB default. THE EXTENSION, designed after the dial stage's first run had printed the older table's 99 cells in seconds, where the older record had spent most of an hour: each column carried past the table to J = 16, under a cap of 2 x 10^7 enumeration nodes per decision and a wall of 240 s per column, the prediction being P4's law at every new cell and the kill K2 unchanged. On a K2 the stage steps the width down to the exact d_k(J) and prints the witness's cyclotomic factors and residual; the step-down runs under the same node cap with no handler, so a capped step-down ends the run with a traceback rather than printing it undecided. FINDINGS (copied from the printed output). F1 THE CONTROLS PASS. C1: 45 cells against brute force, 0 mismatches. C2: all 28 rank-1 cells print C(J, J//2) with h = ph. C3: the dial's route finds a vector of height <= k one width up, where the product sits, at all 99 cells of the default run and 115 with extend (15 more decided and the K2 cell), 0 failed. C4: the split finds all four, 0 off. F2 THE CORNER (P1 holds, K1 never fired): 63 cells, h = ph at all. F3 THE CHART (P2 holds to the cell): 695 cells, h = ph at 612 and h < ph at 83; by depth 3 of 362 (0.8%), 17 of 188 (9.0%) and 63 of 145 (43.4%); over the cells with ph > 1, 1.4% of 209 in the shallow band; at h >= 10, 2 of 61 (3.3%), 17 of 185 (9.2%) and 63 of 145 (43.4%); over ranks 5..18, 3 of 154 (1.9%), 17 of 112 (15.2%) and 63 of 104 (60.6%). The failing ranks are 5..18. The first failures by depth: h(23, 10) = 9 against 10, h(21, 12) = 39 against 44, h(25, 12) = 25 against 28. The cell (15, 6) has h = ph = 3 with parts (1, 1, 1, 2, 3, 3). F4 RANK 1 (P3 holds): the 28 rank-1 cells. F5 THE BOX (P6 holds): it reaches depth 1 at M = 4..8, depth 2 at M = 9..23 and depth 3 at M = 24..40, and every cell reached has h = 1. F6 THE DIAL (P4 holds at every older cell): all 99 cells decided and equal to the older table, the law holding at each. F7 THE FLOOR (P5 holds): below or at every value; equal at 27 cells, all at J <= 5; at k = 1 it reads 1, 3, 5, 7, 8, 11, 14, 15, 18, 25, 26 against the column's 1, 3, 6, 11, 15, 22, 30, 41, 48, 61, 69. F8 THE EXTENSION. Fifteen new cells decided, the law holding at each: d_2(13) = 74, d_2(15) = 95; d_3(12) = 52, d_3(16) = 87; d_4 = 52, 60, 70, 78 at J = 13..16; d_5 = 55, 64, 71 at J = 14..16; d_6 = 52, 60 at J = 14, 15; d_7 = 56, 64 at J = 15, 16. Thirteen decisions ran into the node cap, among them every k = 1 cell past J = 11, whose least pure degrees are 93, 112, 120, 141 and 159 at J = 12..16: upper bounds only. F9 K2 FIRED ONCE, AT k = 8, J = 15: d_8(15) = 55 against the least pure degree 56. The witness has height 8, clears 15 divisions, and splits as Phi_1^15 Phi_2^4 Phi_3^2 Phi_4^2 Phi_5^3 Phi_7 times the residual 2 + 3x + 3x^2 + 3x^3 + 4x^4 + 5x^5 + 4x^6 + 3x^7 + 3x^8 + 3x^9 + 2x^10, of value 35 at 1 and 1 at -1. Its cyclotomic part is the pure product with parts 1^5, 2, 2, 3, 3, 4, 4, 5, 5, 5, 7 (degree 45), so the witness is a pure product with ten degrees of it traded for a non-cyclotomic factor: the pure law on the inverse dial fails at height 8. F10 THE BREAK, DECIDED ON ITS OWN (added at audit, when the claim that no product with MORE parts reaches degree 55 was found to rest on nothing the walk had run): the least 15-part pure degree is 56, parts (1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 7, 9); a height-8 vector below degree 56 (3,090,666 nodes); none below degree 55 (4,984,378 nodes); and nothing of height <= 8 flattened to 16 below degree 56 (12,075 nodes), which covers every product with 16 or more parts. The split reprints F9's. HOLDS. RUN RECORD. chart: 2.7 s, peak working set 13.9 MB. dial: 4.2 s, peak 13.1 MB. extend: estimated up to 32 min (eight columns at a 240 s wall), cost 468.8 s on its first run and 481.1 s on the run carrying the step-down print, peak 15.6 MB; every figure identical across the two. break: 19.3 s; the default run (chart, dial, break) 28.6 s, peak 15.1 MB. All under a 512 MB memory guard. A later reading found the cyclotomic split stopping at the first Phi_n longer than its polynomial, though deg Phi_n is not monotone in n; it now tries every n whose totient fits (C4). The default run reprinted every figure, the break stage's split of the K2 cell among them, the one extend splits. members.py was not rerun: every factor it prints (A, B, Q, R, D, E, F, G, N) has leading coefficient above 1 or roots off the circle, so none is cyclotomic. """ import sys import time from math import comb, log, isqrt # ------------------------------------------------------------- polynomials def pow_xm1(J): """Coefficients of (x - 1)^J, lowest degree first.""" return [comb(J, i) * (-1) ** (J - i) for i in range(J + 1)] def height(p): return max(abs(c) for c in p) def clears(p, J): """True when (x - 1)^J divides p, by J synthetic divisions.""" q = list(p) for _ in range(J): if sum(q) != 0: return False # divide by (x - 1): q = (x - 1) s, s_{n-1} = q_n, s_{i-1} = q_i + s_i s = [0] * (len(q) - 1) acc = 0 for i in range(len(q) - 1, 0, -1): acc += q[i] s[i - 1] = acc q = s return True def basis(M, J): base = pow_xm1(J) rows = [] for i in range(M - J): v = [0] * M v[i:i + J + 1] = base rows.append(v) return rows # ----------------------------------------------------------- integral LLL def lll(B): """Cohen's integral LLL (Algorithm 2.6.7) at delta = 99/100. Returns the reduced basis, the integers d_0..d_n and lam[k][j] (j < k), with B*_k = d_{k+1}/d_k and mu_{k,j} = lam[k][j]/d_{j+1} (0-indexed).""" b = [list(v) for v in B] n = len(b) dot = lambda u, v: sum(x * y for x, y in zip(u, v)) d = [1] + [0] * n # d[i] for i = 0..n; d[i+1] belongs to b_i lam = [[0] * n for _ in range(n)] def redi(k, l): if 2 * abs(lam[k][l]) > d[l + 1]: q = (2 * lam[k][l] + d[l + 1]) // (2 * d[l + 1]) bl, bk = b[l], b[k] b[k] = [x - q * y for x, y in zip(bk, bl)] lam[k][l] -= q * d[l + 1] for i in range(l): lam[k][i] -= q * lam[l][i] def swapi(k): b[k], b[k - 1] = b[k - 1], b[k] for j in range(k - 1): lam[k][j], lam[k - 1][j] = lam[k - 1][j], lam[k][j] la = lam[k][k - 1] Bn = (d[k - 1] * d[k + 1] + la * la) // d[k] for i in range(k + 1, kmax[0] + 1): t = lam[i][k] lam[i][k] = (d[k + 1] * lam[i][k - 1] - la * t) // d[k] lam[i][k - 1] = (Bn * t + la * lam[i][k]) // d[k + 1] d[k] = Bn kmax = [0] d[1] = dot(b[0], b[0]) k = 1 while k < n: if k > kmax[0]: kmax[0] = k for j in range(k + 1): u = dot(b[k], b[j]) for i in range(j): u = (d[i + 1] * u - lam[k][i] * lam[j][i]) // d[i] if j < k: lam[k][j] = u else: if u == 0: raise ValueError("dependent basis") d[k + 1] = u redi(k, k - 1) la = lam[k][k - 1] if 100 * d[k + 1] * d[k - 1] < 99 * d[k] * d[k] - 100 * la * la: swapi(k) k = max(1, k - 1) else: for l in range(k - 2, -1, -1): redi(k, l) k += 1 return b, d, lam # ------------------------------------------------------------ enumeration class NodeCap(Exception): pass def enumerate_min(b, d, lam, M, s0, cap=None, stop_at=None): """Least sup norm of a nonzero lattice vector below s0, by Schnorr-Euchner enumeration of the Euclidean ball of radius (s - 1) sqrt(M), s the best height met so far. Returns (s, v, nodes) with v None when nothing below s0 exists. stop_at: return at the first vector of height <= stop_at.""" n = len(b) Bs = [d[i + 1] / d[i] for i in range(n)] mu = [[lam[k][j] / d[j + 1] for j in range(k)] for k in range(n)] best = [s0, None] R2 = [(s0 - 1) ** 2 * M * (1 + 1e-9) + 1e-9] nodes = [0] x = [0] * n def leaf(): v = [0] * M for i in range(n): xi = x[i] if xi: bi = b[i] for t in range(M): v[t] += xi * bi[t] s = max(abs(c) for c in v) if 0 < s < best[0]: best[0], best[1] = s, v R2[0] = (s - 1) ** 2 * M * (1 + 1e-9) + 1e-9 if stop_at is not None and s <= stop_at: raise StopIteration def rec(i, partial, top): # top: every coordinate above i is zero, so c = 0 and by the # sign symmetry only x_i >= 0 is walked, in order 0, 1, 2, ... # Otherwise the zigzag x0, x0 + s, x0 - s, x0 + 2s, ... is # ordered by distance from the centre c, so the first point # outside the ball ends the level. c = 0.0 if not top: for j in range(i + 1, n): if x[j]: c -= mu[j][i] * x[j] x0 = round(c) s = 1 if c >= x0 else -1 t = 0 while True: if top: xi = t elif t == 0: xi = x0 elif t % 2 == 1: xi = x0 + s * ((t + 1) // 2) else: xi = x0 - s * (t // 2) t += 1 nodes[0] += 1 if cap is not None and nodes[0] > cap: raise NodeCap y = xi - c p = partial + y * y * Bs[i] if p > R2[0]: break x[i] = xi if i == 0: if not (top and xi == 0): leaf() else: rec(i - 1, p, top and xi == 0) x[i] = 0 try: rec(n - 1, 0.0, True) except StopIteration: pass x[:] = [0] * n if best[1] is None: return None, None, nodes[0] return best[0], best[1], nodes[0] def least_height(M, J, cap=None): """h(M, J) exactly with a witness: reduce, start from the reduced basis's own least height, enumerate. Returns (h, v, nodes).""" b, d, lam = lll(basis(M, J)) H = min(height(v) for v in b) s, v, nodes = enumerate_min(b, d, lam, M, H, cap=cap) if s is None: # nothing below H: a reduced basis vector of height H is least v = min(b, key=height) s = H return s, v, nodes def exists_at_most(M, J, k, cap=None): """A nonzero vector of height <= k in the cell's lattice, or None.""" b, d, lam = lll(basis(M, J)) for v in b: if height(v) <= k: return v, 0 s, v, nodes = enumerate_min(b, d, lam, M, k + 1, cap=cap, stop_at=k) return v, nodes # ---------------------------------------------------------- pure products def mul_part(p, e): """p * (x^e - 1).""" q = [0] * (len(p) + e) for i, c in enumerate(p): q[i] -= c q[i + e] += c return q def pure_table(maxdeg): """Least height of a pure product per (part count, degree), with one witness multiset, over every multiset of degree sum <= maxdeg.""" best = {} def walk(p, parts, total, last): if parts: key = (len(parts), total) hgt = height(p) if key not in best or hgt < best[key][0]: best[key] = (hgt, tuple(parts)) for e in range(last, maxdeg - total + 1): parts.append(e) walk(mul_part(p, e), parts, total + e, e) parts.pop() walk([1], [], 0, 1) return best def ph_of(table, M, J): """ph(M, J) and an attaining multiset, or (None, None).""" out = (None, None) for (cnt, deg), (hgt, parts) in table.items(): if cnt >= J and deg < M and (out[0] is None or hgt < out[0]): out = (hgt, parts) return out def least_pure_degree(J, k, cap_deg): """Least degree of a pure product of height <= k with exactly J parts, by the final-coefficient walk (hand attack (4)); None when none has degree <= cap_deg.""" best = [cap_deg + 1, None] def walk(p, parts, total, last): if len(parts) == J: if height(p) <= k and total < best[0]: best[0], best[1] = total, tuple(parts) return need = J - len(parts) for e in range(last, best[0]): if total + need * e >= best[0]: break q = mul_part(p, e) if max(abs(c) for c in q[:e]) > k: continue parts.append(e) walk(q, parts, total + e, e) parts.pop() walk([1], [], 0, 1) if best[1] is None: return None, None return best[0], best[1] # ------------------------------------------------------ floor and the box def prime_powers(limit): out = [] for p in range(2, limit + 1): if all(p % q for q in range(2, isqrt(p) + 1)): q, e = p, 1 while q <= limit: out.append((p, e, q, q - q // p)) q *= p e += 1 return out def cyclotomic_floor(J, k): """Least T with J + (sum of the forced phi(p^u)) <= T.""" pps = prime_powers(4 * J * J + 50) N = J while True: L = log(k * (N + 1)) forced = sum(ph for (p, e, q, ph) in pps if J * log(p) > ph * L) if J + forced <= N: return N N += 1 def box_reach(M): """Largest J with 2^M > prod_{j nxt: prod, J = nxt, J + 1 else: break return J # ------------------------------------------------------ the cyclotomic split def pdivmod(a, b): """Quotient and remainder of integer polynomials, b monic.""" a = list(a) q = [0] * max(1, len(a) - len(b) + 1) for i in range(len(a) - len(b), -1, -1): c = a[i + len(b) - 1] q[i] = c if c: for j, bj in enumerate(b): a[i + j] -= c * bj while len(a) > 1 and a[-1] == 0: a.pop() return q, a def phi(n, _memo={}): if n not in _memo: p = [-1] + [0] * (n - 1) + [1] for d in range(1, n): if n % d == 0: p, _ = pdivmod(p, phi(d)) _memo[n] = p return _memo[n] def totients(N, _memo=[[0]]): """Euler's totient of 0..N, sieved once and kept.""" if len(_memo[0]) <= N: t = list(range(N + 1)) for i in range(2, N + 1): if t[i] == i: for j in range(i, N + 1, i): t[j] -= t[j] // i _memo[0] = t return _memo[0] def cyclotomic_split(v): """{n: multiplicity of Phi_n} and the residual, sign-normalised to a positive leading coefficient, with any power of x removed. Every n with deg Phi_n = totient(n) <= deg p is tried: the degree is not monotone in n (Phi_6 after Phi_5), and totient(n) >= sqrt(n / 2) bounds n by 2 deg^2.""" p = list(v) while p and p[0] == 0: p.pop(0) while p[-1] == 0: p.pop() mult = {} top = 2 * len(p) * len(p) tot = totients(top) for n in range(1, top + 1): if tot[n] >= len(p): continue f = phi(n) while len(f) <= len(p): q, r = pdivmod(p, f) if any(r): break p = q while len(p) > 1 and p[-1] == 0: p.pop() mult[n] = mult.get(n, 0) + 1 if p[-1] < 0: p = [-c for c in p] return mult, p # --------------------------------------------------------------- controls def brute_h(M, J, kmax): """Least height <= kmax of a nonzero flattened-to-J vector over the box [-kmax, kmax]^M, leading entry positive; None if none.""" from itertools import product best = None for v in product(range(-kmax, kmax + 1), repeat=M): nz = [c for c in v if c] if not nz or nz[-1] < 0: continue hv = max(abs(c) for c in v) if best is not None and hv >= best: continue if clears(list(v), J): best = hv return best # ------------------------------------------------------------------ stages def stage_chart(): t0 = time.time() table = pure_table(39) print("pure table: %d (count, degree) keys, %.1f s" % (len(table), time.time() - t0)) H = {} fails = [] kb = 0 nodes_all = 0 cells = [(M, J) for M in range(4, 41) for J in range(1, min(30, M - 1) + 1)] for (M, J) in cells: h, v, nodes = least_height(M, J) nodes_all += nodes p, parts = ph_of(table, M, J) if len(v) != M or height(v) != h or not clears(v, J) or h > p: kb += 1 print("KB bug at", M, J, h, p) H[(M, J)] = (h, p, parts, v) print("chart: %d cells (J = 1 included), %d nodes, %.1f s" % (len(cells), nodes_all, time.time() - t0)) # C1: brute force c1 = c1n = 0 for M in range(4, 7): for J in range(1, M): bh = brute_h(M, J, 3) c1n += 1 h = H[(M, J)][0] if (bh if bh is not None else 99) != (h if h <= 3 else 99): c1 += 1 print("C1 MISMATCH", M, J, h, bh) for M in range(4, 10): for J in range(1, M): bh = brute_h(M, J, 1) c1n += 1 if (bh == 1) != (H[(M, J)][0] == 1): c1 += 1 print("C1 MISMATCH h=1", M, J, H[(M, J)][0], bh) print("C1 brute force: %d cells, %d mismatches" % (c1n, c1)) # C2 and P3: rank 1 r1 = [(M, J) for (M, J) in H if M == J + 1] bad = [(M, J) for (M, J) in r1 if not (H[(M, J)][0] == H[(M, J)][1] == comb(J, J // 2))] print("C2/P3 rank 1: %d cells, %d off C(J, J//2) or h != ph" % (len(r1), len(bad))) # C4: the split tries every Phi_n that fits, deg Phi_n not monotone B = [2, 3, 2] def times(a, b): out = [0] * (len(a) + len(b) - 1) for i, x in enumerate(a): for j, y in enumerate(b): out[i + j] += x * y return out cases = [([1, 2], 6), (B, 14), (B, 22), ([1], 22)] c4 = [(f, n) for f, n in cases if cyclotomic_split(times(f, phi(n))) != ({n: 1}, f)] print("C4 the split finds Phi_n past a longer one, %d cases: %d off" % (len(cases), len(c4))) # P1: the corner corner = [(M, J) for M in range(4, 13) for J in range(1, M)] cf = [(M, J) for (M, J) in corner if H[(M, J)][0] < H[(M, J)][1]] print("P1 corner: %d cells, %d with h < ph%s" % (len(corner), len(cf), (" K1 " + str(cf)) if cf else "")) # P2: the chart chart = [(M, J) for (M, J) in H if J >= 2] fl = [(M, J) for (M, J) in chart if H[(M, J)][0] < H[(M, J)][1]] print("P2 chart: %d cells, h = ph at %d, h < ph at %d" % (len(chart), len(chart) - len(fl), len(fl))) for lo, hi in ((2, 12), (13, 20), (21, 30)): band = [c for c in chart if lo <= c[1] <= hi] bf = [c for c in band if c in fl] can = [c for c in band if H[c][1] > 1] big = [c for c in band if H[c][0] >= 10] bigf = [c for c in big if c in fl] print(" J = %2d..%2d: %3d of %3d fail (%.1f%%); over the %3d with" " ph > 1, %.1f%%; at h >= 10, %d of %d (%.1f%%)" % (lo, hi, len(bf), len(band), 100.0 * len(bf) / len(band), len(can), 100.0 * len(bf) / max(1, len(can)), len(bigf), len(big), 100.0 * len(bigf) / max(1, len(big)))) mid = [c for c in band if 5 <= c[0] - c[1] <= 18] midf = [c for c in mid if c in fl] print(" ranks 5..18: %d of %d (%.1f%%)" % (len(midf), len(mid), 100.0 * len(midf) / len(mid))) ranks = sorted(set(M - J for (M, J) in fl)) print(" failing ranks: %s" % ranks) print(" the first failures by depth:", [(M, J, H[(M, J)][0], H[(M, J)][1]) for (M, J) in sorted(fl, key=lambda c: (c[1], c[0]))[:6]]) # P6: the box bx = 0 reach = [] for M in range(4, 41): R = box_reach(M) reach.append(R) for J in range(1, R + 1): if (M, J) in H and H[(M, J)][0] != 1: bx += 1 print("P6 box: reach by M = 4..40 %s; cells off h = 1: %d" % (reach, bx)) print(" the cell (15, 6): h = %d, ph = %d, attaining parts %s" % (H[(15, 6)][0], H[(15, 6)][1], H[(15, 6)][2])) print("KB bugs: %d" % kb) return kb + c1 + len(bad) + len(cf) + bx + len(c4) OLD_TABLE = { 1: [1, 3, 6, 11, 15, 22, 30, 41, 48, 61, 69], 2: [1, 2, 4, 7, 11, 16, 22, 26, 34, 45, 51], 3: [1, 2, 3, 5, 8, 11, 17, 20, 27, 33, 43], 4: [1, 2, 3, 5, 8, 10, 14, 19, 22, 29, 34, 44], 5: [1, 2, 3, 5, 6, 9, 14, 18, 22, 27, 30, 38, 46], 6: [1, 2, 3, 4, 6, 9, 11, 16, 21, 23, 30, 36, 46], 7: [1, 2, 3, 4, 6, 9, 11, 15, 20, 23, 25, 32, 37, 49], 8: [1, 2, 3, 4, 6, 9, 11, 15, 19, 23, 25, 31, 37, 47], } def stage_dial(cap, jmax, colwall): t0 = time.time() rows = {} k2 = k2_known = kb = und = c3 = c3n = 0 for k in range(1, 9): out = [] J = 0 tcol = time.time() while True: J += 1 known = J <= len(OLD_TABLE[k]) if not known and (J > jmax or time.time() - tcol > colwall): break t1 = time.time() dp, parts = least_pure_degree(J, k, 400) if dp is None: out.append("none") continue # the lattice decides whether anything of degree < dp exists: # degree <= dp - 1 is the cell M = dp if dp == J: v, nodes = None, 0 # rank 0 below: nothing shorter else: try: v, nodes = exists_at_most(dp, J, k, cap=cap) except NodeCap: und += 1 out.append("%d?" % dp) print(" k=%d J=%2d: pure %d %s; lattice UNDECIDED at" " %d nodes (%.1f s)" % (k, J, dp, parts, cap, time.time() - t1)) continue # C3: the same route must find a vector of height <= k one # width up, where the pure product sits c3n += 1 try: w, _ = exists_at_most(dp + 1, J, k, cap=cap) except NodeCap: w = None if w is None or height(w) > k or not clears(w, J): c3 += 1 print(" C3 FAILED at k=%d J=%d" % (k, J)) if v is not None: k2 += 1 k2_known += known # step down to the exact d_k(J), then read the witness M = dp while True: w, _ = exists_at_most(M - 1, J, k, cap=cap) if w is None: break v, M = w, M - 1 mult, res = cyclotomic_split(v) print(" K2 k=%d J=%d: d_k(J) = %d against the least pure" " degree %d; witness height %d, cleared %s;" " cyclotomic part %s; residual %s, value %d at 1," " %d at -1" % (k, J, M - 1, dp, height(v), clears(v, J), sorted(mult.items()), res, sum(res), sum(c * (-1) ** i for i, c in enumerate(res)))) out.append("%d<%d" % (M - 1, dp)) continue fl = cyclotomic_floor(J, k) if fl > dp: kb += 1 print(" KB floor %d above %d at k=%d J=%d" % (fl, dp, k, J)) out.append(dp) rows[(k, J)] = (dp, parts, fl, nodes) print("k=%d: %s (%.1f s)" % (k, out, time.time() - t0)) old = [(k, J) for (k, J) in rows if J <= len(OLD_TABLE[k])] same = sum(1 for (k, J) in old if rows[(k, J)][0] == OLD_TABLE[k][J - 1]) print("C3 the route finds height <= k one width up: %d cells, %d failed" % (c3n, c3)) print("P4 dial: %d cells decided, %d of the %d older cells equal, %d K2," " %d undecided" % (len(rows), same, len(old), k2, und)) for k in range(1, 9): ext = [(J, rows[(k, J)][0], rows[(k, J)][1]) for J in range(1, 40) if (k, J) in rows and J > len(OLD_TABLE[k])] print(" past the older table, k=%d: %s" % (k, ext)) ties = [(k, J) for (k, J), r in rows.items() if r[2] == r[0]] print("P5 floor: KB %d; floor = value at %d cells: %s" % (kb, len(ties), sorted(ties))) for k in (1, 2): print(" floors k=%d: %s" % (k, [rows[(k, J)][2] for J in range(1, 12) if (k, J) in rows])) print(" the deepest decisions: %s" % sorted( ((r[3], k, J) for (k, J), r in rows.items()), reverse=True)[:5]) print(" attaining multisets, k = 2 and 3: %s" % [ (k, J, rows[(k, J)][1]) for k in (2, 3) for J in range(6, 12) if (k, J) in rows]) return kb + k2_known + c3 + (len(old) - same) def stage_break(cap): """The inverse dial's break at k = 8, J = 15, decided on its own: a height-8 vector below degree 56, nothing below degree 55, the walk's least 15-part degree, and no height-8 vector flattened to 16 below degree 56, which covers every product with more parts.""" t0 = time.time() k, J = 8, 15 dp, parts = least_pure_degree(J, k, 400) v, n1 = exists_at_most(56, J, k, cap=cap) none55, n2 = exists_at_most(55, J, k, cap=cap) deeper, n3 = exists_at_most(56, J + 1, k, cap=cap) ok = (dp == 56 and v is not None and none55 is None and deeper is None and height(v) <= k and clears(v, J)) mult, res = cyclotomic_split(v) if v is not None else ({}, []) print("break k=8 J=15: least 15-part pure degree %s %s; a height-%s" " vector of degree < 56 (%d nodes); none of degree < 55 (%s, %d" " nodes); none flattened to 16 below degree 56 (%s, %d nodes)" % (dp, parts, height(v) if v else None, n1, none55 is None, n2, deeper is None, n3)) print(" cyclotomic part %s, residual %s; %s (%.1f s)" % (sorted(mult.items()), res, "HOLDS" if ok else "FAILED", time.time() - t0)) return 0 if ok else 1 if __name__ == "__main__": t_all = time.time() stages = sys.argv[1:] or ["chart", "dial", "break"] bad = 0 if "chart" in stages: bad += stage_chart() if "dial" in stages: bad += stage_dial(cap=20000000, jmax=0, colwall=0) if "break" in stages: bad += stage_break(cap=200000000) if "extend" in stages: bad += stage_dial(cap=20000000, jmax=16, colwall=240) print("wall %.1f s" % (time.time() - t_all)) sys.exit(1 if bad else 0)