flake.py
8.0 kB · python · 214 lines
1import sys2import time3from fractions import Fraction4from math import lcm, sqrt56import numpy as np78CORNERS = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]], dtype=np.int64)9TINY = 1e-131011def cells(level):12 pts = np.zeros((1, 3), dtype=np.int64)13 for _ in range(level):14 pts = (2 * pts[:, None, :] + CORNERS[None, :, :]).reshape(-1, 3)15 return pts1617def face_pairs(pts, level):18 bits = level + 119 key = (pts[:, 0] << (2 * bits)) | (pts[:, 1] << bits) | pts[:, 2]20 rank = np.argsort(key)21 ordered = key[rank]22 src, dst = [], []23 for step in (1 << (2 * bits), 1 << bits, 1):24 want = key + step25 slot = np.clip(np.searchsorted(ordered, want), 0, key.size - 1)26 hit = ordered[slot] == want27 src.append(np.nonzero(hit)[0])28 dst.append(rank[slot[hit]])29 return np.concatenate(src), np.concatenate(dst)3031def layers(count, indptr, nbr):32 seen = np.zeros(count, dtype=bool)33 seen[0] = True34 front = np.array([0], dtype=np.int64)35 groups, parents = [front], [np.array([-1], dtype=np.int64)]36 while True:37 span = indptr[front + 1] - indptr[front]38 total = int(span.sum())39 if total == 0:40 break41 base = np.repeat(indptr[front], span)42 step = np.arange(total) - np.repeat(np.cumsum(span) - span, span)43 cand = nbr[base + step]44 came = np.repeat(front, span)45 fresh = ~seen[cand]46 cand, came = cand[fresh], came[fresh]47 if cand.size == 0:48 break49 seen[cand] = True50 groups.append(cand)51 parents.append(came)52 front = cand53 return groups, parents, int(seen.sum())5455def rooted(level):56 pts = cells(level)57 count = pts.shape[0]58 src, dst = face_pairs(pts, level)59 tail = np.concatenate([src, dst])60 deg = np.bincount(tail, minlength=count)61 indptr = np.zeros(count + 1, dtype=np.int64)62 np.cumsum(deg, out=indptr[1:])63 nbr = np.concatenate([dst, src])[np.argsort(tail, kind="stable")]64 groups, parents, reached = layers(count, indptr, nbr)65 label = np.concatenate(groups)66 place = np.empty(count, dtype=np.int64)67 place[label] = np.arange(count)68 offset = np.zeros(len(groups) + 1, dtype=np.int64)69 np.cumsum([g.size for g in groups], out=offset[1:])70 local = [np.array([], dtype=np.int64)]71 for d in range(1, len(groups)):72 local.append(place[parents[d]] - offset[d - 1])73 return {74 "count": count,75 "deg": deg[label],76 "offset": offset,77 "local": local,78 "tree": reached == count and int(src.size) == count - 1,79 "edges": (place[src], place[dst]),80 }8182def below(tree, shift):83 offset, degree = tree["offset"], tree["deg"].astype(np.float64)84 acc = np.zeros(tree["count"])85 neg = flat = 086 for d in range(offset.size - 2, -1, -1):87 lo, hi = offset[d], offset[d + 1]88 piv = degree[lo:hi] - shift - acc[lo:hi]89 small = np.abs(piv) < TINY90 hits = int(small.sum())91 if hits:92 flat += hits93 piv = np.where(small, np.where(piv >= 0.0, TINY, -TINY), piv)94 neg += int((piv < 0.0).sum())95 if d:96 acc[offset[d - 1]:offset[d]] += np.bincount(97 tree["local"][d], weights=1.0 / piv, minlength=int(lo - offset[d - 1])98 )99 return neg, flat100101def edge(tree, target, lo, hi):102 for _ in range(200):103 mid = 0.5 * (lo + hi)104 if mid <= lo or mid >= hi:105 break106 if below(tree, mid)[0] >= target:107 hi = mid108 else:109 lo = mid110 return hi111112def exact_pivots(tree, shift):113 offset, degree = tree["offset"], tree["deg"]114 acc = [Fraction(0)] * tree["count"]115 piv = [Fraction(0)] * tree["count"]116 neg = 0117 for d in range(offset.size - 2, -1, -1):118 lo, hi = int(offset[d]), int(offset[d + 1])119 up, home = int(offset[d - 1]), tree["local"][d]120 for v in range(lo, hi):121 here = Fraction(int(degree[v])) - shift - acc[v]122 piv[v] = here123 if here < 0:124 neg += 1125 if d:126 if here == 0:127 return neg, piv, False128 acc[up + int(home[v - lo])] += Fraction(1) / here129 return neg, piv, True130131def null_vector(tree, piv):132 offset = tree["offset"]133 val = [Fraction(0)] * tree["count"]134 val[0] = Fraction(1)135 for d in range(1, offset.size - 1):136 lo, hi = int(offset[d]), int(offset[d + 1])137 up, home = int(offset[d - 1]), tree["local"][d]138 for v in range(lo, hi):139 val[v] = val[up + int(home[v - lo])] / piv[v]140 scale = 1141 for x in val:142 scale = lcm(scale, x.denominator)143 return [int(x * scale) for x in val]144145def residual(tree, vec):146 res = [(int(tree["deg"][v]) - 4) * x for v, x in enumerate(vec)]147 for a, b in zip(*tree["edges"]):148 res[a] -= vec[b]149 res[b] -= vec[a]150 return max(abs(x) for x in res)151152def main():153 top = int(sys.argv[1]) if len(sys.argv) > 1 else 11154 exact_top = int(sys.argv[2]) if len(sys.argv) > 2 else 6155 vector_top = int(sys.argv[3]) if len(sys.argv) > 3 else 4156 print(f"flake band gap, code 23 base 2, levels 1..{top}")157 print("\n L N tree #(eig<2) #(eig<4) target flat lo(L)")158 los, his, above = {}, {}, {}159 for level in range(1, top + 1):160 clock = time.time()161 tree = rooted(level)162 target = 3 * 4 ** (level - 1)163 b2, t2 = below(tree, 2.0)164 b4, t4 = below(tree, 4.0)165 los[level] = edge(tree, target, 1.0, 3.0)166 his[level] = edge(tree, tree["count"], 4.0, 8.0)167 above[level] = tree["count"] - b4168 good = tree["tree"] and b2 == b4 == target and t2 == 0 and t4 == 1169 print(f"{level:5} {tree['count']:10} {'yes' if tree['tree'] else 'NO':>4} {b2:8} "170 f"{b4:8} {target:8} {t2 + t4:4} {los[level]:.15f} "171 f"[{time.time() - clock:.1f}s] {'PASS' if good else 'FAIL'}")172 print("\n L hi(L) #(eig>=4) 4^(L-1)")173 for level, upper in his.items():174 split = above[level] == 4 ** (level - 1)175 print(f"{level:5} {upper:.10f} {above[level]:9} {4 ** (level - 1):9} "176 f"{'PASS' if split else 'FAIL'}")177 if 2 in his:178 print(f"hi(2) against 3+sqrt(5): {his[2] - (3.0 + sqrt(5.0)):.1e}")179180 print("\n L 2-lo(L) ratio (2-lo)*8^L")181 prev, ratios = None, []182 for level, lo in los.items():183 defect = 2.0 - lo184 ratio = None if prev is None else prev / defect185 ratios += [] if ratio is None else [ratio]186 shown = "-" if ratio is None else f"{ratio:.4f}"187 print(f"{level:5} {defect:.12f} {shown:>6} {defect * 8.0 ** level:11.6f}")188 prev = defect189 climb = all(a < b < 8.0 for a, b in zip(ratios, ratios[1:]))190 print(f"ratios rise monotonically and stay under 8: {'PASS' if climb else 'FAIL'}")191 print(f"(2-lo)*8^L at L={top}: {(2.0 - los[top]) * 8.0 ** top:.6f}")192 print("\nexact rational elimination")193 print(" L N root pivot at 4 earlier zeros #(eig<2) #(eig<4-d) #(eig<4+d) mult")194 delta = Fraction(1, 10 ** 9)195 for level in range(1, exact_top + 1):196 tree = rooted(level)197 piv, clean = exact_pivots(tree, Fraction(4))[1:]198 n2 = exact_pivots(tree, Fraction(2))[0]199 low = exact_pivots(tree, Fraction(4) - delta)[0]200 high = exact_pivots(tree, Fraction(4) + delta)[0]201 good = clean and piv[0] == 0 and n2 == low == 3 * 4 ** (level - 1) and high - low == 1202 print(f"{level:5} {tree['count']:5} {str(piv[0]):>15} {'none' if clean else 'SOME':>13}"203 f" {n2:8} {low:10} {high:10} {high - low:4} {'PASS' if good else 'FAIL'}")204 print("\ninteger null vector at 4")205 print(" L N digits max|(Lap-4I)v|")206 for level in range(1, vector_top + 1):207 tree = rooted(level)208 vec = null_vector(tree, exact_pivots(tree, Fraction(4))[1])209 worst = residual(tree, vec)210 digits = len(str(max(abs(x) for x in vec)))211 print(f"{level:5} {tree['count']:5} {digits:6} {worst:14} "212 f"{'PASS' if worst == 0 else 'FAIL'}")213214main()