occupancy.py
10.8 kB · python · 306 lines
1import argparse2import time3from math import ceil, floor, gcd, log45import numpy as np67LOW = 0.44759788HIGH = 0.6402122910def up(x, d):11 return ceil(x * 10 ** d) / 10 ** d1213def down(x, d):14 return floor(x * 10 ** d) / 10 ** d1516# GASKET CHUNKS1718def base_pair(m, fix0):19 u = np.zeros(1, dtype=np.int64)20 v = np.zeros(1, dtype=np.int64)21 start = 022 if fix0:23 v = np.ones(1, dtype=np.int64)24 start = 125 for i in range(start, m):26 p = 3 ** i27 u = np.concatenate([u, u + p, u])28 v = np.concatenate([v, v, v + p])29 return u, v3031def top_offsets(lo, hi):32 out = [(0, 0)]33 for i in range(lo, hi):34 p = 3 ** i35 out = [(a, b) for a, b in out] + [(a + p, b) for a, b in out] + [(a, b + p) for a, b in out]36 return out3738def chunks(n, fix0, m):39 m = min(m, n)40 bu, bv = base_pair(m, fix0)41 for du, dv in top_offsets(m, n):42 yield bu + du, bv + dv4344# RATIO SET4546def ratio_values(k, m):47 mod = 3 ** k48 for u, v in chunks(k, True, m):49 keep = u > 050 u = u[keep]51 v = v[keep]52 vs = np.unique(v)53 iv = np.array([pow(int(x), -1, mod) for x in vs], dtype=np.int64)54 yield (u * iv[np.searchsorted(vs, v)]) % mod5556def ratio_stats(k, m, full):57 mod = 3 ** k58 if full:59 hist = np.zeros(mod, dtype=np.int32)60 for r in ratio_values(k, m):61 hist += np.bincount(r, minlength=mod).astype(np.int32)62 size = int(np.count_nonzero(hist))63 h = hist.astype(np.int64)64 return size, int((h * h).sum()), int(h.max())65 parts = []66 for r in ratio_values(k, m):67 parts.append(np.unique(r))68 if len(parts) > 24:69 parts = [np.unique(np.concatenate(parts))]70 return int(np.unique(np.concatenate(parts)).size), 0, 07172def cmd_ratios(args):73 print("k R_k sigma_k growth c_k k*c_k M2/4^k mean CSlower top secs")74 prev = 075 band = []76 for k in range(2, args.kmax + 1):77 t = time.time()78 size, second, top = ratio_stats(k, args.mem, k <= args.full)79 pairs = 3 ** (k - 1) - 2 ** (k - 1)80 growth = size / prev if prev else 0.081 ck = 1 - log(growth) / log(3) if growth else 0.082 if growth:83 band.append((k, ck))84 prev = size85 row = [k, size, round(size / 3 ** k, 6), round(growth, 4), round(ck, 7),86 round(k * ck, 7)]87 if second:88 row += [round(second / 4 ** k, 4), round(pairs / size, 4),89 round(pairs * pairs / second / 3 ** k, 6), top]90 print(*row, round(time.time() - t, 1))91 if len(band) > 1:92 prod = [k * c for k, c in band[-6:]]93 print("decay band k*c_k", down(min(prod), 4), up(max(prod), 4),94 "over k", band[-6:][0][0], band[-1][0])9596# CENSUS9798POW3 = np.array([3 ** i for i in range(39)], dtype=np.int64)99100def census(n, m, cut_pow, alphas, caps=()):101 thr = np.array([3.0 ** (a * n) for a in alphas] + [3.0 ** j for j in caps])102 keys = []103 pts = 0104 fibre = 0105 weighted = np.zeros(n + 2, dtype=np.int64)106 heavy = np.zeros(thr.size, dtype=np.int64)107 cut = max(3.0 ** cut_pow, float(thr.max()) if thr.size else 0.0)108 for u, v in chunks(n, False, m):109 ok = (u > 0) & (v > 0)110 u = u[ok]111 v = v[ok]112 fibre += int(ok.size - u.size)113 g = np.gcd(u, v)114 z1 = u // g115 z2 = v // g116 h = np.maximum(z1, z2)117 pts += h.size118 oct_ = np.searchsorted(POW3, h, side="right") - 1119 weighted += np.bincount(oct_, minlength=n + 2)[: n + 2]120 for i, t in enumerate(thr):121 heavy[i] += int(np.count_nonzero(h <= t))122 sel = h <= cut123 if sel.any():124 keys.append(np.unique((z1[sel] << np.int64(32)) + z2[sel]))125 if len(keys) > 48:126 keys = [np.unique(np.concatenate(keys))]127 keys = np.unique(np.concatenate(keys)) if keys else np.zeros(0, dtype=np.int64)128 kh = np.maximum(keys >> np.int64(32), keys & np.int64((1 << 32) - 1))129 occ = [int(np.count_nonzero(kh <= t)) for t in thr]130 return pts, fibre, weighted, heavy, occ, keys, kh131132def cmd_census(args):133 alphas = args.alphas134 caps = [j for j in args.caps]135 print("convention height max(z1,z2), window height <= 3^(alpha n), octave floor(log_3 h),"136 " desk octave = octave + 1 above h = 1, ray totals exclude the two fibre rays")137 print("n alpha A F meanM expA expF theta box A/box")138 seen = {}139 for n in args.levels:140 t = time.time()141 cutp = min(n, int(args.cut * n) + 1)142 pts, fibre, weighted, heavy, occ, keys, kh = census(n, args.mem, cutp, alphas, caps)143 for i, a in enumerate([str(x) for x in alphas] + ["3^%d" % j for j in caps]):144 x = 3.0 ** (alphas[i] * n) if i < len(alphas) else 3.0 ** caps[i - len(alphas)]145 e_a = log(occ[i]) / (n * log(3)) if occ[i] else 0.0146 th = log(occ[i]) / log(x) if occ[i] else 0.0147 box = int(x) ** 2148 seen.setdefault(a, []).append((n, e_a, th))149 print(n, a, occ[i], int(heavy[i]), round(heavy[i] / occ[i], 4) if occ[i] else 0.0,150 round(e_a, 4), round(log(int(heavy[i])) / (n * log(3)), 4) if heavy[i] else 0.0,151 round(th, 4), box, round(occ[i] / box, 6))152 print("level", n, "nonfibre", pts, "fibre", fibre, "cut", cutp, "keys", keys.size,153 "total_rays", keys.size if cutp >= n else 0,154 keys.size + 2 if cutp >= n else 0, round(time.time() - t, 1))155 print("octaves", n, *[int(x) for x in weighted[: n + 1]])156 for a, rows in seen.items():157 print("band", a, "expA", down(min(r[1] for r in rows), 4), up(max(r[1] for r in rows), 4),158 "theta", down(min(r[2] for r in rows), 4), up(max(r[2] for r in rows), 4),159 "over n", rows[0][0], rows[-1][0])160161# LARGE PRIME SUM162163def spf_sieve(limit):164 spf = np.zeros(limit + 1, dtype=np.int32)165 spf[2::2] = 2166 for p in range(3, int(limit ** 0.5) + 1, 2):167 if spf[p] == 0:168 spf[p * p:: 2 * p] = np.where(spf[p * p:: 2 * p] == 0, p, spf[p * p:: 2 * p])169 odd = np.arange(3, limit + 1, 2)170 spf[3::2] = np.where(spf[3::2] == 0, odd, spf[3::2])171 spf[1] = 1172 return spf173174def gcd_hist(n, m):175 hist = np.zeros(3 ** n, dtype=np.int32)176 for u, v in chunks(n, False, m):177 g = np.gcd(u, v)178 hist += np.bincount(g, minlength=3 ** n).astype(np.int32)179 hist[0] = 0180 return hist181182def cmd_sieve(args):183 print("n beta primesum F bound ratio")184 for n in args.levels:185 t = time.time()186 hist = gcd_hist(n, args.mem)187 spf = spf_sieve(3 ** n)188 vals = np.nonzero(hist)[0]189 mult = hist[vals].astype(np.int64)190 for beta in args.betas:191 thr = 3.0 ** (beta * n)192 work = vals.copy()193 big = np.zeros(work.size, dtype=np.int64)194 last = np.zeros(work.size, dtype=np.int64)195 while True:196 live = work > 1197 if not live.any():198 break199 p = spf[work].astype(np.int64)200 new = live & (p > thr) & (p != last)201 big += new202 last = np.where(live, p, last)203 work = np.where(live, work // np.maximum(p, 1), work)204 total = int((big * mult).sum())205 _, _, _, heavy, occ, _, _ = census(n, args.mem, 1, [1.0 - beta])206 bound = (int(heavy[0]) + 2 ** (n + 1)) / beta207 print(n, beta, total, int(heavy[0]), round(bound, 1), round(total / bound, 4))208 print("level", n, "secs", round(time.time() - t, 1))209210# CHECKS211212def brute_ratios(k):213 mod = 3 ** k214 seen = set()215 for lab in range(3 ** k):216 u = v = 0217 w = lab218 for i in range(k):219 d = w % 3220 w //= 3221 if d == 1:222 u += 3 ** i223 elif d == 2:224 v += 3 ** i225 if u and v % 3:226 seen.add(u * pow(v, -1, mod) % mod)227 return len(seen)228229def brute_rays(n, xcap):230 out = {}231 for lab in range(3 ** n):232 u = v = 0233 w = lab234 for i in range(n):235 d = w % 3236 w //= 3237 if d == 1:238 u += 3 ** i239 elif d == 2:240 v += 3 ** i241 if u and v:242 g = gcd(u, v)243 z = (u // g, v // g)244 if max(z) <= xcap:245 out[z] = out.get(z, 0) + 1246 return out247248def cmd_constants(args):249 c = log(4 / 3) / log(3)250 print("window edges", LOW, HIGH)251 print("first moment moves the window above alpha", up(1 - HIGH, 7))252 print("first moment closes the window at alpha", up(1 - LOW, 7))253 print("congruence decay cap c <=", up(c, 7))254 print("congruence alpha cap <=", up(1 / (2 - c), 6))255 print("trivial box C threshold reading", 1, "octave reading", 9, "delta 1 - 2 alpha")256 print("O needs theta <", down(1 / 0.5533, 4), "at alpha 0.5533, box theta 2")257258def cmd_check(args):259 print("k brute numpy CSlower ok")260 for k in range(2, 10):261 a = brute_ratios(k)262 b, second, _ = ratio_stats(k, 13, True)263 floor = (3 ** (k - 1) - 2 ** (k - 1)) ** 2 / second264 print(k, a, b, round(floor, 2), a == b and b >= floor)265 reg = census(9, 13, min(9, int(0.62 * 9) + 1), [], [7])266 print("regression A(9, 3^7)", reg[4][0], len(brute_rays(9, 3 ** 7)),267 reg[4][0] == 2818 == len(brute_rays(9, 3 ** 7)))268 print("n bruteA A bruteF F onecoord3 weight3 nonfibre")269 for n in range(4, 11):270 cutp = min(n, int(0.62 * n) + 1)271 a = 0.5272 rays = brute_rays(n, int(3.0 ** (a * n)))273 pts, fibre, weighted, heavy, occ, keys, kh = census(n, 13, cutp, [a])274 div3 = all((z[0] % 3 == 0) != (z[1] % 3 == 0) for z in rays)275 w3 = all((z[0] + z[1]) % 3 for z in rays)276 print(n, len(rays), occ[0], sum(rays.values()), int(heavy[0]), div3, w3,277 pts == 3 ** n - 2 ** (n + 1) + 1)278279def main():280 p = argparse.ArgumentParser()281 s = p.add_subparsers(dest="cmd", required=True)282 r = s.add_parser("ratios")283 r.add_argument("kmax", type=int)284 r.add_argument("--full", type=int, default=15)285 r.add_argument("--mem", type=int, default=13)286 r.set_defaults(fn=cmd_ratios)287 c = s.add_parser("census")288 c.add_argument("levels", type=int, nargs="+")289 c.add_argument("--cut", type=float, default=0.62)290 c.add_argument("--mem", type=int, default=13)291 c.add_argument("--alphas", type=float, nargs="+", default=[0.45, 0.5, 0.5533, 0.6])292 c.add_argument("--caps", type=int, nargs="*", default=[5, 6, 7])293 c.set_defaults(fn=cmd_census)294 q = s.add_parser("sieve")295 q.add_argument("levels", type=int, nargs="+")296 q.add_argument("--betas", type=float, nargs="+", default=[0.45, 0.5])297 q.add_argument("--mem", type=int, default=13)298 q.set_defaults(fn=cmd_sieve)299 n = s.add_parser("constants")300 n.set_defaults(fn=cmd_constants)301 k = s.add_parser("check")302 k.set_defaults(fn=cmd_check)303 a = p.parse_args()304 a.fn(a)305306main()