flatten.py
Checks: FLATTEN: The least height is a lattice minimum property; FLATTEN: At rank 1 the least height is the central binomial property; FLATTEN: In the corner the least height is a pure product's rule; FLATTEN: On the chart the products fail more often the deeper rule; FLATTEN: Through its table the inverse dial is a pure-product degree rule; FLATTEN: A failing cell is a break on the dial property; FLATTEN: At heights 8, 9 and 11 the dial breaks rule; FLATTEN: Two free bounds property.
Cited on: FLATTEN.
To run it: save flatten.py; standard library only.
"""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<J} (C(M, j + 1) + 1), two
vectors share their first J moments and their difference is a
height-1 vector flattened to J, so h(M, J) = 1 (property).
THE SLATE, frozen before the engine ran. Every figure below is a
TRANSPLANT from an older record of the same question and is read here
as a prediction, never as a fact.
P1 THE CORNER. At all 63 cells M = 4..12, J = 1..M-1, h = ph.
P2 THE CHART. Over the 695 cells M = 4..40, J = 2..min(30, M - 1),
h = ph at 612 and h < ph at 83; the failures fall 3 of 362 at
J <= 12, 17 of 188 at J = 13..20 and 63 of 145 at 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<J} (C(M, j+1) + 1), 0 if none."""
prod, J = 1, 0
while J < M:
nxt = prod * (comb(M, J + 1) + 1)
if 2 ** M > 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)