design_meter.py
22.2 kB · python · 613 lines
1import sys23import numpy as np4from mpmath import mp, findroot, mpc, mpf, power, zeta, zetazero5from scipy.ndimage import median_filter6from scipy.stats import spearmanr78SAMPLES = 327689BAND = (4.0, 60.0)10SCORE = 8.011FLOOR_WIDTH = 10112HEIGHT = 62.013ZEROS = 2014LINES = 1415SHIFTS = 400016SEED = 2026090717LADDER = (14, 16, 18, 20, 22)18LADDER_CENSUS = {14: (11, 105), 16: (149, 173), 18: (-30, 312), 20: (496, 539), 22: (1009, 1089)}19DRAWS = 2420SIGNS = 821DETAIL = ("3-01", "10-x9")22FIELD = ((100000000, 10, (0, 1, 2, 3, 4, 5, 6, 7, 8)), (129140163, 3, (0, 1)))23ESCAPE = ((10000000, 16, (0, 1)), (10000000, 10, (0, 1)))24ESCAPE_DECADES = (3, 4, 5, 6, 7)2526DESIGNS = {27 "3-01": (3, (0, 1), 20, 3),28 "3-02": (3, (0, 2), 20, 3),29 "3-12": (3, (1, 2), 18, 3),30 "4-01": (4, (0, 1), 18, 4),31 "4-012": (4, (0, 1, 2), 12, 4),32 "5-01": (5, (0, 1), 16, 5),33 "5-012": (5, (0, 1, 2), 12, 5),34 "9-012": (9, (0, 1, 2), 11, 3),35 "9-0123": (9, (0, 1, 2, 3), 10, 3),36 "9-0134": (9, (0, 1, 3, 4), 10, 3),37 "10-x9": (10, (9,), 100000000, 5),38}3940CENSUS = {"3-01": (496, 539), "3-02": (-382, 485), "3-12": (-1461, 1582), "10-x9": (2181, 5234)}4142QUADRATIC = {3: {1: 1, 2: -1}, 4: {1: 1, 3: -1}, 5: {1: 1, 2: -1, 3: -1, 4: 1}}4344CONTROL = 10000000045CONTROL_ANCHOR = 1928464748def primes_upto(limit):49 if limit < 2:50 return np.zeros(0, dtype=np.int64)51 sieve = np.ones(limit + 1, dtype=bool)52 sieve[:2] = False53 for p in range(2, int(limit**0.5) + 1):54 if sieve[p]:55 sieve[p * p :: p] = False56 return np.flatnonzero(sieve).astype(np.int64)575859def is_prime(n):60 for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):61 if n % p == 0:62 return n == p63 d = n - 164 r = 065 while d % 2 == 0:66 d //= 267 r += 168 for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):69 x = pow(a, d, n)70 if x == 1 or x == n - 1:71 continue72 for _ in range(r - 1):73 x = x * x % n74 if x == n - 1:75 break76 else:77 return False78 return True798081def mobius_values(values):82 mu = np.ones(values.size, dtype=np.int8)83 rem = values.astype(np.int64).copy()84 alive = np.flatnonzero(rem > 1)85 limit = int(round(float(values.max()) ** (1.0 / 3.0))) + 286 for p in primes_upto(limit):87 if alive.size == 0:88 break89 hit = alive[rem[alive] % p == 0]90 if hit.size:91 rem[hit] //= p92 mu[hit] = -mu[hit]93 again = hit[rem[hit] % p == 0]94 if again.size:95 mu[again] = 096 while again.size:97 rem[again] //= p98 again = again[rem[again] % p == 0]99 block = rem[alive]100 prime = alive[(block > 1) & (block < p * p)]101 mu[prime] = -mu[prime]102 alive = alive[rem[alive] >= p * p]103 for i in alive:104 c = int(rem[i])105 root = int(c**0.5)106 while root * root > c:107 root -= 1108 while (root + 1) * (root + 1) <= c:109 root += 1110 if root * root == c:111 mu[i] = 0112 elif is_prime(c):113 mu[i] = -mu[i]114 return mu115116117def elements(q, digits, length):118 lead = np.array([d for d in digits if d > 0], dtype=np.int64)119 tail = np.array(digits, dtype=np.int64)120 current = lead121 out = [lead]122 for _ in range(length - 1):123 current = (current[:, None] * q + tail[None, :]).ravel()124 out.append(current)125 values = np.concatenate(out)126 values.sort()127 return values128129130def design(q, digits, length):131 values = elements(q, digits, length)132 return values, mobius_values(values)133134135def mobius_sieve(n):136 rem = np.arange(n + 1, dtype=np.int32)137 mu = np.ones(n + 1, dtype=np.int8)138 mu[0] = 0139 for p in primes_upto(int(n**0.5) + 1):140 rem[p::p] //= p141 mu[p::p] = -mu[p::p]142 square = p * p143 if square <= n:144 mu[square::square] = 0145 power_of = square146 while power_of <= n:147 rem[power_of::power_of] //= p148 power_of *= p149 tail = rem > 1150 mu[tail] = -mu[tail]151 del rem152 return mu153154155def control_series(n, samples):156 mu = mobius_sieve(n)157 grid = np.exp(np.linspace(np.log(2.0), np.log(float(n)), samples))158 index = np.floor(grid).astype(np.int64)159 uniq, inverse = np.unique(index, return_inverse=True)160 starts = np.concatenate((np.ones(1, dtype=np.int64), uniq[:-1] + 1))161 running = np.cumsum(np.add.reduceat(mu, starts).astype(np.int64))162 del mu163 log_x = np.log(grid)164 return log_x, running[inverse] / np.exp(0.5 * log_x), int(running[-1])165166167def design_series(values, running, exponent, samples, window=None):168 lo = np.log(float(values[0])) if window is None else window[0]169 hi = np.log(float(values[-1])) if window is None else window[1]170 log_x = np.linspace(lo, hi, samples)171 slot = np.searchsorted(values, np.exp(log_x), side="right") - 1172 inside = slot >= 0173 out = np.zeros(samples, dtype=np.float64)174 out[inside] = running[slot[inside]]175 return log_x, out / np.exp(exponent * log_x)176177178def power_spectrum(log_x, series):179 step = (log_x[-1] - log_x[0]) / (log_x.size - 1)180 centred = series - series.mean()181 amplitude = np.fft.rfft(centred * np.hanning(log_x.size))182 gamma = np.fft.rfftfreq(log_x.size, d=step) * 2.0 * np.pi183 return gamma, np.abs(amplitude) ** 2184185186def scored(gamma, power):187 floor = median_filter(power, size=FLOOR_WIDTH, mode="nearest")188 return power / np.maximum(floor, np.finfo(float).tiny)189190191def window_max(score):192 return np.maximum.reduce([np.roll(score, 1), score, np.roll(score, -1)])193194195def band_slice(gamma):196 return np.flatnonzero((gamma > BAND[0]) & (gamma < BAND[1]))197198199def shift_test(gamma, score, targets, rng):200 inside = band_slice(gamma)201 lo, hi = float(gamma[inside[0]]), float(gamma[inside[-1]])202 span = hi - lo203 picked = np.array([t for t in targets if lo < t < hi], dtype=np.float64)204 if picked.size == 0:205 return 0, 0.0, 0.0, 1.0206 best = window_max(score)207208 def statistic(offset):209 moved = lo + np.mod(picked - lo + offset, span)210 slot = np.clip(np.searchsorted(gamma, moved), 1, gamma.size - 2)211 return float(np.mean(np.log10(best[slot])))212213 value = statistic(0.0)214 null = np.array([statistic(o) for o in rng.uniform(0.0, span, SHIFTS)])215 p = float((null >= value).mean())216 return picked.size, value, float(null.mean()), max(p, 1.0 / SHIFTS)217218219def peaks(gamma, score):220 inside = (gamma > BAND[0]) & (gamma < BAND[1])221 rise = np.zeros(score.size, dtype=bool)222 rise[1:-1] = (score[1:-1] > score[:-2]) & (score[1:-1] > score[2:])223 return np.flatnonzero(inside & rise & (score > SCORE))224225226def base_rate(gamma, score):227 inside = band_slice(gamma)228 return float((window_max(score)[inside] > SCORE).mean())229230231def zeta_zeros():232 mp.dps = 15233 out = []234 n = 1235 while True:236 g = float(zetazero(n).imag)237 if g > HEIGHT:238 return out239 out.append(g)240 n += 1241242243def dirichlet_l(s, q):244 table = QUADRATIC[q]245 return power(q, -s) * sum(c * zeta(s, mpf(a) / q) for a, c in table.items())246247248def l_zeros(q):249 mp.dps = 20250 grid = np.arange(0.5, HEIGHT, 0.05)251 mods = [abs(dirichlet_l(mpc(0.5, t), q)) for t in grid]252 out = []253 for i in range(1, grid.size - 1):254 if mods[i] < mods[i - 1] and mods[i] < mods[i + 1] and mods[i] < 0.3:255 root = findroot(lambda s: dirichlet_l(s, q), mpc(0.5, grid[i]))256 out.append(float(root.imag))257 return out258259260def lattice(q):261 return [2.0 * np.pi * j / np.log(q) for j in range(1, LINES + 1)]262263264def families(q, alpha, conductor, zeros, cache):265 if conductor not in cache:266 cache[conductor] = l_zeros(conductor)267 lz = cache[conductor]268 return (269 ("zeta", zeros),270 ("zeta+2", [g + 2.0 for g in zeros]),271 ("zetaflip", [64.0 - g for g in zeros]),272 ("alphazeta", [alpha * g for g in zeros]),273 (f"L{conductor}", lz),274 ("lattice", lattice(q)),275 )276277278def split_correlation(values, running, exponent):279 lo = np.log(float(values[0]))280 hi = np.log(float(values[-1]))281 mid = 0.5 * (lo + hi)282 halves = []283 for window in ((lo, mid), (mid, hi)):284 log_x, series = design_series(values, running, exponent, SAMPLES // 2, window)285 gamma, power = power_spectrum(log_x, series)286 halves.append((gamma, scored(gamma, power)))287 inside = band_slice(halves[0][0])288 rho, _ = spearmanr(halves[0][1][inside], halves[1][1][inside])289 return float(rho), float(halves[0][0][1])290291292def split_null(values, support, exponent, rng, draws):293 out = []294 for _ in range(draws):295 signs = np.zeros(values.size, dtype=np.int32)296 signs[support] = rng.choice([-1, 1], int(support.sum()))297 out.append(split_correlation(values, np.cumsum(signs, dtype=np.int64), exponent)[0])298 return float(np.mean(out)), float(np.std(out))299300301def digit_mask(n, base, missing):302 mask = np.ones(n + 1, dtype=bool)303 place = 1304 while place <= n:305 block = base * place306 size = ((n + 1) // block) * block307 view = mask[:size].reshape(-1, base, place)308 view[:, missing, :] = False309 tail = mask[size:]310 if tail.size:311 offset = missing * place312 if offset < tail.size:313 tail[offset : offset + place] = False314 place = block315 mask[0] = False316 return mask317318319def digit_set_mask(n, base, digits):320 mask = np.ones(n + 1, dtype=bool)321 for d in range(base):322 if d not in digits:323 mask &= digit_mask(n, base, d)324 mask[0] = False325 return mask326327328def grid_running(weights, grid):329 index = np.floor(grid).astype(np.int64)330 uniq, inverse = np.unique(index, return_inverse=True)331 starts = np.concatenate((np.ones(1, dtype=np.int64), uniq[:-1] + 1))332 out = []333 for w in weights:334 out.append(np.cumsum(np.add.reduceat(w, starts).astype(np.int64))[inverse])335 return out336337338def mean_field(n, base, digits, samples):339 alpha = np.log(len(digits)) / np.log(base)340 mu = mobius_sieve(n)341 mask = digit_set_mask(n, base, digits)342 grid = np.exp(np.linspace(np.log(2.0), np.log(float(n)), samples))343 log_x = np.log(grid)344 uniq, inverse = np.unique(np.floor(grid).astype(np.int64), return_inverse=True)345 bounds = np.concatenate((np.ones(1, dtype=np.int64), uniq + 1))346 echo = np.empty(uniq.size)347 total = 0.0348 seen = 0349 for j in range(uniq.size):350 lo, hi = bounds[j], bounds[j + 1]351 run = seen + np.cumsum(mask[lo:hi], dtype=np.int64)352 line = np.arange(lo, hi, dtype=np.float64)353 total += float(np.dot(mu[lo:hi].astype(np.float64), run / line))354 seen = int(run[-1]) if run.size else seen355 echo[j] = total356 meter = np.cumsum(np.add.reduceat(mu * mask, bounds[:-1]))357 del mu, mask358 scale = np.exp(alpha * 0.5 * log_x)359 return log_x, alpha, meter[inverse] / scale, echo[inverse] / scale, (meter[inverse] - echo[inverse]) / scale, int(meter[-1]), seen360361362def upper_rms(series):363 half = series[series.size // 2 :]364 return float(np.sqrt(np.mean(half * half)))365366367def kempner(n, base, missing):368 mu = mobius_sieve(n)369 mask = digit_mask(n, base, missing)370 values = np.flatnonzero(mask).astype(np.int64)371 weights = mu[values].copy()372 del mu, mask373 return values, weights374375376def line(tag, name, count, value, null, p):377 print(f" {tag:8s} {name:8s} n {count:3d} logscore {value:6.3f} null {null:6.3f} p {p:.4f}")378379380def verb_sieve():381 for name in DESIGNS:382 q, alpha, _, values, mu = sources(name)383 running = np.cumsum(mu, dtype=np.int64)384 top = int(running[-1])385 peak = int(np.abs(running).max())386 span = np.log(float(values[-1])) - np.log(float(values[0]))387 print(388 f"design {name:8s} q {q:2d} alpha {alpha:.6f} A {values.size:8d}"389 f" M {top:6d} Mmax {peak:6d} thetamax {np.log(peak) / np.log(values.size):.4f}"390 f" support {int((mu != 0).sum()):8d} logspan {span:.4f} bin {2 * np.pi / span:.4f}"391 )392 want = CENSUS.get(name)393 if want:394 print(f"anchor {name:8s} census {want} meter {(top, peak)} match {want == (top, peak)}")395 del values, mu, running396 lhs, lmu = design(9, (0, 1, 3, 4), 10)397 rhs, rmu = design(3, (0, 1), 20)398 same = bool(np.array_equal(lhs, rhs)) and bool(np.array_equal(lmu, rmu))399 print(f"identity base 9 F 0134 equals base 3 F 01 at 3^20 {same}")400 _, _, anchor = control_series(CONTROL, 1024)401 print(f"control base 10 full set M(10^8) {anchor} A084237 {CONTROL_ANCHOR} match {anchor == CONTROL_ANCHOR}")402403404def sources(name):405 q, digits, deep, conductor = DESIGNS[name]406 if q == 10:407 values, mu = kempner(deep, 10, digits[0])408 alpha = np.log(9.0) / np.log(10.0)409 else:410 values, mu = design(q, digits, deep)411 alpha = np.log(len(digits)) / np.log(q)412 return q, alpha, conductor, values, mu413414415def census_row(name, rng, zeros, cache, detail=False):416 q, alpha, conductor, values, mu = sources(name)417 support = mu != 0418 running = np.cumsum(mu, dtype=np.int64)419 top = int(running[-1])420 peak = int(np.abs(running).max())421 log_x, series = design_series(values, running, alpha * 0.5, SAMPLES)422 gamma, power = power_spectrum(log_x, series)423 score = scored(gamma, power)424 rho, _ = split_correlation(values, running, alpha * 0.5)425 mean, sigma = split_null(values, support, alpha * 0.5, rng, DRAWS if values.size < 5000000 else 6)426 z = (rho - mean) / sigma if sigma > 0 else 0.0427 cells = []428 for label, targets in families(q, alpha, conductor, zeros, cache):429 _, value, null, p = shift_test(gamma, score, targets, rng)430 cells.append(f"{label} {value:6.3f}/{null:5.3f} p {p:.4f}")431 print(432 f"census {name:8s} q {q:2d} alpha {alpha:.6f} A {values.size:8d}"433 f" M {top:6d} Mmax {peak:6d} bin {gamma[1]:.4f} rate {base_rate(gamma, score):.4f}"434 f" split {rho:7.4f} z {z:5.2f}"435 )436 print(" meter " + " ".join(cells))437 if detail:438 for tag, weight in (("count", np.ones(values.size, dtype=np.int64)), ("random", None)):439 if weight is None:440 weight = np.zeros(values.size, dtype=np.int32)441 weight[support] = rng.choice([-1, 1], int(support.sum()))442 run = np.cumsum(weight)443 exponent = alpha if tag == "count" else alpha * 0.5444 lx, sr = design_series(values, run, exponent, SAMPLES)445 g, pw = power_spectrum(lx, sr)446 sc = scored(g, pw)447 cells = []448 for label, targets in families(q, alpha, conductor, zeros, cache):449 _, value, null, p = shift_test(g, sc, targets, rng)450 cells.append(f"{label} {value:6.3f}/{null:5.3f} p {p:.4f}")451 print(f" {tag:8s} " + " ".join(cells))452 print(f" {tag:8s} split {split_correlation(values, run, exponent)[0]:7.4f}")453 found = peaks(gamma, score)454 for i in sorted(found[np.argsort(-score[found])][:6]):455 row = f" peak {gamma[i]:8.4f} score {score[i]:8.2f}"456 for label, targets in families(q, alpha, conductor, zeros, cache):457 near = min(targets, key=lambda z: abs(z - gamma[i]))458 row += f" {label} {near:7.3f} off {abs(near - gamma[i]):5.3f}"459 print(row)460 return top, peak461462463def escape_line(n, base, digits):464 log_x, alpha, _, echo, _, top, count = mean_field(n, base, digits, SAMPLES)465 raw = echo * np.exp(alpha * 0.5 * log_x)466 print(467 f"escape base {base} F {''.join(str(d) for d in digits)} N {n} alpha {alpha:.6f}"468 f" A_F {count} M_F {top} echo rms {upper_rms(echo):.6f}"469 )470 for decade in ESCAPE_DECADES:471 j = min(int(np.searchsorted(log_x, decade * np.log(10.0))), log_x.size - 1)472 old_bound = float(np.exp((alpha - 0.5) * log_x[j]))473 print(474 f" echo x 1e{decade} sum mu A_F over n {float(raw[j]):8.4f}"475 f" old bound x^(alpha-1/2) {old_bound:7.4f} ratio {abs(float(raw[j])) / old_bound:6.2f}"476 )477 print()478479480def verb_spectrum():481 rng = np.random.default_rng(SEED)482 zeros = zeta_zeros()483 cache = {}484 log_x, series, anchor = control_series(CONTROL, SAMPLES)485 gamma, power = power_spectrum(log_x, series)486 score = scored(gamma, power)487 print(488 f"control base 10 full set N {CONTROL} M {anchor} A084237 {CONTROL_ANCHOR}"489 f" match {anchor == CONTROL_ANCHOR} logspan {log_x[-1] - log_x[0]:.4f}"490 f" bin {gamma[1]:.4f} rate {base_rate(gamma, score):.4f}"491 )492 line("control", "zeta", *shift_test(gamma, score, zeros, rng))493 found = peaks(gamma, score)494 order = found[np.argsort(-score[found])][:10]495 matched = 0496 for i in sorted(order):497 near = min(zeros, key=lambda z: abs(z - gamma[i]))498 matched += abs(near - gamma[i]) < gamma[1]499 print(f" peak {gamma[i]:8.4f} score {score[i]:11.1f} zeta {near:8.4f} off {abs(near - gamma[i]):5.3f}")500 print(f" control top 10 peaks within one bin of a zeta zero {matched} of 10")501 print()502 for name in DESIGNS:503 got = census_row(name, rng, zeros, cache, detail=name in DETAIL)504 want = CENSUS.get(name)505 if want:506 print(f"anchor {name:8s} census {want} meter {got} match {want == got}")507 print()508 for n, base, digits in FIELD:509 log_x, alpha, meter, echo, residual, top, count = mean_field(n, base, digits, SAMPLES)510 print(511 f"field base {base} F {''.join(str(d) for d in digits)} N {n} alpha {alpha:.6f}"512 f" A_F {count} M_F {top} bin {2 * np.pi / (log_x[-1] - log_x[0]):.4f}"513 )514 for tag, series in (("meter", meter), ("echo", echo), ("residual", residual)):515 g, pw = power_spectrum(log_x, series)516 sc = scored(g, pw)517 _, value, null, p = shift_test(g, sc, zeros, rng)518 top6 = peaks(g, sc)519 top6 = top6[np.argsort(-sc[top6])][:6]520 near = sum(1 for i in top6 if min(abs(z - g[i]) for z in zeros) < g[1])521 print(522 f" {tag:9s} zeta {value:6.3f}/{null:5.3f} p {p:.4f}"523 f" rms {upper_rms(series):8.4f} top6 within one bin {near} of {top6.size}"524 )525 cut = log_x.size526 base_ratio = None527 for fraction in (0.6, 0.8, 1.0):528 take = int(cut * fraction)529 ratio = upper_rms(echo[:take]) / upper_rms(meter[:take])530 span = log_x[take - 1]531 if base_ratio is None:532 base_ratio, base_span = ratio, span533 rate = min(alpha, 1.0 - alpha)534 predicted = base_ratio * np.exp(-0.5 * rate * (span - base_span))535 print(536 f" decay logx {span:7.4f} echo over meter {ratio:.6f}"537 f" predicted {predicted:.6f} ratio {ratio / predicted:.4f}"538 )539 print()540 for n, base, digits in ESCAPE:541 escape_line(n, base, digits)542 for deep in LADDER:543 values, mu = design(3, (0, 1), deep)544 alpha = np.log(2.0) / np.log(3.0)545 running = np.cumsum(mu, dtype=np.int64)546 log_x, series = design_series(values, running, alpha * 0.5, SAMPLES)547 gamma, power = power_spectrum(log_x, series)548 score = scored(gamma, power)549 _, value, null, p = shift_test(gamma, score, zeros, rng)550 rho, _ = split_correlation(values, running, alpha * 0.5)551 mean, sigma = split_null(values, mu != 0, alpha * 0.5, rng, DRAWS)552 got = (int(running[-1]), int(np.abs(running).max()))553 want = LADDER_CENSUS.get(deep)554 print(555 f"ladder 3-01 L {deep:2d} A {values.size:8d} bin {gamma[1]:.4f}"556 f" zeta {value:6.3f}/{null:5.3f} p {p:.4f}"557 f" split {rho:7.4f} z {(rho - mean) / sigma if sigma > 0 else 0.0:5.2f}"558 f" census {want} meter {got} match {want == got}"559 )560561562def verb_family():563 rng = np.random.default_rng(SEED)564 zeros = zeta_zeros()565 cache = {}566 for left, right in (("4-01", "9-012"), ("3-01", "9-0123"), ("3-01", "9-0134"), ("3-01", "3-02"), ("10-x9", "3-01")):567 rows = []568 for name in (left, right):569 q, alpha, conductor, values, mu = sources(name)570 log_x, series = design_series(values, np.cumsum(mu, dtype=np.int64), alpha * 0.5, SAMPLES)571 gamma, power = power_spectrum(log_x, series)572 score = scored(gamma, power)573 rows.append((name, q, alpha, gamma, score, peaks(gamma, score), values, mu))574 (an, aq, aa, ag, asc, ap, av, amu), (bn, bq, ba, bg, bsc, bp, bv, bmu) = rows575 width = max(float(ag[1]), float(bg[1]))576 shared = sum(1 for i in ap if bp.size and np.min(np.abs(bg[bp] - ag[i])) < width)577 cover = min(1.0, bp.size * 2.0 * width / (BAND[1] - BAND[0]))578 print(579 f"family {an:8s} alpha {aa:.6f} peaks {ap.size} against {bn:8s} alpha {ba:.6f}"580 f" peaks {bp.size} shared {shared} expect {cover * ap.size:.2f} bin {width:.4f}"581 )582 grid = np.linspace(BAND[0], BAND[1], 4096)583 rho, pv = spearmanr(np.interp(grid, ag, asc), np.interp(grid, bg, bsc))584 print(f" spectra rho {float(rho):7.4f} p {float(pv):.3e}")585 for name, alpha, values, mu in ((an, aa, av, amu), (bn, ba, bv, bmu)):586 if values.size >= 5000000:587 continue588 support = mu != 0589 counts = []590 for _ in range(SIGNS):591 weight = np.zeros(values.size, dtype=np.int32)592 weight[support] = rng.choice([-1, 1], int(support.sum()))593 lx, sr = design_series(values, np.cumsum(weight, dtype=np.int64), alpha * 0.5, SAMPLES)594 g, pw = power_spectrum(lx, sr)595 counts.append(int(peaks(g, scored(g, pw)).size))596 print(f" random {name:8s} peaks {min(counts)} to {max(counts)} over {SIGNS} sign draws")597 for name, q, alpha, gamma, score in ((an, aq, aa, ag, asc), (bn, bq, ba, bg, bsc)):598 conductor = DESIGNS[name][3]599 cells = []600 for label, targets in families(q, alpha, conductor, zeros, cache):601 _, value, null, p = shift_test(gamma, score, targets, rng)602 cells.append(f"{label} {value:6.3f}/{null:5.3f} p {p:.4f}")603 print(f" {name:8s} " + " ".join(cells))604 print()605606607def main():608 verb = sys.argv[1] if len(sys.argv) > 1 else "sieve"609 {"sieve": verb_sieve, "spectrum": verb_spectrum, "family": verb_family}[verb]()610611612if __name__ == "__main__":613 main()