#!/usr/bin/env python3 """Smart Mathematical Analysis & Reverse-Engineering for Bitcoin Challenge Puzzles. Reads the 160-puzzle Bitcoin Challenge dataset (BTC-Solved-Unsolved.txt if present, otherwise btc-puzzles.json) and runs the four forensic steps on the SOLVED private keys: 1. Mathematical gap & ratio profiling (offsets, linear, recurrence fits) 2. Bit-level & masking distribution (density, popcount, nibble masks) 3. Predictive formula derivation (OLS fits + out-of-sample check) 4. Actionable lightweight strategy (what, if anything, shrinks a range) Optionally predicts the honest sub-space of an unsolved target: python3 math_check.py # full report, target #71 python3 math_check.py --predict 121 # any unsolved puzzle python3 math_check.py --data my.txt # parse the raw challenge dump After the four steps it runs five advanced non-linear / structural checks on the solved keys only (bit-position bias matrix, spectral & multi-lag autocorrelation, cross-key algebraic relations, nibble/byte entropy, Benford's Law), with family-wise alpha = 0.01 per check. It then runs five deep cryptographic / dynamic / information-theoretic checks (pairwise Hamming matrix, NIST runs + LZ78 complexity, bit-segment mutual information, Collatz / Rule-30 dynamics, secp256k1 field residue), also with family-wise alpha = 0.01 per check. Pure standard library, no GPU, no network, no dependencies. The first four steps complete in milliseconds; the deep section (checks 11-15) additionally simulates its Monte-Carlo nulls with a fixed xorshift32 stream and enumerates the exact LZ78 null for bit lengths up to 20, so the whole report finishes in about ten seconds. The report only prints numbers measured from the dataset itself. If the data contains no deterministic structure the script says so with statistics instead of inventing a formula. """ import argparse import json import math import os import re import sys # -------------------------------------------------------------------------- # dataset loading # -------------------------------------------------------------------------- ROW_RE = re.compile( r"\|\s*([0-9a-fA-F]+)\s*:\s*([0-9a-fA-F]+)\s*\|\s*(\S+)\s*\|\s*" r"(SOLVED|UNSOLVED)\s*\|([^|]*)\|([^|]*)\|" ) def load_from_text(path): """Parse the raw `| lo:hi | address | STATUS | pubkey | privkey |` dump.""" rows = [] with open(path, "r", encoding="utf-8", errors="replace") as fh: for line in fh: m = ROW_RE.search(line) if not m: continue lo_hex, hi_hex, address, status, _pub, key = ( m.group(1), m.group(2), m.group(3), m.group(4), m.group(5), m.group(6), ) lo, hi = int(lo_hex, 16), int(hi_hex, 16) # puzzle number is the bit length of the range: 2^(n-1) <= k < 2^n n = hi.bit_length() entry = {"n": n, "solved": status == "SOLVED", "lo": lo, "hi": hi, "address": address.strip()} if entry["solved"]: key = key.strip() if key and key.upper() != "NOT SOLVED": entry["key"] = int(key, 16) rows.append(entry) rows.sort(key=lambda r: r["n"]) return rows def load_from_json(path): with open(path, "r", encoding="utf-8") as fh: doc = json.load(fh) rows = [] for p in doc["puzzles"]: n = p["n"] lo, hi = 1 << (n - 1), (1 << n) - 1 entry = {"n": n, "solved": bool(p.get("solved")), "lo": lo, "hi": hi, "address": p.get("address", "")} if entry["solved"]: entry["key"] = int(p["key"], 16) rows.append(entry) rows.sort(key=lambda r: r["n"]) return rows def load_dataset(explicit=None): here = os.path.dirname(os.path.abspath(__file__)) if explicit: path = explicit else: path = None for name in ("BTC-Solved-Unsolved.txt", "btc-puzzles.json"): cand = os.path.join(here, name) if os.path.exists(cand): path = cand break if os.path.exists(name): path = name break if path is None: sys.exit("error: no dataset found (BTC-Solved-Unsolved.txt or " "btc-puzzles.json next to this script, or use --data)") return (load_from_text(path) if path.endswith((".txt", ".csv")) else load_from_json(path)), path # -------------------------------------------------------------------------- # small statistics helpers (stdlib only) # -------------------------------------------------------------------------- def nsum(xs): """Plain sequential IEEE-754 summation. CPython's built-in sum() uses compensated summation (Python 3.12+), while the math.html page adds the same terms with a naive loop. Summing the same way here keeps this script's output bit-identical to the page instead of differing in the final double (the compensation shifts sxx by one ulp, which changes the fitted median's low hex digits). """ s = 0.0 for x in xs: s += x return s def mean(xs): return nsum(xs) / len(xs) def sd(xs): m = mean(xs) return math.sqrt(nsum((x - m) * (x - m) for x in xs) / len(xs)) def pearson(xs, ys): mx, my = mean(xs), mean(ys) num = nsum((x - mx) * (y - my) for x, y in zip(xs, ys)) den = math.sqrt(nsum((x - mx) * (x - mx) for x in xs) * nsum((y - my) * (y - my) for y in ys)) return num / den if den else 0.0 def ols(xs, ys): """Ordinary least squares: returns (slope, intercept, r_squared).""" mx, my = mean(xs), mean(ys) sxx = nsum((x - mx) * (x - mx) for x in xs) sxy = nsum((x - mx) * (y - my) for x, y in zip(xs, ys)) slope = sxy / sxx intercept = my - slope * mx syy = nsum((y - my) * (y - my) for y in ys) r2 = (sxy * sxy) / (sxx * syy) if sxx and syy else 0.0 return slope, intercept, r2 def popcount(x): return bin(x).count("1") # -------------------------------------------------------------------------- # the four forensic steps # -------------------------------------------------------------------------- def analyze(rows, out=print): solved = [r for r in rows if r["solved"] and "key" in r] n_rows = len(rows) n_solved = len(solved) out("=" * 74) out("SMART MATHEMATICAL ANALYSIS — BITCOIN CHALLENGE PUZZLES") out("=" * 74) out(f"puzzles loaded : {n_rows}") out(f"solved with key: {n_solved} unsolved: {n_rows - n_solved}") in_range = all(r["lo"] <= r["key"] <= r["hi"] for r in solved) bits_ok = all(r["key"].bit_length() == r["n"] for r in solved) out(f"range check : {in_range} (every key inside its own 2^(n-1)..2^n-1)") out(f"bit check : {bits_ok} (every key exactly n significant bits)") # ---- step 1: gap & ratio profiling ------------------------------- out("") out("-" * 74) out("STEP 1 — MATHEMATICAL GAP & RATIO PROFILING") out("-" * 74) gaps = [solved[i]["key"] - solved[i - 1]["key"] for i in range(1, n_solved)] # round to double first — the same path as the page's Number(key) division ratios = [float(solved[i]["key"]) / float(solved[i - 1]["key"]) for i in range(1, n_solved)] gapf = [float(g) for g in gaps] # same rounding as the page's Number(gap) cv = sd(gapf) / mean(gapf) out(f"consecutive pairs : {len(gaps)}") out(f"raw gap k[i+1]-k[i] : {gaps[0]} .. {gaps[-1]} " f"(grows with the ranges themselves)") out(f"gap coefficient of variation: {cv:.3f} " f"(0.000 = constant offset, ~8 = pure scale growth)") out(f"ratio k[i+1]/k[i] : mean {mean(ratios):.4f} " f"sd {sd(ratios):.4f} min {min(ratios):.3f} max {max(ratios):.3f}") out(f"constant offset test : {'REJECTED' if cv > 0.5 else 'not rejected'}" f" (CV={cv:.3f}; a creator offset would be 0)") best_a, best_hits = 0, 0 for a in range(1, 9): hits = sum(1 for i in range(1, n_solved) if abs(solved[i]["key"] - solved[i - 1]["key"] * a) * 10 <= solved[i]["key"]) if hits > best_hits: best_a, best_hits = a, hits pct = 100.0 * best_hits / (n_solved - 1) out(f"multiplicative recurrence : best k[i+1]=a*k[i] is a={best_a}, " f"explains {best_hits}/{n_solved - 1} pairs ({pct:.1f}%)" f" <- a<=100% is no relation") # double-rounding path identical to the page's Number(k-lo) / Number(lo) pos = [float(r["key"] - r["lo"]) / float(r["lo"]) for r in solved] ns = [r["n"] for r in solved] out("") out("normalised position inside each range: p = (k - 2^(n-1)) / 2^(n-1)") out(f" mean {mean(pos):.4f} sd {sd(pos):.4f} " f"(uniform on [0,1) expects 0.5000 / 0.2887)") out(f" min {min(pos):.4f} (#{ns[pos.index(min(pos))]}) " f"max {max(pos):.4f} (#{ns[pos.index(max(pos))]})") r_pos = pearson(ns, pos) crit = 1.96 / math.sqrt(len(pos)) out(f" Pearson r(position, puzzle#) = {r_pos:+.4f} " f"(95% critical value {crit:.3f}) -> " f"{'NO significant trend' if abs(r_pos) < crit else 'significant trend'}") bins = [0] * 10 for p in pos: bins[min(9, int(p * 10))] += 1 exp = len(pos) / 10 chi2 = nsum((o - exp) * (o - exp) / exp for o in bins) verdict = "consistent with uniform" if chi2 < 16.92 else "STRUCTURE PRESENT" out(f" decile counts : {bins} (expect {exp:.1f} each)") out(f" chi-square(9 df) = {chi2:.2f} vs 16.92 (5%) / 21.67 (1%) -> {verdict}") lag1 = nsum((pos[i] - mean(pos)) * (pos[i - 1] - mean(pos)) for i in range(1, len(pos))) lag1 /= (len(pos) - 1) * (nsum((y - mean(pos)) * (y - mean(pos)) for y in pos) / len(pos)) out(f" lag-1 autocorrelation of p = {lag1:+.3f} (white noise ~0)") # ---- step 2: bit level & masking --------------------------------- out("") out("-" * 74) out("STEP 2 — BIT-LEVEL & MASKING DISTRIBUTION") out("-" * 74) pr = [popcount(r["key"]) / r["n"] for r in solved] odds = sum(1 for r in solved if r["key"] & 1) out(f"popcount(key)/n : mean {mean(pr):.4f} sd {sd(pr):.4f}" f" (fair coin expects 0.5000)") out(f"odd keys : {odds}/{n_solved} " f"(fair expectation {n_solved / 2:.1f})") zeros = total = 0 for r in solved: hexlen = (r["n"] + 3) // 4 for ch in format(r["key"], f"0{hexlen}x"): total += 1 zeros += ch == "0" out(f"zero nibbles inside keys : {zeros}/{total} = {zeros / total:.4f}" f" (uniform hex expects 0.0625)") mask_ok = True for r in solved: first = format(r["key"], "x")[0] allowed = {1: "123456789abcdef", 2: "234567", 3: "4567", 0: "89abcdef"}[r["n"] % 4] if first not in allowed: mask_ok = False out(f"first hex digit inside the mask the RANGE forces : {mask_ok}") out(" (the mask is the puzzle range itself: 2^(n-1) <= k < 2^n forces the") out(" leading bit and the leading nibble — it is not a creator pattern)") dens = [0.0] * 4 cnts = [0] * 4 for r in solved: hexlen = (r["n"] + 3) // 4 for i, ch in enumerate(format(r["key"], f"0{hexlen}x")): d = int(ch, 16) slot = i % 4 for b in range(4): dens[slot] += (d >> (3 - b)) & 1 cnts[slot] += 1 out("bit density by nibble slot : " + " ".join( f"slot{i}={dens[i] / cnts[i]:.3f}" for i in range(4)) + " (random expects 0.500)") # ---- step 3: predictive formula ---------------------------------- out("") out("-" * 74) out("STEP 3 — PREDICTIVE FORMULA DERIVATION") out("-" * 74) logs = [math.log2(float(r["key"])) for r in solved] slope, intercept, r2 = ols(ns, logs) out(f"fit A: log2(k) = {slope:.6f} * n + {intercept:.6f} R^2 = {r2:.6f}") out(" (slope 1 / intercept 0 is forced by the range definition —") out(" this fit describes the boundaries, not the creator)") resid = [logs[i] - (slope * ns[i] + intercept) for i in range(n_solved)] out(f" residual sd = {sd(resid):.4f} bits") slope_p, intercept_p, r2_p = ols(ns, pos) sse = 0.0 for i in range(n_solved): e = pos[i] - (slope_p * ns[i] + intercept_p) sse += e * e resid_p = math.sqrt(sse / max(1, n_solved - 2)) out(f"fit B: position p = {slope_p:.6f} * n + {intercept_p:.6f} " f"R^2 = {r2_p:.4f} residual sd {resid_p:.4f}") band = min(100.0, 3.92 * resid_p * 100.0) out(f" 95% band (+-1.96 sd) covers {band:.0f}% of the half-range " f"-> search-space reduction: {'NONE' if band >= 100 else 'partial'}") train = [r for r in solved if r["n"] <= 130] if any(r["n"] == 135 for r in solved) and len(train) >= 3: tpos = [float(r["key"] - r["lo"]) / float(r["lo"]) for r in train] tns = [r["n"] for r in train] ts, ti, _ = ols(tns, tpos) p135 = next(r for r in solved if r["n"] == 135) actual = float(p135["key"] - p135["lo"]) / float(p135["lo"]) predicted = ts * 135 + ti out(f"out-of-sample check : fit #1-#{train[-1]['n']}, predict #135 " f"-> p={predicted:.3f}, actual {actual:.3f}") out(f" : error {abs(predicted - actual):.3f} " f"vs residual sd {sd(tpos):.3f} -> prediction misses by " f"{abs(predicted - actual) / sd(tpos):.1f} sd") out("") out("VERDICT: no deterministic key-generation formula is present in this") out("dataset. Positions are uniform (chi-square not rejected), trends are") out("below the significance threshold, and the fitted band covers the whole") out("range. Any 'formula' would be fitted noise, not a creator's algorithm.") return solved, slope_p, intercept_p, resid_p # -------------------------------------------------------------------------- # step 4: actionable strategy for an unsolved target # -------------------------------------------------------------------------- def predict(solved_rows, target_n, slope_p, intercept_p, resid_p, out=print): lo, hi = 1 << (target_n - 1), (1 << target_n) - 1 span = hi - lo + 1 out("") out("-" * 74) out(f"STEP 4 — ACTIONABLE STRATEGY: PUZZLE #{target_n}") out("-" * 74) out(f"range (exact) : 2^({target_n - 1}) .. 2^{target_n} - 1") out(f" {format(lo, 'x')} .. {format(hi, 'x')}") out(f"keys in range : {span} (2^{target_n - 1} = {span:.3e})") p_hat = min(1.0, max(0.0, slope_p * target_n + intercept_p)) median_key = lo + int(p_hat * lo) band_pct = min(100.0, 3.92 * resid_p * 100.0) band_keys = int(span * min(1.0, 3.92 * resid_p)) out(f"model median guess : {format(median_key, 'x')} (position {p_hat:.3f})") out(f"model 95% band : {band_pct:.0f}% of the range " f"({band_keys:.3e} keys)") out("") if band_pct >= 100: out("HONEST SUB-SPACE: the full range. The model's own uncertainty") out("already spans all of 2^(n-1), so it removes exactly zero keys:") out("search-space reduction = 0. No sub-space claim is supported by") out("this data, and none is asserted here.") else: out(f"HONEST SUB-SPACE: the {band_pct:.0f}% band — but see the") out("out-of-sample error above before trusting it on a new puzzle.") out("") out("zero-GPU actions that are actually available:") out(f" 1. small targets (n<=~35) are minutes-to-hours on one CPU core;") out(" 2. everything above that is range-limited, not math-limited: no") out(" analysis of these 83 keys shrinks 2^(n-1) by a single bit;") out(" 3. run this script again when new puzzles are solved — a genuine") out(" creator formula would collapse the residual sd toward 0, and") out(" the tests above would reject 'uniform' instead of failing to.") # -------------------------------------------------------------------------- # shared special functions (ported verbatim to math-adv.js) # -------------------------------------------------------------------------- # Bonferroni two-sided z thresholds, family-wise alpha = 0.01: # inv_normal(1 - 0.01 / (2 * comparisons)); the same literals appear in # math-adv.js so the page and this script print identical numbers. Z_A1_REL = 3.7809198112947677 # 64 relative bit buckets Z_A1_ABS = 4.003167571450172 # 160 absolute bit positions Z_A2_ACF = 3.764823649533927 # 60 value-series lags Z_A2_STR = 4.42229885489002 # 1024 bit-stream lags _LG = (0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7) def lgamma(x): """Lanczos g=7, n=9 — identical coefficients and order in math-adv.js.""" if x < 0.5: return math.log(math.pi / math.sin(math.pi * x)) - lgamma(1.0 - x) z = x - 1.0 a = _LG[0] t = z + 7.5 for i in range(1, 9): a += _LG[i] / (z + i) return 0.5 * math.log(2.0 * math.pi) + (z + 0.5) * math.log(t) - t + math.log(a) def gammainc_upper(a, x): """Regularized upper incomplete gamma Q(a, x) (series / continued fraction).""" if x <= 0.0: return 1.0 if x < a + 1.0: ap, s, d = a, 1.0 / a, 1.0 / a for _ in range(500): ap += 1.0 d *= x / ap s += d if abs(d) < abs(s) * 1e-16: break return 1.0 - s * math.exp(-x + a * math.log(x) - lgamma(a)) tiny = 1e-300 b = x + 1.0 - a c = 1.0 / tiny d = 1.0 / b h = d for i in range(1, 500): an = -i * (i - a) b += 2.0 d = an * d + b if abs(d) < tiny: d = tiny c = b + an / c if abs(c) < tiny: c = tiny d = 1.0 / d de = d * c h *= de if abs(de - 1.0) < 1e-16: break return math.exp(-x + a * math.log(x) - lgamma(a)) * h def chi2sf(x, k): """Upper tail of the chi-square distribution with k degrees of freedom.""" return gammainc_upper(k / 2.0, x / 2.0) def poisson_tail(k, lam): """P(X >= k) for X ~ Poisson(lam) — null probability of >= k hits.""" if k <= 0: return 1.0 if lam <= 0.0: return 0.0 term = math.exp(-lam) s = term for i in range(1, k): term *= lam / i s += term return max(0.0, min(1.0, 1.0 - s)) def hex_entropy(cnt, total): """Plug-in Shannon entropy (bits) + observed category count.""" h, m = 0.0, 0 for c in cnt: if c: p = c / total h -= p * math.log2(p) m += 1 return h, m # -------------------------------------------------------------------------- # five advanced non-linear & structural checks (solved keys only) # -------------------------------------------------------------------------- def advanced(solved, out=print): """Checks 1-5 of the advanced forensic protocol; returns measured results. Everything printed here is computed from the solved keys themselves. Each check carries a family-wise alpha of 0.01 (Bonferroni where it scans many candidates); the verdict line states whether any of them flags. """ N = len(solved) ns = [r["n"] for r in solved] keys = [r["key"] for r in solved] hexlens = [(n + 3) // 4 for n in ns] maxn = max(ns) flags = [] res = {} # ---- check 1: bit-position bias matrix --------------------------- out("") out("-" * 74) out("ADVANCED CHECK 1 — BIT-POSITION BIAS MATRIX (positional probability)") out("-" * 74) BINS = 64 ones = [0] * BINS tot = [0] * BINS msb_pinned = 0 total_bits = 0 for r in solved: n, k = r["n"], r["key"] total_bits += n for p in range(n): b = (k >> p) & 1 rel = 1.0 if n == 1 else p / (n - 1) bi = min(BINS - 1, int(rel * BINS)) if p == n - 1: # the leading bit: forced to 1 by 2^(n-1) <= k msb_pinned += b else: tot[bi] += 1 ones[bi] += b chi = 0.0 nb = 0 zmax, zmax_b = 0.0, -1 zflags = 0 zs = [None] * BINS for bi in range(BINS): if not tot[bi]: continue nb += 1 e = 0.5 * tot[bi] o = ones[bi] chi += (o - e) * (o - e) / e z = (o / tot[bi] - 0.5) / math.sqrt(0.25 / tot[bi]) zs[bi] = z if abs(z) > abs(zmax): zmax, zmax_b = z, bi if abs(z) > Z_A1_REL: zflags += 1 chi_p = chi2sf(chi, nb - 1) out(f"relative buckets (r = p/(n-1)), MSB excluded : {nb} buckets, " f"{sum(tot)} bits ({total_bits} total - {N} pinned MSBs)") out(f"global chi-square({nb - 1} df) = {chi:.2f} p = {chi_p:.4f} -> " f"{'REJECTED' if chi_p < 0.01 else 'not rejected'} at 1%") if zmax_b >= 0: out(f"max |z| = {zmax:+.3f} at bucket {zmax_b} " f"(r {zmax_b / BINS:.3f}-{(zmax_b + 1) / BINS:.3f}), " f"freq {ones[zmax_b] / tot[zmax_b]:.4f}") out(f"buckets beyond Bonferroni |z| > {Z_A1_REL:.3f} (fw 0.01): " f"{zflags}/{nb}") out(f"MSB pinning: {msb_pinned}/{N} keys have bit n-1 = 1 -> forced by " f"2^(n-1) <= k (range arithmetic, not a generator footprint)") abs_rows = [] for p in range(160): np_, o = 0, 0 for r in solved: if r["n"] >= p + 2: # bit p present AND not that key's pinned MSB np_ += 1 o += (r["key"] >> p) & 1 if np_ >= 40: abs_rows.append((p, np_, o, (o / np_ - 0.5) / math.sqrt(0.25 / np_))) abs_zmax = max(abs_rows, key=lambda t: abs(t[3])) abs_flags = sum(1 for t in abs_rows if abs(t[3]) > Z_A1_ABS) out(f"absolute free-bit positions (>=40 keys each) : {len(abs_rows)} " f"(p {abs_rows[0][0]}..{abs_rows[-1][0]})") out(f" max |z| = {abs_zmax[3]:+.3f} at p={abs_zmax[0]} " f"({abs_zmax[2]}/{abs_zmax[1]} ones); flags beyond |z| > " f"{Z_A1_ABS:.3f}: {abs_flags}/{len(abs_rows)}") if chi_p < 0.01 or zflags or abs_flags: flags.append("A1 bit-position bias") res["a1"] = {"bins": BINS, "ones": ones, "tot": tot, "zs": zs, "chi": chi, "chip": chi_p, "zmax": zmax, "zmaxB": zmax_b, "zflags": zflags, "msb": msb_pinned, "absRows": abs_rows, "absZ": abs_zmax[3], "absP": abs_zmax[0], "absFlags": abs_flags, "bitsExcl": sum(tot), "totalBits": total_bits} # ---- check 2: spectral & frequency domain ------------------------ out("") out("-" * 74) out("ADVANCED CHECK 2 — SPECTRAL & FREQUENCY DOMAIN (hidden periodicities)") out("-" * 74) u = [float(k) / float(1 << n) - 0.5 for k, n in zip(keys, ns)] LAGS = 60 acf = [pearson(u[:N - L], u[L:]) for L in range(1, LAGS + 1)] acf_max_i = max(range(LAGS), key=lambda i: abs(acf[i])) acf_crit = math.tanh(Z_A2_ACF / math.sqrt((N - acf_max_i - 1) - 3)) acf_flags = 0 acf_rows = [] for i, r in enumerate(acf): nL = N - (i + 1) crit = math.tanh(Z_A2_ACF / math.sqrt(nL - 3)) hit = abs(r) > crit acf_flags += hit acf_rows.append((i + 1, nL, r, crit, hit)) out(f"value series u = k/2^n - 0.5, {N} points") out(f"multi-lag autocorrelation 1..{LAGS}: max |r| = {abs(acf[acf_max_i]):.3f} " f"at lag {acf_max_i + 1} (crit {acf_crit:.3f}, n_L {N - acf_max_i - 1})") out(f" lags beyond Bonferroni |r| > tanh(z/sqrt(n_L-3)) (fw 0.01): " f"{acf_flags}/{LAGS}") m_freq = (N - 1) // 2 powers = [] for k in range(1, m_freq + 1): re, im = 0.0, 0.0 w = 2.0 * math.pi * k / N for t, v in enumerate(u): re += v * math.cos(w * t) im -= v * math.sin(w * t) powers.append(re * re + im * im) p_mean = nsum(powers) / m_freq p_max = max(powers) k_peak = powers.index(p_max) + 1 fisher_y = p_max / p_mean fisher_p = 1.0 - (1.0 - math.exp(-fisher_y)) ** m_freq out(f"Fisher periodogram: {m_freq} ordinates, peak k={k_peak} " f"(period n/k = {N / k_peak:.2f}), max/mean power = {fisher_y:.3f}") out(f" Fisher spike probability p = {fisher_p:.4f} -> " f"{'SPECTRAL SPIKE' if fisher_p < 0.01 else 'no frequency spike'} at 1%") # concatenated bit stream (keys in order, MSB-first) total = sum(ns) stream, start = 0, 0 for k, n in zip(keys, ns): stream |= k << (total - start - n) start += n all_bits = bin(stream).count("1") SL = min(1024, total - 2) str_max, str_max_lag = 0.0, 0 str_flags = 0 blocks = [{"max": 0.0, "lag": 0, "crit": 0.0} for _ in range(8)] blk_flags = [False] * 8 for L in range(1, SL + 1): both = bin(stream & (stream >> L)).count("1") ones_f = all_bits - bin(stream & ((1 << L) - 1)).count("1") ones_s = bin(stream & ((1 << (total - L)) - 1)).count("1") nL = total - L s1 = (2 * ones_f - nL) / nL s2 = (2 * ones_s - nL) / nL mp = (4 * both - 2 * ones_f - 2 * ones_s + nL) / nL cov = mp - s1 * s2 r = cov / math.sqrt((1.0 - s1 * s1) * (1.0 - s2 * s2)) crit = math.tanh(Z_A2_STR / math.sqrt(nL - 3)) hit = abs(r) > crit str_flags += hit if abs(r) > abs(str_max): str_max, str_max_lag = r, L blk = (L - 1) // 128 if hit: blk_flags[blk] = True if abs(r) > abs(blocks[blk]["max"]): blocks[blk] = {"max": r, "lag": L, "crit": crit} for _b in range(8): blocks[_b]["flag"] = blk_flags[_b] out(f"bit stream: {total} bits ({N} keys, MSB-first), lags 1..{SL}") out(f" max |r| = {abs(str_max):.3f} at lag {str_max_lag} " f"(crit {math.tanh(Z_A2_STR / math.sqrt(total - str_max_lag - 3)):.3f}, " f"n_L {total - str_max_lag})") out(f" lags beyond Bonferroni |r| > tanh(z/sqrt(n_L-3)) (fw 0.01): " f"{str_flags}/{SL}") if acf_flags or fisher_p < 0.01 or str_flags: flags.append("A2 spectral / periodicity") res["a2"] = {"acf": acf, "acfRows": acf_rows, "acfMax": acf[acf_max_i], "acfMaxLag": acf_max_i + 1, "acfFlags": acf_flags, "mFreq": m_freq, "powers": powers, "kPeak": k_peak, "fisherY": fisher_y, "fisherP": fisher_p, "streamBits": total, "streamLags": SL, "strMax": str_max, "strMaxLag": str_max_lag, "strFlags": str_flags, "blocks": blocks} # ---- check 3: cross-key algebraic multipliers & relations -------- out("") out("-" * 74) out("ADVANCED CHECK 3 — CROSS-KEY ALGEBRAIC MULTIPLIERS & RELATIONS") out("-" * 74) t1_pairs, t1_exp, t1_hits = 0, 0.0, [] for i in range(N): for j in range(i + 1, N): if min(ns[i], ns[j]) < 16: continue # tiny keys make small multiples arithmetic t1_pairs += 1 lo_k, hi_k = (keys[i], keys[j]) if keys[i] < keys[j] \ else (keys[j], keys[i]) nhi = max(ns[i], ns[j]) lo_r, hi_r = 1 << (nhi - 1), (1 << nhi) - 1 qmin = max(2, -((-lo_r) // lo_k)) # ceil(lo_r / lo_k) qmax = min(4096, hi_r // lo_k) cand = max(0, qmax - qmin + 1) t1_exp += min(1.0, cand / float(hi_r - lo_r + 1)) if hi_k % lo_k == 0 and 2 <= hi_k // lo_k <= 4096: t1_hits.append((min(ns[i], ns[j]), nhi, hi_k // lo_k)) t1_p = poisson_tail(len(t1_hits), t1_exp) out(f"exact multiplier K_big = m*K_small, m in [2,4096], both >= 2^15:") out(f" pairs {t1_pairs}, hits {len(t1_hits)}, expected under H0 {t1_exp:.3g}, " f"P(X >= hits) = {t1_p:.4f}") for h in t1_hits[:5]: out(f" HIT: n={h[0]} key * {h[2]} = n={h[1]} key") t2 = {} for M, need, label in ((1 << 32, 33, "2^32"), (1 << 64, 65, "2^64")): tested, hits, excl_hits = 0, [], 0 for i in range(N): xi = keys[i] % M if xi & 1: inv = pow(xi, -1, M) for j in range(N): if j == i: continue m = (keys[j] % M) * inv % M if ns[j] < need: # K_j confined below 2^M: residue not uniform, the # uniform-m null would manufacture hits here if 1 <= m <= 65536: excl_hits += 1 continue tested += 1 if 1 <= m <= 65536: hits.append((ns[i], ns[j], m)) exp = tested * 65536 / M p_poi = poisson_tail(len(hits), exp) t2[label] = {"tested": tested, "hits": hits, "exp": exp, "p": p_poi, "excl": excl_hits, "need": need} out(f"modular multiplier K_j = m*K_i mod {label}, m in [1,65536]:") out(f" null-valid scan (n_j >= {need}): ordered pairs {tested}, " f"hits {len(hits)}, expected {exp:.3g}, P(X >= hits) = {p_poi:.4f}") out(f" excluded (n_j < {need}, confined residue): {excl_hits} flags — " f"trivial small-key relations, not evidence") for h in hits[:5]: out(f" HIT: #{h[0]} * m={h[2]} = #{h[1]} (mod {label})") keyset = set(keys) t3_pairs, t3_hits, t3_exp = 0, [], 0.0 for i in range(N): for j in range(i + 1, N): t3_pairs += 1 na, nb_ = ns[i], ns[j] s_lo = (1 << (na - 1)) + (1 << (nb_ - 1)) s_hi = (1 << na) + (1 << nb_) - 2 cand = 0 for k in keys: if k < s_lo: continue if k > s_hi: break cand += 1 t3_exp += cand / float(s_hi - s_lo + 1) if keys[i] + keys[j] in keyset: t3_hits.append((na, nb_)) t3_p = poisson_tail(len(t3_hits), t3_exp) out(f"sum triples K_a + K_b = K_c (a,b,c distinct):") out(f" pairs {t3_pairs}, hits {len(t3_hits)}, expected {t3_exp:.3g}, " f"P(X >= hits) = {t3_p:.4f}") for h in t3_hits[:5]: out(f" HIT: n={h[0]} + n={h[1]} equals another solved key") pivots = {} for k in keys: v = k while v: p = v.bit_length() - 1 if p in pivots: v ^= pivots[p] else: pivots[p] = v break rank = len(pivots) out(f"GF(2) linear rank of the {N} key vectors (width {maxn}): " f"{rank}/{N} -> " + ("linearly independent" if rank == N else "LINEAR DEPENDENCY EXISTS")) if (t1_p < 0.01 or any(v["p"] < 0.01 for v in t2.values()) or t3_p < 0.01 or rank < N): flags.append("A3 algebraic relation") res["a3"] = {"t1Pairs": t1_pairs, "t1Hits": t1_hits, "t1Exp": t1_exp, "t1P": t1_p, "t2": t2, "t3Pairs": t3_pairs, "t3Hits": t3_hits, "t3Exp": t3_exp, "t3P": t3_p, "rank": rank, "width": maxn} # ---- check 4: suffix, nibble & byte entropy profiling ------------ out("") out("-" * 74) out("ADVANCED CHECK 4 — SUFFIX, NIBBLE & BYTE ENTROPY PROFILING") out("-" * 74) maxhl = max(hexlens) nib = [] for j in range(maxhl): cnt = [0] * 16 np_ = 0 for r, hl in zip(solved, hexlens): if hl > j: h = format(r["key"], "0{}x".format(hl)) cnt[int(h[len(h) - 1 - j], 16)] += 1 np_ += 1 nib.append((np_, cnt)) pool_js = [j for j, (np_, _) in enumerate(nib) if np_ >= 64] pool_cnt = [0] * 16 pool_n = 0 for j in pool_js: np_, cnt = nib[j] pool_n += np_ for v in range(16): pool_cnt[v] += cnt[v] pool_e = pool_n / 16 pool_chi = nsum(((pool_cnt[v] - pool_e) * (pool_cnt[v] - pool_e)) / pool_e for v in range(16)) pool_p = chi2sf(pool_chi, 15) out(f"LSB-aligned nibble positions j=0..{maxhl - 1} " f"(j=0 = least significant hex digit)") out(f" pooled j={pool_js[0]}..{pool_js[-1]} ({pool_n} samples, " f"E {pool_e:.1f}/category): chi-square(15) = {pool_chi:.2f}, " f"p = {pool_p:.4f} -> {'REJECTED' if pool_p < 0.01 else 'not rejected'}") np0, cnt0 = nib[0] e0 = np0 / 16 chi0 = nsum(((cnt0[v] - e0) * (cnt0[v] - e0)) / e0 for v in range(16)) p0 = chi2sf(chi0, 15) h0, m0 = hex_entropy(cnt0, np0) out(f" position j=0 (LSB nibble): n={np0}, E={e0:.1f}/category, " f"chi-square(15) = {chi0:.2f}, p = {p0:.4f} -> " f"{'REJECTED' if p0 < 0.01 else 'not rejected'}; entropy {h0:.3f} bits") power_ok = sum(1 for np_, _ in nib if np_ / 16 >= 5.0) out(f" per-position chi-square with E >= 5: {power_ok} of {len(nib)} " f"positions (rest reported with n and entropy, low power)") top_cnt = [0] * 16 top_ok = 0 for r, n in zip(solved, ns): h = format(r["key"], "x") d = int(h[0], 16) top_cnt[d] += 1 allowed = {0: "89abcdef", 1: "1", 2: "23", 3: "4567"}[n % 4] top_ok += h[0] in allowed top_h, _ = hex_entropy(top_cnt, N) out(f"MSB significant nibble: {top_ok}/{N} inside the tight range-forced " f"set (n mod 4), entropy {top_h:.3f} bits — range mask, not a seed") bytes_ = [] for b in range(32): cnt = {} zeros = 0 for r in solved: h = format(r["key"], "064x") v = int(h[62 - 2 * b:64 - 2 * b], 16) cnt[v] = cnt.get(v, 0) + 1 zeros += v == 0 h_b, _ = hex_entropy([cnt.get(v, 0) for v in range(256)], N) bytes_.append((zeros, h_b)) pad_from = next((b for b in range(32) if all(bytes_[bb][0] == N for bb in range(b, 32))), None) out(f"byte view (32 big-endian bytes): LSB byte entropy " f"{bytes_[0][1]:.3f} bits (plug-in, n={N}); every byte b >= {pad_from} " f"is zero in all {N} keys -> key-length padding (max n = {maxn} bits)") if pool_p < 0.01 or p0 < 0.01 or top_ok < N: flags.append("A4 nibble/byte entropy") res["a4"] = {"nib": nib, "poolJs": pool_js, "poolN": pool_n, "poolChi": pool_chi, "poolP": pool_p, "j0": {"n": np0, "cnt": cnt0, "chi": chi0, "p": p0, "h": h0}, "topCnt": top_cnt, "topOk": top_ok, "topH": top_h, "bytes": bytes_, "padFrom": pad_from} # ---- check 5: Benford's law -------------------------------------- out("") out("-" * 74) out("ADVANCED CHECK 5 — BENFORD'S LAW COMPLIANCE") out("-" * 74) dec = [0] * 9 for k in keys: dec[int(str(k)[0]) - 1] += 1 ben_d = [math.log10(1.0 + 1.0 / d) for d in range(1, 10)] chi_d = nsum(((dec[d] - N * ben_d[d]) * (dec[d] - N * ben_d[d])) / (N * ben_d[d]) for d in range(9)) p_d = chi2sf(chi_d, 8) mad_d = nsum(abs(dec[d] / N - ben_d[d]) for d in range(9)) / 9 out(f"decimal leading digit (1-9): counts {dec}") out(f" chi-square(8) = {chi_d:.2f}, p = {p_d:.4f} -> " f"{'REJECTED' if p_d < 0.01 else 'not rejected'} at 1%; " f"MAD = {mad_d:.4f} (uniform-random MAD would be ~0.0556)") hx = [0] * 15 for k in keys: hx[int(format(k, "x")[0], 16) - 1] += 1 ben_h = [math.log(1.0 + 1.0 / d) / math.log(16.0) for d in range(1, 16)] chi_h = nsum(((hx[d - 1] - N * ben_h[d - 1]) * (hx[d - 1] - N * ben_h[d - 1])) / (N * ben_h[d - 1]) for d in range(1, 16)) p_h = chi2sf(chi_h, 14) mad_h = nsum(abs(hx[d - 1] / N - ben_h[d - 1]) for d in range(1, 16)) / 15 min_e = min(N * b for b in ben_h) out(f"hex leading digit (1-f): counts {hx}") out(f" chi-square(14) = {chi_h:.2f}, p = {p_h:.4f} -> " f"{'REJECTED' if p_h < 0.01 else 'not rejected'} at 1%; " f"MAD = {mad_h:.4f}; min expected {min_e:.2f} (low power caveat)") if p_d < 0.01 or p_h < 0.01: flags.append("A5 Benford") res["a5"] = {"dec": dec, "chiD": chi_d, "pD": p_d, "madD": mad_d, "hx": hx, "chiH": chi_h, "pH": p_h, "madH": mad_h, "minE": min_e} # ---- advanced verdict -------------------------------------------- out("") out("-" * 74) out("ADVANCED VERDICT (family-wise alpha = 0.01 per check)") out("-" * 74) if flags: out(f"FLAGS: {len(flags)} of 5 checks flagged -> " + ", ".join(flags)) out("Inspect the flagged check above before trusting any structure claim;") out("a single hit at alpha 0.01 can still be a false positive.") else: out("FLAGS: 0 of 5 checks flagged.") out("The null hypothesis (true cryptographic pseudo-randomness) holds") out("across all five advanced metrics: no bit-position pinning beyond") out("the range-forced MSB, no periodic or spectral spike beyond the") out("Bonferroni thresholds, no cross-key algebraic relation beyond the") out("chance expectation, nibble/byte entropy at chance level, and") out("Benford compliance within sampling noise. No generator footprint") out("and no non-random structural anomaly is detected.") res["flags"] = flags return res # -------------------------------------------------------------------------- # five low-level bitwise / algebraic / residue checks (solved keys only) # -------------------------------------------------------------------------- def norm_two_sided_p(z): """Two-sided normal p-value: erfc(|z|/sqrt2) = Q(1/2, z^2/2).""" return gammainc_upper(0.5, z * z / 2.0) def z_for_p(p): """Bonferroni z threshold for a per-comparison two-sided alpha p.""" lo, hi = 0.0, 40.0 for _ in range(80): mid = 0.5 * (lo + hi) if norm_two_sided_p(mid) > p: lo = mid else: hi = mid return 0.5 * (lo + hi) def binom_pmf(n, k, p=0.5): if k < 0 or k > n: return 0.0 lg = lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1) if p == 0.5: return math.exp(lg - n * math.log(2.0)) return math.exp(lg + k * math.log(p) + (n - k) * math.log(1.0 - p)) def binom_tail(k, n, p=0.5): """P(X >= k) for X ~ Binomial(n, p).""" if k <= 0: return 1.0 if k > n: return 0.0 s = 0.0 for i in range(k, n + 1): s += binom_pmf(n, i, p) return min(1.0, s) def bm_lc(s): """Berlekamp-Massey linear complexity of a GF(2) sequence (LFSR order).""" C, B, L, m = [1], [1], 0, 1 for n in range(len(s)): d = s[n] for i in range(1, L + 1): d ^= C[i] & s[n - i] if d: T = C[:] if len(C) < len(B) + m: C += [0] * (len(B) + m - len(C)) for j in range(len(B)): C[j + m] ^= B[j] if 2 * L <= n: L, B, m = n + 1 - L, T, 1 else: m += 1 else: m += 1 return L def max_bit_run(x, nbits): """Longest run of identical bits in the nbits-bit pattern of x.""" best, run, prev = 0, 0, -1 for p in range(nbits): b = (x >> p) & 1 run = run + 1 if b == prev else 1 prev = b if run > best: best = run return best def pmaxrun_ge(L, r): """P(longest run of identical bits >= r) for L fair (Bernoulli 1/2) bits.""" if r <= 1 or L < r: return 1.0 if (r <= 1 and L >= 1) else 0.0 cur = [0.0] * r cur[1] = 1.0 for _ in range(1, L): nxt = [0.0] * r for j in range(1, r): nxt[1] += cur[j] * 0.5 if j + 1 < r: nxt[j + 1] += cur[j] * 0.5 cur = nxt return 1.0 - sum(cur[1:]) def rotate_bits(x, nbits, t): """Rotate the nbits-bit pattern of x left by t positions.""" mask = (1 << nbits) - 1 t %= nbits return ((x << t) | (x >> (nbits - t))) & mask def lowlevel(solved, out=print): """Checks 6-10 of the low-level forensic protocol; returns measured results. Same contract as advanced(): every number is computed from the solved keys, each check carries a family-wise alpha of 0.01, and the verdict states whether any of them flags. """ N = len(solved) flags = [] res = {} rows = sorted(solved, key=lambda r: r["n"]) keys = [r["key"] for r in rows] ns = [r["n"] for r in rows] # ---- check 6: bitwise XOR gap analysis --------------------------- out("") out("-" * 74) out("LOW-LEVEL CHECK 6 — BITWISE XOR GAP ANALYSIS K[i+1] XOR K[i]") out("-" * 74) GAPS = N - 1 gaps = [] for i in range(GAPS): ka, na = keys[i], ns[i] kb, nb = keys[i + 1], ns[i + 1] x = ka ^ kb L = max(na, nb) O = min(na, nb) same = 1 if na == nb else 0 ext = L - O free = O - same e = 0.5 * free + (0.5 * (ext - 1) + 1.0 if ext > 0 else 0.0) v = 0.25 * free + (0.25 * (ext - 1) if ext > 0 else 0.0) pop = popcount(x) z = (pop - e) / math.sqrt(v) if v > 0 else 0.0 gaps.append({"i": i, "na": na, "nb": nb, "L": L, "x": x, "pop": pop, "exp": e, "var": v, "z": z, "run": max_bit_run(x, L)}) zsum = 0.0 vsum = 0.0 for g in gaps: zsum += g["pop"] - g["exp"] vsum += g["var"] ztot = zsum / math.sqrt(vsum) ptot = norm_two_sided_p(ztot) zc6 = z_for_p(0.01 / GAPS) zflags = [g for g in gaps if abs(g["z"]) > zc6] nflags = len(zflags) zmaxg = max(gaps, key=lambda g: abs(g["z"])) out(f"consecutive gaps (n-sorted solved keys) : {GAPS}") out(f"mean XOR popcount {mean([g['pop'] for g in gaps]):.3f} vs " f"null expectation {mean([g['exp'] for g in gaps]):.3f}") out(f"summed z = {ztot:+.3f} (two-sided p = {ptot:.4f}) -> " f"{'REJECTED' if ptot < 0.01 else 'not rejected'} at 1%") out(f"max |z| = {zmaxg['z']:+.3f} at gap {zmaxg['i']} " f"(n {zmaxg['na']}->{zmaxg['nb']}, popcount {zmaxg['pop']} vs " f"{zmaxg['exp']:.2f}); gaps beyond Bonferroni |z| > {zc6:.3f}: " f"{nflags}/{GAPS}") runs = [g["run"] for g in gaps] runmax = max(runs) runmax_g = max(gaps, key=lambda g: g["run"]) exp_runs = sum(pmaxrun_ge(g["L"], runmax) for g in gaps) notany = 1.0 for g in gaps: notany *= 1.0 - pmaxrun_ge(g["L"], runmax) p_runmax = 1.0 - notany out(f"longest identical-bit run inside an XOR: {runmax} bits " f"(gap {runmax_g['i']}, n {runmax_g['na']}->{runmax_g['nb']}); " f"null: {exp_runs:.2f} gaps expected to reach it, P(any) = " f"{p_runmax:.4f}") big = [g for g in gaps if g["L"] >= 16] per_hits, rot_hits = [], [] for g in big: x, L = g["x"], g["L"] found = None for t in range(1, L): if rotate_bits(x, L, t) == x: found = t break if found is not None: rot_hits.append({"i": g["i"], "t": found}) if L % found == 0: per_hits.append({"i": g["i"], "period": found}) out(f"XOR values with a nontrivial cyclic period (L>=16): {len(per_hits)}/" f"{len(big)}; rotation-symmetric: {len(rot_hits)}/{len(big)}") PLANES = 14 planes = [] for p in range(PLANES): seq = [(g["x"] >> p) & 1 for g in gaps] lc = bm_lc(seq) planes.append({"p": p, "lc": lc, "n": len(seq), "ratio": lc / float(len(seq))}) lc_min = min(planes, key=lambda r: r["ratio"]) lc_lim = 0.30 lc_flags = [r for r in planes if r["ratio"] <= lc_lim] out(f"LFSR / Berlekamp-Massey linear complexity of the XOR bit planes " f"p=0..{PLANES - 1} ({GAPS} samples): min {lc_min['lc']}/{GAPS} " f"(ratio {lc_min['ratio']:.3f} at p={lc_min['p']}); planes with " f"LC <= {lc_lim:.2f}L: {len(lc_flags)}/{PLANES}") a6on = ptot < 0.01 or nflags > 0 or len(lc_flags) > 0 if a6on: flags.append("L6 XOR gap analysis") res["l6"] = {"gaps": gaps, "n": GAPS, "ztot": ztot, "ptot": ptot, "zcrit": zc6, "zflags": nflags, "zmax": zmaxg["z"], "zmaxI": zmaxg["i"], "meanPop": mean([g["pop"] for g in gaps]), "meanExp": mean([g["exp"] for g in gaps]), "runMax": runmax, "runMaxI": runmax_g["i"], "expRuns": exp_runs, "pRun": p_runmax, "periodHits": per_hits, "rotHits": rot_hits, "big": len(big), "planes": planes, "lcMin": lc_min["lc"], "lcMinP": lc_min["p"], "lcMinRatio": lc_min["ratio"], "lcFlags": len(lc_flags), "lcLim": lc_lim, "flag": a6on} # ---- check 7: GCD gap profiling --------------------------------- out("") out("-" * 74) out("LOW-LEVEL CHECK 7 — GCD GAP PROFILING (shared factor in the diffs)") out("-" * 74) diffs = [keys[i + 1] - keys[i] for i in range(GAPS)] g_all = 0 for d in diffs: g_all = math.gcd(g_all, d) g_keys = 0 for k in keys: g_keys = math.gcd(g_keys, k) succ = [math.gcd(diffs[i], diffs[i + 1]) for i in range(GAPS - 1)] c_succ = sum(1 for g in succ if g > 1) p_coprime = 1.0 - 6.0 / (math.pi * math.pi) e_succ = len(succ) * p_coprime v_succ = len(succ) * p_coprime * (1.0 - p_coprime) z_succ = (c_succ - e_succ) / math.sqrt(v_succ) p_succ = norm_two_sided_p(z_succ) even = sum(1 for d in diffs if d % 2 == 0) z_even = (even - GAPS * 0.5) / math.sqrt(GAPS * 0.25) p_even = norm_two_sided_p(z_even) out(f"gcd of all {GAPS} consecutive key differences : {g_all}") out(f"gcd of all {N} solved keys : {g_keys}") out(f"consecutive-difference gcd > 1 : {c_succ}/{len(succ)} " f"(null expectation {e_succ:.2f}, z {z_succ:+.3f}, two-sided p " f"{p_succ:.4f}) -> {'REJECTED' if p_succ < 0.01 else 'not rejected'}" f" at 1% (approximate: neighbouring pairs share one difference)") out(f"even differences : {even}/{GAPS} (z {z_even:+.3f}, p {p_even:.4f})") divrows = [] for m in range(2, 201): c = sum(1 for d in diffs if d % m == 0) e = GAPS / float(m) v = GAPS * (1.0 / m) * (1.0 - 1.0 / m) divrows.append({"m": m, "cnt": c, "exp": e, "z": (c - e) / math.sqrt(v) if v > 0 else 0.0}) zc7 = z_for_p(0.01 / len(divrows)) dflags = [r for r in divrows if abs(r["z"]) > zc7] dmax = max(divrows, key=lambda r: abs(r["z"])) out(f"divisibility scan m = 2..200 over the differences: max |z| " f"{dmax['z']:+.3f} at m={dmax['m']} ({dmax['cnt']} vs " f"{dmax['exp']:.2f} expected); beyond Bonferroni |z| > {zc7:.3f}: " f"{len(dflags)}/{len(divrows)}") g_succ_max = max(succ) if succ else 0 succ_top = sorted(enumerate(succ), key=lambda t: -t[1])[:10] a7on = p_succ < 0.01 or len(dflags) > 0 or p_even < 0.01 if a7on: flags.append("L7 GCD gap profile") res["l7"] = {"gAll": g_all, "gKeys": g_keys, "succN": len(succ), "n": GAPS, "succHits": c_succ, "succExp": e_succ, "zSucc": z_succ, "succP": p_succ, "succGcd": g_succ_max, "succTop": [{"i": i, "g": g} for i, g in succ_top], "even": even, "zEven": z_even, "pEven": p_even, "divRows": divrows, "zcrit": zc7, "dFlags": len(dflags), "dMaxZ": dmax["z"], "dMaxM": dmax["m"], "flag": a7on} # ---- check 8: popcount distribution shape ----------------------- out("") out("-" * 74) out("LOW-LEVEL CHECK 8 — FULL HAMMING WEIGHT (POPCOUNT) DISTRIBUTION") out("-" * 74) pops = [popcount(k) for k in keys] pz = [(pp - n / 2.0) / math.sqrt(n / 4.0) for pp, n in zip(pops, ns)] mz = mean(pz) sz = sd(pz) sk = sum((x - mz) ** 3 for x in pz) / (N * sz ** 3) ku = sum((x - mz) ** 4 for x in pz) / (N * sz ** 4) - 3.0 BINS8 = 10 obs8 = [0] * BINS8 exp8 = [0.0] * BINS8 for pp, n in zip(pops, ns): obs8[min(BINS8 - 1, (pp * BINS8) // n)] += 1 for k in range(n + 1): exp8[min(BINS8 - 1, (k * BINS8) // n)] += binom_pmf(n, k) CUT = (-1.0, 0.0, 1.0) NCAT = len(CUT) + 1 def zcat(v, n): zz = (v - n / 2.0) / math.sqrt(n / 4.0) c = 0 while c < len(CUT) and zz > CUT[c]: c += 1 return c obs8c = [0] * NCAT exp8c = [0.0] * NCAT exp_tail8 = 0.0 for pp, n in zip(pops, ns): obs8c[zcat(pp, n)] += 1 for k in range(n + 1): pm = binom_pmf(n, k) exp8c[zcat(k, n)] += pm if abs((k - n / 2.0) / math.sqrt(n / 4.0)) > 2.0: exp_tail8 += pm chi8 = 0.0 for c in range(NCAT): chi8 += (obs8c[c] - exp8c[c]) ** 2 / exp8c[c] p8 = chi2sf(chi8, NCAT - 1) out(f"standardised popcount z=(pop-n/2)/sqrt(n/4): mean {mz:+.4f} " f"(expect 0), sd {sz:.4f} (expect 1), skew {sk:+.4f} (0), " f"excess kurtosis {ku:+.4f} (0)") out(f"standardised-popcount bands (|z|<=1 core) vs binomial B(n,1/2): " f"chi-square({NCAT - 1} df over 4 z-bands) = {chi8:.2f} " f"p = {p8:.4f} -> {'REJECTED' if p8 < 0.01 else 'not rejected'} " f"at 1% (obs {obs8c} vs exp " f"{[round(x, 2) for x in exp8c]})") z2 = sum(1 for x in pz if abs(x) > 2.0) z3 = sum(1 for x in pz if abs(x) > 3.0) p_z2 = poisson_tail(z2, exp_tail8) out(f"keys with |z| > 2 : {z2}/{N} (exact binomial expectation " f"{exp_tail8:.2f}, P {p_z2:.4f}); |z| > 3 : {z3}/{N}") evn = sum(1 for pp in pops if pp % 2 == 0) z_par = (evn - N * 0.5) / math.sqrt(N * 0.25) p_par = norm_two_sided_p(z_par) out(f"even popcounts : {evn}/{N} (z {z_par:+.3f}, p {p_par:.4f}) -> " f"{'parity bias' if p_par < 0.01 else 'no parity bias'}") a8on = p8 < 0.01 or p_par < 0.01 or p_z2 < 0.01 if a8on: flags.append("L8 popcount shape") res["l8"] = {"meanZ": mz, "sdZ": sz, "skew": sk, "kurt": ku, "obs": obs8, "exp": exp8, "chi": chi8, "p": p8, "df": NCAT - 1, "obsCat": obs8c, "expCat": exp8c, "z2": z2, "z3": z3, "expZ2": exp_tail8, "pZ2": p_z2, "even": evn, "zPar": z_par, "pPar": p_par, "pops": pops, "pz": pz, "flag": a8on} # ---- check 9: decimal (base-10) artifacts ----------------------- out("") out("-" * 74) out("LOW-LEVEL CHECK 9 — DECIMAL (BASE-10) CONVERSION ARTIFACTS") out("-" * 74) lastd = [0] * 10 for k in keys: lastd[k % 10] += 1 e9 = N / 10.0 chi9 = sum((c - e9) ** 2 / e9 for c in lastd) p9 = chi2sf(chi9, 9) dlast = [0] * 10 for d in diffs: dlast[d % 10] += 1 e9b = GAPS / 10.0 chi9b = sum((c - e9b) ** 2 / e9b for c in dlast) p9b = chi2sf(chi9b, 9) lz = [k % 100 for k in keys] distinct2 = len(set(lz)) exp_distinct = 100.0 * (1.0 - (99.0 / 100.0) ** N) zeros = lastd[0] bigrams = 0 exp_bigram = 0.0 triples = 0 exp_triple = 0.0 for d in diffs: s = str(d) Ld = len(s) # tandem block repeat: the same 2-digit block twice in a row if any(s[j:j + 2] == s[j + 2:j + 4] for j in range(max(0, Ld - 3))): bigrams += 1 if Ld >= 4: exp_bigram += 1.0 - 0.99 ** (Ld - 3) # three identical digits in a row if any(s[j] == s[j + 1] == s[j + 2] for j in range(max(0, Ld - 2))): triples += 1 if Ld >= 3: exp_triple += 1.0 - 0.99 ** (Ld - 2) p_big = poisson_tail(bigrams, exp_bigram) p_tri = poisson_tail(triples, exp_triple) ascii_ok = [] exp_ascii = 0.0 for k, n in zip(keys, ns): nb = (n + 7) // 8 exp_ascii += (95.0 / 256.0) ** nb raw = k.to_bytes(nb, "big") if all(32 <= c < 127 for c in raw): ascii_ok.append({"n": n, "key": k}) big_ascii = [a for a in ascii_ok if a["n"] >= 32] p_ascii = poisson_tail(len(ascii_ok), exp_ascii) out(f"last decimal digit of the keys : chi-square(9) = {chi9:.2f} " f"p = {p9:.4f} -> {'REJECTED' if p9 < 0.01 else 'not rejected'} " f"at 1% (counts {lastd})") out(f"last decimal digit of the gaps : chi-square(9) = {chi9b:.2f} " f"p = {p9b:.4f} -> {'REJECTED' if p9b < 0.01 else 'not rejected'} " f"at 1% (counts {dlast})") out(f"last two digits: {distinct2} distinct of 100 (null " f"{exp_distinct:.2f}); trailing zero: {zeros}/{N} (expect {e9:.2f})") out(f"tandem digit repeats inside the gaps: repeated 2-digit block " f"(s[j:j+2] = s[j+2:j+4]) in {bigrams} vs {exp_bigram:.2f} expected " f"(P {p_big:.4f}); three identical digits in a row {triples} vs " f"{exp_triple:.2f} (P {p_tri:.4f})") out(f"keys whose byte image is fully printable ASCII: {len(ascii_ok)}/" f"{N} vs {exp_ascii:.2f} expected from (95/256)^bytes " f"(P {p_ascii:.4f}); with n >= 32: {len(big_ascii)} -> " f"{'text-to-number seed signature' if big_ascii else 'none'}") a9on = (p9 < 0.01 or p9b < 0.01 or p_big < 0.01 or p_tri < 0.01 or len(big_ascii) > 0 or p_ascii < 0.01) if a9on: flags.append("L9 decimal artifacts") res["l9"] = {"lastd": lastd, "chi": chi9, "p": p9, "dlast": dlast, "chiB": chi9b, "pB": p9b, "distinct2": distinct2, "expDistinct": exp_distinct, "zeros": zeros, "bigrams": bigrams, "expBigram": exp_bigram, "pBigram": p_big, "triples": triples, "expTriple": exp_triple, "pTriple": p_tri, "ascii": len(ascii_ok), "expAscii": exp_ascii, "pAscii": p_ascii, "asciiBig": len(big_ascii), "asciiKeys": ascii_ok, "flag": a9on} # ---- check 10: forbidden residues across small primes ------------ out("") out("-" * 74) out("LOW-LEVEL CHECK 10 — FORBIDDEN RESIDUES ACROSS SMALL PRIMES") out("-" * 74) PRIMES = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31] mrows = [] tot_buckets = sum(PRIMES) tot_empty = 0 tot_exp_empty = 0.0 forbidden = [] all_rz = [] for M in PRIMES: cnt = [0] * M for k in keys: cnt[k % M] += 1 e = N / float(M) chi = sum((c - e) ** 2 / e for c in cnt) p = chi2sf(chi, M - 1) emp = sum(1 for c in cnt if c == 0) ee = M * (1.0 - 1.0 / M) ** N pe = poisson_tail(emp, ee) zs = [(c - e) / math.sqrt(e * (1.0 - 1.0 / M)) for c in cnt] zmin = min(zs) empties = [r for r in range(M) if cnt[r] == 0] for r in empties: forbidden.append({"M": M, "r": r}) for r in range(M): all_rz.append({"M": M, "r": r, "z": zs[r], "cnt": cnt[r]}) tot_empty += emp tot_exp_empty += ee mrows.append({"M": M, "cnt": cnt, "chi": chi, "p": p, "df": M - 1, "empty": emp, "expEmpty": ee, "pEmpty": pe, "zMin": zmin, "zMinR": zs.index(zmin), "empties": empties, "minCnt": min(cnt)}) out(f"mod {M:>2} : chi-square({M - 1:>2} df) = {chi:5.2f} " f"p = {p:.4f} empty buckets {emp}/{M} (null {ee:.2f}, " f"P {pe:.4f}) lowest z {zmin:+.2f} at r={zs.index(zmin)}") zc10 = z_for_p(0.01 / tot_buckets) rz_flags = [r for r in all_rz if abs(r["z"]) > zc10] rz_min = min(all_rz, key=lambda r: r["z"]) p_empty_tot = poisson_tail(tot_empty, tot_exp_empty) pmin = min(r["p"] for r in mrows) fstr = " ".join("%d:%d" % (f["M"], f["r"]) for f in forbidden) or "none" out(f"forbidden residues (count 0 over all {N} keys) : {len(forbidden)}" f" [{fstr}]") out(f"empty buckets summed over the 10 moduli : {tot_empty}/" f"{tot_buckets} (null {tot_exp_empty:.2f}, P {p_empty_tot:.4f})") out(f"residue-bucket z scan ({tot_buckets} buckets): max |z| " f"{abs(rz_min['z']):.3f} at mod {rz_min['M']} r={rz_min['r']}; beyond " f"Bonferroni |z| > {zc10:.3f}: {len(rz_flags)}/{tot_buckets}") out(f"smallest per-modulus p = {pmin:.4f} " f"(Bonferroni alpha {0.01 / len(PRIMES):.4f})") a10on = (pmin < 0.01 / len(PRIMES) or p_empty_tot < 0.01 or len(rz_flags) > 0) if a10on: flags.append("L10 forbidden residues") res["l10"] = {"primes": PRIMES, "rows": mrows, "total": tot_buckets, "empty": tot_empty, "expEmpty": tot_exp_empty, "pEmpty": p_empty_tot, "zcrit": zc10, "rzFlags": len(rz_flags), "rzMinZ": rz_min["z"], "rzMinM": rz_min["M"], "rzMinR": rz_min["r"], "forbidden": forbidden, "pMin": pmin, "flag": a10on} out("") out("=" * 74) out("LOW-LEVEL VERDICT (family-wise alpha = 0.01 per check)") out("=" * 74) if flags: out(f"FLAGS: {len(flags)} of 5 checks flagged: {', '.join(flags)}") out("At least one low-level statistic exceeded its null expectation;") out("inspect the flagged row before asserting a generator artifact.") else: out("FLAGS: 0 of 5 checks flagged.") out("No bitwise pattern, shared GCD factor, forbidden residue, digit") out("periodicity or popcount-shape anomaly was uncovered: the XOR gaps") out("and popcounts are binomial-consistent, the gap gcds carry no common") out("step factor, the decimal digit streams are uniform, and the only") out("empty residue classes are as few as chance predicts.") out("The null hypothesis (true cryptographic pseudo-randomness) holds") out("across all five low-level metrics — no deterministic generator") out("footprint is detected.") res["flags"] = flags return res # -------------------------------------------------------------------------- # five deep cryptographic / dynamic / information-theoretic checks # -------------------------------------------------------------------------- def f(v, d): """Fixed-point string with d decimals — the same digits math-deep.js prints through MC.f (page and script stay visually identical).""" return "%.*f" % (d, v) def sg(v, d): """Signed fixed-point string: a leading + on non-negative values.""" return "%+.*f" % (d, v) P256 = (1 << 256) - (1 << 32) - 977 # secp256k1 field prime SECP_ORDER = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141 RUNS_MIN_N = 16 # NIST runs applicability floor for a single key R_LEN = 200 # Monte-Carlo null draws per key bit length R_POOL = 500 # ... for the aggregated bitstream (key-structured null) R_BLOCK = 300 # ... for the block mutual-information null R_DATA = 500 # ... whole-dataset draws for check 11's matrix statistics R_IMG = 500 # ... whole-dataset draws for check 15's field images EXACT_LZ_N = 20 # largest bit length whose LZ78 null is enumerated exactly def erfc(x): """erfc(x) = Q(1/2, x^2) — the same gammainc_upper math-deep.js uses.""" if x <= 0.0: return 1.0 if math.isinf(x) or x != x: return 0.0 return gammainc_upper(0.5, x * x) def bits_msb(x, n): """The n-bit pattern of x, most significant bit first (the bitstream).""" return [(x >> p) & 1 for p in range(n - 1, -1, -1)] def nist_runs(bits): """NIST SP 800-22 runs test, exactly as the reference runs.c computes it. Returns (ok, pi, V, p). ok is False when step 2 (|pi-0.5| >= 2/sqrt n) makes the test not applicable, where NIST sets p = 0. """ n = len(bits) if n == 0: return False, 0.0, 0, 0.0 s = 0 for b in bits: s += b pi = s / n V = 1 for i in range(1, n): if bits[i] != bits[i - 1]: V += 1 if abs(pi - 0.5) >= 2.0 / math.sqrt(n): return False, pi, V, 0.0 if pi * (1.0 - pi) == 0.0: return True, pi, V, 0.0 arg = abs(V - 2.0 * n * pi * (1.0 - pi)) / (2.0 * pi * (1.0 - pi) * math.sqrt(2.0 * n)) return True, pi, V, erfc(arg) def lz78(bits): """Phrases in the LZ78 parse; the trailing partial phrase counts as one. The dictionary is keyed on the packed integer (state << 1) | bit, which is bijective with the (state, bit) pair -- the same parse, roughly twice as fast in the exact enumeration below, and mirrored in math-deep.js. """ d = {} nxt = 1 cur = 0 ph = 0 for b in bits: key = (cur << 1) | b hit = d.get(key) if hit is not None: cur = hit else: d[key] = nxt nxt += 1 ph += 1 cur = 0 if cur != 0: ph += 1 return ph def xs32(s): """xorshift32 (period 2^32-1) — bit-for-bit identical in math-deep.js.""" s &= 0xFFFFFFFF s ^= (s << 13) & 0xFFFFFFFF s ^= s >> 17 s ^= (s << 5) & 0xFFFFFFFF return s & 0xFFFFFFFF def rand_int(state, nbits): """A uniform n-bit integer with the MSB set (the key shape).""" v = 0 for _ in range(nbits - 1): state = xs32(state) v = (v << 1) | ((state >> 24) & 1) return v | (1 << (nbits - 1)), state def collatz_stats(x): """(total stopping time, odd-step proportion, parity-sequence runs).""" t = odd = runs = 0 prev = -1 while x != 1: p = x & 1 if p: odd += 1 if p != prev: runs += 1 prev = p x = 3 * x + 1 if p else x >> 1 t += 1 return t, (odd / t if t else 0.0), runs def rule30_stats(x, nbits, nrows): """Rule-30 spacetime: (mean 1-density, repeated rows) over `nrows` rows. Cell d is bit d (MSB at the left), left neighbour = cell d+1, so one row is (row >> 1) ^ (row | (row << 1)) with zero boundaries. """ mask = (1 << nbits) - 1 row = x & mask ones = cells = 0 seen = set() for _ in range(nrows): row = ((row >> 1) ^ (row | (row << 1))) & mask pc = row.bit_count() ones += pc cells += nbits seen.add(row) return (ones / cells if cells else 0.0), nrows - len(seen) def window_collisions(bits, width=8): """(colliding pairs of width-bit windows, exact mean under independence).""" L = len(bits) if L < width: return 0, 0.0 m = L - width + 1 cnt = {} w = 0 for j in range(width): w = (w << 1) | bits[j] cnt[w] = 1 for j in range(1, m): w = ((w << 1) | bits[j + width - 1]) & 255 cnt[w] = cnt.get(w, 0) + 1 c = 0 for v in cnt.values(): c += v * (v - 1) // 2 return c, (m * (m - 1) / 2.0) / 256.0 def mi_bits(tbl, ntot): """Plugin Shannon mutual information (bits) of a two-way count table.""" rows, cols = len(tbl), len(tbl[0]) rs = [0.0] * rows cs = [0.0] * cols for i in range(rows): for j in range(cols): rs[i] += tbl[i][j] cs[j] += tbl[i][j] mi = 0.0 for i in range(rows): for j in range(cols): o = tbl[i][j] if o: mi += (o / ntot) * math.log2((o * ntot) / (rs[i] * cs[j])) return mi def binom_two_sided(k, n): """Exact two-sided p for X ~ Binomial(n, 1/2): P(X<=k)+P(X>=n-k).""" if 2 * k == n: return 1.0 if 2 * k < n: lo = nsum([binom_pmf(n, i, 0.5) for i in range(0, k + 1)]) hi = nsum([binom_pmf(n, i, 0.5) for i in range(n - k, n + 1)]) else: lo = nsum([binom_pmf(n, i, 0.5) for i in range(0, n - k + 1)]) hi = nsum([binom_pmf(n, i, 0.5) for i in range(k, n + 1)]) return min(1.0, lo + hi) def deep_null_len(n, r): """Monte-Carlo null for one bit length: (mean, sd) of the six statistics that have no closed form (Collatz length, odd share, parity runs, Rule-30 density, repeated Rule-30 rows, LZ78 phrases). One deterministic stream per length, drawn identically by math-deep.js.""" state = (0x00133701 + n * 7919) & 0xFFFFFFFF ts, ops, prs, dns, rps, lzs = [], [], [], [], [], [] for _ in range(r): k, state = rand_int(state, n) t, op, pr = collatz_stats(k) ts.append(t) ops.append(op) prs.append(pr) dn, rp = rule30_stats(k, n, n) dns.append(dn) rps.append(rp) lzs.append(lz78(bits_msb(k, n))) def ms(a): return nsum(a) / r, sd(a) return {"T": ms(ts), "odd": ms(ops), "pr": ms(prs), "dens": ms(dns), "rep": ms(rps), "lz": ms(lzs)} def lz_exact_null(n): """Exact LZ78 null for one bit length: (mean, sd, support, counts, total). Enumerates all 2^(n-1) keys of length n. Used for n <= EXACT_LZ_N, where the phrase-count law is discrete and left-skewed enough that a normal approximation to a Monte-Carlo z misstates the tail by an order of magnitude or more (checked against the enumeration for n = 16 and 20). """ top = 1 << (n - 1) total = top hist = {} for v in range(top, top << 1): ph = lz78(bits_msb(v, n)) hist[ph] = hist.get(ph, 0) + 1 mean = 0.0 for ph, c in hist.items(): mean += ph * c mean /= total var = 0.0 for ph, c in hist.items(): d = ph - mean var += c * d * d var /= total support = sorted(hist) return mean, math.sqrt(var), support, [hist[k] for k in support], total def lz_exact_p(null, k): """Exact P(X <= k) from an lz_exact_null() table (lower tail).""" mean, sig, support, counts, total = null s = 0 for i in range(len(support)): if support[i] > k: break s += counts[i] return s / total def deep_null_stream(ns, r): """Monte-Carlo null for the aggregated bitstream's LZ78 phrase count. Draws whole keys -- one uniform MSB-set integer per observed bit length -- and concatenates them in dataset order. The null therefore inherits the real stream's key boundaries and length profile; a flat i.i.d. bit stream would not, and would place the expected phrase count in the wrong place. Returns (mean, sd, draws). """ state = (0x0BADC0DE + sum(ns)) & 0xFFFFFFFF vals = [] for _ in range(r): bits = [] for n in ns: k, state = rand_int(state, n) bits.extend(bits_msb(k, n)) vals.append(lz78(bits)) return nsum(vals) / r, sd(vals), vals def deep_z(v, mu, sig): """Standardised score with a degenerate-null guard (sd 0 -> z 0).""" if sig <= 0.0: return 0.0 return (v - mu) / sig def deep(solved, out=print): """Checks 11-15 of the deep cryptographic / information-theoretic protocol. Same contract as advanced() and lowlevel(): every number is measured from the solved keys, the Monte-Carlo nulls use one deterministic xorshift32 stream (identical seeds and draws in math-deep.js), and each check carries a family-wise alpha of 0.01 with Bonferroni correction over the candidates it scans. Returns the measured results. Three null models are used because the statistics differ in kind. Check 11 resamples whole datasets: its pairwise distances share keys, so the analytic binomial bands and the pooled z are not independent and their nominal p-values are too small. Checks 12 and 15 keep their Monte-Carlo draws (key-structured for the pooled bitstream, length-matched for the secp256k1 field image) and read the null's own spread rather than a uniform law. Where a null is discrete and small-n -- the LZ78 phrase count and the Collatz / Rule-30 statistics -- the Monte-Carlo resolution is coarser than the Bonferroni threshold, so per-key z's are reported descriptively and the inference is carried by aggregate tests, except for LZ78 below EXACT_LZ_N where the null is enumerated in full. """ N = len(solved) flags = [] res = {"N": N} rows = sorted(solved, key=lambda r: r["n"]) keys = [r["key"] for r in rows] ns = [r["n"] for r in rows] res["ns"] = ns blist = [bits_msb(k, n) for k, n in zip(keys, ns)] # shared per-key Monte-Carlo nulls (used by checks 12 and 14) nulls = [deep_null_len(n, R_LEN) for n in ns] # ================= check 11: pairwise Hamming distance matrix ========== out("") out("-" * 74) out("DEEP CHECK 11 — PAIRWISE BINARY HAMMING DISTANCE MATRIX") out("-" * 74) n_pairs = N * (N - 1) // 2 q1, q2, q3 = (norm_two_sided_p(1.0), norm_two_sided_p(2.0), norm_two_sided_p(3.0)) band_exp = [(1.0 - q1) * n_pairs, (q1 - q2) * n_pairs, (q2 - q3) * n_pairs, q3 * n_pairs] def pair_stats(ks): """The five matrix statistics of a key set -- one code path shared by the observed keys and every whole-dataset null draw, so both sides are computed identically.""" zv, msum, mcnt = [], [0.0] * N, [0] * N for i in range(N): ki = ks[i] for j in range(i + 1, N): h = (ki ^ ks[j]).bit_count() if ns[i] == ns[j]: free, det = ns[i], 0 else: free, det = ns[j] - 1, 1 z = (h - (free / 2.0 + det)) / math.sqrt(free / 4.0) zv.append(z) msum[i] += z mcnt[i] += 1 msum[j] += z mcnt[j] += 1 bands = [0, 0, 0, 0] for z in zv: a = abs(z) bands[0 if a <= 1.0 else (1 if a <= 2.0 else (2 if a <= 3.0 else 3))] += 1 chi = 0.0 for b in range(4): d = bands[b] - band_exp[b] chi += (d * d) / band_exp[b] zmax = 0.0 for z in zv: if abs(z) > zmax: zmax = abs(z) mmax = 0.0 for t in range(N): v = abs(msum[t] / mcnt[t]) if v > mmax: mmax = v return (nsum(zv) / n_pairs, sd(zv), zmax, mmax, chi, bands, zv, msum, mcnt) (zmean, zsd, zmax, mmax, chi_band, band_obs, zs, msum, mcnt) = pair_stats(keys) mkey = [msum[t] / mcnt[t] for t in range(N)] zpool = zmean * math.sqrt(n_pairs) meta = [] for i in range(N): for j in range(i + 1, N): h = (keys[i] ^ keys[j]).bit_count() free, det = (ns[i], 0) if ns[i] == ns[j] else (ns[j] - 1, 1) meta.append((i, j, h, free / 2.0 + det)) # Pairs share keys, so the analytic chi2(3) and pooled-z variances # understate the real spread: two pairs sharing a key co-vary on every # position above max(n_j, n_k). Re-measure all five statistics on whole # random datasets with the observed length profile instead. st11 = 0x11C0DE11 sim11 = [] for _ in range(R_DATA): ks = [] for n in ns: k, st11 = rand_int(st11, n) ks.append(k) sim11.append(pair_stats(ks)[:5]) sim_mu, sim_sd, sim_z, sim_p = [], [], [], [] for b in range(5): arr = [s[b] for s in sim11] mu = nsum(arr) / R_DATA sig = sd(arr) z = deep_z((zmean, zsd, zmax, mmax, chi_band)[b], mu, sig) sim_mu.append(mu) sim_sd.append(sig) sim_z.append(z) sim_p.append(norm_two_sided_p(z)) p_naive_band = chi2sf(chi_band, 3) zc_pair = z_for_p(0.01 / n_pairs) zc_key = z_for_p(0.01 / N) n_pair_flags = 0 i_max = 0 for t in range(n_pairs): if abs(zs[t]) > zc_pair: n_pair_flags += 1 if abs(zs[t]) > abs(zs[i_max]): i_max = t k_max = 0 n_key_flags = 0 for t in range(N): if abs(mkey[t]) > zc_key: n_key_flags += 1 if abs(mkey[t]) > abs(mkey[k_max]): k_max = t a11on = (sim_p[0] < 0.01 or sim_p[1] < 0.01 or sim_p[2] < 0.01 or sim_p[3] < 0.01 or sim_p[4] < 0.01 or n_pair_flags > 0 or n_key_flags > 0) i_a, i_b, h_a, mu_a = meta[i_max] out(f"pairs compared : {n_pairs} " f"(N={N} -> N(N-1)/2 combinations)") out(f"mean pairwise z (null 0) : {sg(zmean, 4)} " f"(null {sg(sim_mu[0], 4)} +- {f(sim_sd[0], 4)} from {R_DATA} " f"draws -> p {f(sim_p[0], 4)})") out(f"sd of the pairwise z (null 1) : {f(zsd, 4)} " f"(null {f(sim_mu[1], 4)} +- {f(sim_sd[1], 4)} -> " f"p {f(sim_p[1], 4)})") out(f"band counts |z|<=1 /2 /3 />3 : {band_obs[0]} / " f"{band_obs[1]} / {band_obs[2]} / {band_obs[3]} (binomial " f"expectation {band_exp[0]:.1f} / {band_exp[1]:.1f} / " f"{band_exp[2]:.1f} / {band_exp[3]:.1f})") out(f"band chi2(3) : {f(chi_band, 2)} " f"(null {f(sim_mu[4], 2)} +- {f(sim_sd[4], 2)} -> " f"p {f(sim_p[4], 4)}; the naive chi2 p {f(p_naive_band, 4)} " f"assumes independent pairs and is not used)") out(f"largest pairwise |z| : {sg(abs(zs[i_max]), 3)} " f"(n {ns[i_a]} vs {ns[i_b]}, H {h_a} vs {mu_a:.1f} expected) " f"(null {f(sim_mu[2], 3)} +- {f(sim_sd[2], 3)} -> " f"p {f(sim_p[2], 4)})") out(f"Bonferroni |z| threshold : {f(zc_pair, 3)} " f"({n_pair_flags} of {n_pairs} pairs beyond)") out(f"key-level mean distance z: max |m| : " f"{sg(abs(mkey[k_max]), 3)} at n={ns[k_max]} (threshold " f"{f(zc_key, 3)}, {n_key_flags} of {N} keys beyond) (null max " f"{f(sim_mu[3], 3)} +- {f(sim_sd[3], 3)} -> p {f(sim_p[3], 4)})") if a11on: flags.append("D11 pairwise Hamming") res["d11"] = {"nPairs": n_pairs, "zMean": zmean, "zSd": zsd, "zPool": zpool, "bandObs": band_obs, "bandExp": band_exp, "chi": chi_band, "pNaive": p_naive_band, "simMu": sim_mu, "simSd": sim_sd, "simZ": sim_z, "simP": sim_p, "rData": R_DATA, "zcPair": zc_pair, "nPairFlags": n_pair_flags, "iMax": i_max, "zMax": zs[i_max], "zcKey": zc_key, "nKeyFlags": n_key_flags, "kMax": k_max, "mMax": mkey[k_max], "zs": zs, "mkey": mkey, "top": sorted(range(n_pairs), key=lambda t: -abs(zs[t]))[:8], "meta": meta, "flag": a11on} # ================= check 12: binary runs & LZ78 complexity ============ out("") out("-" * 74) out("DEEP CHECK 12 — BINARY RUNS (NIST SP 800-22) & LEmpel-Ziv PROFILE") out("-" * 74) run_keys = [t for t in range(N) if ns[t] >= RUNS_MIN_N] run_crit = 0.01 / len(run_keys) lz_key_crit = 0.01 / N run_min_p, run_min_n = 1.0, -1 run_flags = 0 lz_key_flags = 0 rep_flags = 0 run_rows = [] lz_rows = [] lz_cache = {} lz_dev = 0.0 lz_var_sum = 0.0 lz_min_p, lz_min_n = 1.0, -1 lz_max_z, lz_max_zt = 0.0, 0 lz_exact_used = 0 for t in run_keys: ok, pi, V, p = nist_runs(blist[t]) run_rows.append((t, pi, V, p, ok)) if p < run_min_p: run_min_p, run_min_n = p, t if p < run_crit: run_flags += 1 for t in range(N): n = ns[t] lz = lz78(blist[t]) if n <= EXACT_LZ_N: if n not in lz_cache: lz_cache[n] = lz_exact_null(n) mu, sig = lz_cache[n][0], lz_cache[n][1] p_lz = lz_exact_p(lz_cache[n], lz) lz_exact_used += 1 if p_lz < lz_key_crit: lz_key_flags += 1 if p_lz < lz_min_p: lz_min_p, lz_min_n = p_lz, t else: mu, sig = nulls[t]["lz"] p_lz = None zlz = deep_z(lz, mu, sig) if abs(zlz) > lz_max_z: lz_max_z, lz_max_zt = abs(zlz), t lz_dev += lz - mu lz_var_sum += sig * sig c, e = window_collisions(blist[t], 8) rp = 1.0 if e <= 0.0 else poisson_tail(c, e) lz_rows.append((t, lz, mu, zlz, p_lz, c, e, rp)) if e > 0.0 and rp < 0.01 / N: rep_flags += 1 z_lz_agg = lz_dev / math.sqrt(lz_var_sum) if lz_var_sum > 0.0 else 0.0 p_lz_agg = norm_two_sided_p(z_lz_agg) stream = [] for t in range(N): stream.extend(blist[t]) s_ok, s_pi, s_V, s_p = nist_runs(stream) s_lz = lz78(stream) s_mu, s_sd, s_draws = deep_null_stream(ns, R_POOL) s_zlz = deep_z(s_lz, s_mu, s_sd) s_n_le = 0 for v in s_draws: if v <= s_lz: s_n_le += 1 s_lz_p = (1.0 + s_n_le) / (R_POOL + 1.0) s_lz_pnorm = norm_two_sided_p(s_zlz) s_c, s_e = window_collisions(stream, 8) s_prep = 1.0 if s_e <= 0.0 else poisson_tail(s_c, s_e) a12on = (s_p < 0.01 or run_flags > 0 or s_lz_p < 0.01 or p_lz_agg < 0.01 or lz_key_flags > 0 or s_prep < 0.01 or rep_flags > 0) out(f"NIST runs test applied to keys : n >= {RUNS_MIN_N} bits " f"({len(run_keys)} of {N}; shorter keys are degenerate for it)") out(f"pooled bitstream runs (n={len(stream)}) : pi {f(s_pi, 4)}, " f"V {s_V}, p {f(s_p, 4)}") out(f"weakest per-key runs p : {f(run_min_p, 4)} at " f"n={ns[run_min_n]} (Bonferroni {run_crit:.5f}, {run_flags} beyond)") out(f"LZ78 phrases, pooled stream : {s_lz} vs " f"{s_mu:.1f} +- {s_sd:.1f} null ({R_POOL} key-structured draws) -> " f"z {sg(s_zlz, 3)}; {s_n_le} of {R_POOL} draws at or below, " f"empirical p {f(s_lz_p, 4)} (normal p {f(s_lz_pnorm, 6)})") out(f"per-key LZ78, aggregate : {sg(z_lz_agg, 3)} sd " f"over {N} keys -> p {f(p_lz_agg, 4)} (largest |z| " f"{f(lz_max_z, 3)} at n={ns[lz_max_zt]}; {lz_exact_used} keys with " f"an exact null, {lz_key_flags} beyond Bonferroni " f"{lz_key_crit:.5f}, smallest exact p {f(lz_min_p, 4)} at " f"n={ns[lz_min_n]})") out(f"repeated 8-bit windows, pooled : {s_c} vs {s_e:.2f} " f"expected (Poisson P {f(s_prep, 4)})") out(f"keys with excess 8-bit repeats : {rep_flags} of {N} " f"(Bonferroni {0.01 / N:.5f})") if a12on: flags.append("D12 runs / LZ complexity") res["d12"] = {"runKeys": len(run_keys), "runMinP": run_min_p, "runMinN": run_min_n, "runFlags": run_flags, "runCrit": run_crit, "lzKeyCrit": lz_key_crit, "streamN": len(stream), "sPi": s_pi, "sV": s_V, "sP": s_p, "sOk": s_ok, "sLz": s_lz, "sLzMu": s_mu, "sLzSd": s_sd, "sLzZ": s_zlz, "sLzP": s_lz_p, "sLzPNorm": s_lz_pnorm, "sLzNLe": s_n_le, "rPool": R_POOL, "lzAggZ": z_lz_agg, "lzAggP": p_lz_agg, "lzMaxZ": lz_max_z, "lzMaxZT": lz_max_zt, "lzExact": lz_exact_used, "lzKeyFlags": lz_key_flags, "lzMinP": lz_min_p, "lzMinN": lz_min_n, "sRepC": s_c, "sRepE": s_e, "sRepP": s_prep, "repFlags": rep_flags, "runRows": run_rows, "lzRows": lz_rows, "flag": a12on} # ================= check 13: bit-segment mutual information =========== out("") out("-" * 74) out("DEEP CHECK 13 — BIT-SEGMENT MUTUAL INFORMATION (NON-LINEAR LINKS)") out("-" * 74) elig = [t for t in range(N) if ns[t] >= 32] lsb = {t: [(keys[t] >> d) & 1 for d in range(32)] for t in elig} msb = {t: [(keys[t] >> (ns[t] - 1 - d)) & 1 for d in range(32)] for t in elig} ln2 = math.log(2.0) mi_pairs = [] def push_mi(fam, a, b, ga, gb, excl): tbl = [[0, 0], [0, 0]] n = 0 for t in elig: if excl is not None and excl(t, a, b): continue tbl[ga[t][a]][gb[t][b]] += 1 n += 1 if n < 4: return mi = mi_bits(tbl, n) g2 = 2.0 * n * mi * ln2 mi_pairs.append({"f": fam, "a": a, "b": b, "n": n, "mi": mi, "g2": g2, "p": chi2sf(g2, 1)}) for a in range(32): for b in range(a + 1, 32): push_mi("L", a, b, lsb, lsb, None) for a in range(1, 32): for b in range(a + 1, 32): push_mi("M", a, b, msb, msb, None) for a in range(1, 32): for b in range(32): push_mi("X", a, b, msb, lsb, lambda t, aa, bb: ns[t] == aa + bb + 1) n_mi = len(mi_pairs) for r in mi_pairs: r["ic"] = r["mi"] - 1.0 / (2.0 * r["n"] * ln2) mi_crit = 0.01 / n_mi n_mi_flags = 0 n_mi_05 = 0 g2_sum = 0.0 for t in range(n_mi): g2_sum += mi_pairs[t]["g2"] if mi_pairs[t]["p"] < mi_crit: n_mi_flags += 1 if mi_pairs[t]["p"] < 0.05: n_mi_05 += 1 j_best = 0 for t in range(1, n_mi): if mi_pairs[t]["ic"] > mi_pairs[j_best]["ic"]: j_best = t top_mi = sorted(range(n_mi), key=lambda t: -mi_pairs[t]["g2"])[:10] # block-level MI: quartiles of the MSB-32 and LSB-32 popcounts bx = [(keys[t] >> (ns[t] - 32)).bit_count() for t in elig] by = [(keys[t] & 0xFFFFFFFF).bit_count() for t in elig] sx, sy = sorted(bx), sorted(by) cuts_x = [sx[k * len(sx) // 4 - 1] for k in (1, 2, 3)] cuts_y = [sy[k * len(sy) // 4 - 1] for k in (1, 2, 3)] def quart(v, cuts): q = 0 for c in cuts: if v > c: q += 1 return q def block_mi(xs, ys): tbl = [[0] * 4 for _ in range(4)] for u, v in zip(xs, ys): tbl[quart(u, cuts_x)][quart(v, cuts_y)] += 1 return mi_bits(tbl, len(xs)) mi_block = block_mi(bx, by) state = (0x1CEB00DA + len(elig)) & 0xFFFFFFFF sims = [] for _ in range(R_BLOCK): gx, gy = [], [] for t in elig: k, state = rand_int(state, ns[t]) gx.append((k >> (ns[t] - 32)).bit_count()) gy.append((k & 0xFFFFFFFF).bit_count()) sims.append(block_mi(gx, gy)) n_ge = 0 for v in sims: if v >= mi_block: n_ge += 1 p_block = (1 + n_ge) / (R_BLOCK + 1.0) block_ic = mi_block - 9.0 / (2.0 * len(elig) * ln2) def ent_pos(vals): c1 = 0 for v in vals: c1 += v c0 = len(vals) - c1 h = 0.0 for c in (c0, c1): if c: p = c / len(vals) h -= p * math.log2(p) m = (1 if c0 else 0) + (1 if c1 else 0) return h + (m - 1) / (2.0 * len(vals) * ln2) ent_l = [ent_pos([lsb[t][d] for t in elig]) for d in range(32)] ent_m = [ent_pos([msb[t][d] for t in elig]) for d in range(1, 32)] a13on = n_mi_flags > 0 or p_block < 0.01 best = mi_pairs[j_best] cnt_l = sum(1 for r in mi_pairs if r["f"] == "L") cnt_m = sum(1 for r in mi_pairs if r["f"] == "M") out(f"keys with both 32-bit segments : {len(elig)} (n >= 32; " f"MSB-aligned d=0 is the fixed key MSB)") out(f"bit-position pairs tested : {n_mi} (LSB-LSB {cnt_l}, " f"MSB-MSB {cnt_m}, cross {n_mi - cnt_l - cnt_m})") out(f"largest bias-corrected MI : {sg(best['ic'], 4)} bits " f"at {best['f']}({best['a']},{best['b']}) — G2 {f(best['g2'], 2)}, " f"p {f(best['p'], 4)}") out(f"pairs beyond Bonferroni : {n_mi_flags} of {n_mi} " f"(alpha {mi_crit:.6f}); {n_mi_05} with p < 0.05 vs " f"{0.05 * n_mi:.1f} expected") out(f"sum of G2 over pairs : {f(g2_sum, 1)} vs {n_mi} " f"expected (descriptive: pairs share keys)") out(f"block MI (quartile MSB32 x LSB32) : {f(mi_block, 4)} bits " f"(corrected {f(block_ic, 4)}), null {f(mean(sims), 4)} from " f"{R_BLOCK} draws -> p {f(p_block, 4)} (upper tail: excess " f"dependence)") out(f"segment entropy (corrected) : LSB mean " f"{f(nsum(ent_l) / len(ent_l), 4)} (min {f(min(ent_l), 4)}), MSB " f"mean {f(nsum(ent_m) / len(ent_m), 4)} (min {f(min(ent_m), 4)}), " f"fair bit = 1") if a13on: flags.append("D13 segment mutual information") res["d13"] = {"nPairs": n_mi, "nFlags": n_mi_flags, "nP05": n_mi_05, "g2Sum": g2_sum, "jBest": j_best, "miCrit": mi_crit, "top": top_mi, "pairs": mi_pairs, "block": mi_block, "blockIc": block_ic, "blockMu": mean(sims), "blockP": p_block, "entL": ent_l, "entM": ent_m, "nElig": len(elig), "flag": a13on} # ================= check 14: Collatz / Rule-30 dynamics =============== out("") out("-" * 74) out("DEEP CHECK 14 — CHAOTIC TRAJECTORY & STATE-SPACE DYNAMICS") out("-" * 74) zt, zo, zp, zd, reps, dens = [], [], [], [], [], [] for t in range(N): nul = nulls[t] ts, op, pr = collatz_stats(keys[t]) dn, rp = rule30_stats(keys[t], ns[t], ns[t]) zt.append(deep_z(ts, *nul["T"])) zo.append(deep_z(op, *nul["odd"])) zp.append(deep_z(pr, *nul["pr"])) zd.append(deep_z(dn, *nul["dens"])) dens.append(dn) reps.append(rp) zc14 = z_for_p(0.01 / (4 * N)) n_zt = n_zo = n_zp = n_zd = 0 for t in range(N): if abs(zt[t]) > zc14: n_zt += 1 if abs(zo[t]) > zc14: n_zo += 1 if abs(zp[t]) > zc14: n_zp += 1 if abs(zd[t]) > zc14: n_zd += 1 i_t = i_o = i_p = i_d = 0 for t in range(1, N): if abs(zt[t]) > abs(zt[i_t]): i_t = t if abs(zo[t]) > abs(zo[i_o]): i_o = t if abs(zp[t]) > abs(zp[i_p]): i_p = t if abs(zd[t]) > abs(zd[i_d]): i_d = t # A R_LEN-draw null resolves tail probabilities no finer than about # 1/(R_LEN+1), far coarser than the Bonferroni threshold 0.01/(4N) ~ 3e-5, # and the small-n null laws are discrete and skewed enough that a normal z # overstates their tails. So the per-key z's are reported descriptively # and the inference is carried by aggregate tests: keys are independent, # so the mean of the N standardised z's has sd 1/sqrt(N), and the count of # |z| > 2 is Binomial(N, P(|z|>2)) with P = 0.0455. fams = (("T", "Collatz stopping time", zt), ("odd", "odd-step share", zo), ("pr", "parity-sequence runs", zp), ("dens", "Rule-30 1-density", zd)) ixmap = {"T": i_t, "odd": i_o, "pr": i_p, "dens": i_d} p_z2 = norm_two_sided_p(2.0) agg = {} for k14, _lbl, arr in fams: mz = nsum(arr) / N zagg = mz * math.sqrt(N) c2 = 0 for v in arr: if abs(v) > 2.0: c2 += 1 agg[k14] = (mz, zagg, norm_two_sided_p(zagg), c2, binom_tail(c2, N, p_z2)) rep_obs = 0 for v in reps: rep_obs += v rep_exp = 0.0 for t in range(N): rep_exp += nulls[t]["rep"][0] rep_p = 1.0 if rep_exp <= 0.0 else poisson_tail(rep_obs, rep_exp) idx14 = list(range(N)) idx14.sort(key=lambda t: -max(abs(zt[t]), abs(zo[t]), abs(zp[t]), abs(zd[t]))) i14 = idx14[0] a14on = (agg["T"][2] < 0.01 or agg["odd"][2] < 0.01 or agg["pr"][2] < 0.01 or agg["dens"][2] < 0.01 or agg["T"][4] < 0.01 or agg["odd"][4] < 0.01 or agg["pr"][4] < 0.01 or agg["dens"][4] < 0.01 or rep_p < 0.01) out(f"null model : {R_LEN} xorshift32 " f"draws of each key's bit length (Collatz, parity, Rule-30, LZ)") for k14, lbl, arr in fams: mz, zagg, pagg, c2, pcnt = agg[k14] tmax = ixmap[k14] out(f"{lbl + ': max |z|':<38}: {sg(abs(arr[tmax]), 3)} at " f"n={ns[tmax]} (aggregate mean z {sg(mz, 4)} -> p {f(pagg, 4)}; " f"{c2} of {N} beyond |z|=2 vs {N * p_z2:.2f} expected, " f"p {f(pcnt, 4)})") out(f"per-key |z| beyond Bonferroni {zc14:.3f} : {n_zt} / {n_zo} / " f"{n_zp} / {n_zd} of {N} (descriptive: a {R_LEN}-draw null " f"resolves only ~1/{R_LEN + 1}, coarser than {0.01 / (4 * N):.1e})") out(f"repeated Rule-30 rows across keys : {rep_obs} vs " f"{rep_exp:.2f} expected (Poisson P {f(rep_p, 4)})") if a14on: flags.append("D14 Collatz / Rule 30") res["d14"] = {"zc": zc14, "nZt": n_zt, "nZo": n_zo, "nZp": n_zp, "nZd": n_zd, "iT": i_t, "iO": i_o, "iP": i_p, "iD": i_d, "repObs": rep_obs, "repExp": rep_exp, "repP": rep_p, "pZ2": p_z2, "agg": {k14: {"meanZ": v[0], "z": v[1], "p": v[2], "c2": v[3], "pCnt": v[4]} for k14, v in agg.items()}, "top": idx14[:10], "i14": i14, "zt": zt, "zo": zo, "zp": zp, "zd": zd, "reps": reps, "dens": dens, "flag": a14on} # ================= check 15: secp256k1 geometric residue ============= out("") out("-" * 74) out("DEEP CHECK 15 — secp256k1 ELLIPTIC-CURVE FIELD RESIDUE GEOMETRY") out("-" * 74) resid = 0 zero_a = 0 bad_keys = 0 root_ok = 0 us = [] ys = [] bin16 = [0] * 16 for k in keys: if k == 0 or k >= SECP_ORDER: bad_keys += 1 a = (pow(k, 3, P256) + 7) % P256 bin16[(a * 16) // P256] += 1 us.append(a / P256) leg = pow(a, (P256 - 1) // 2, P256) if leg == 1: resid += 1 y = pow(a, (P256 + 1) // 4, P256) if y * y % P256 == a: root_ok += 1 ys.append(y / P256) elif leg == 0: zero_a += 1 p_leg = binom_two_sided(resid, N) chi16 = 0.0 exp16 = N / 16.0 for b in range(16): d = bin16[b] - exp16 chi16 += (d * d) / exp16 mean_u = nsum(us) / N mean_y = nsum(ys) / len(ys) if ys else 0.0 # (K^3+7) mod p is uniform on [0, p) only for keys that fill the field. # While 3n < 256 no modular reduction ever happens, so both the 16-bin # image and mean u inherit the key-length profile; their null is a # simulation of the same lengths, not U[0, 1). st15 = 0x5EC12345 sim_chi, sim_u = [], [] for _ in range(R_IMG): b16 = [0] * 16 u2 = [] for n in ns: k, st15 = rand_int(st15, n) a = (pow(k, 3, P256) + 7) % P256 b16[(a * 16) // P256] += 1 u2.append(a / P256) c = 0.0 for b in range(16): d = b16[b] - exp16 c += (d * d) / exp16 sim_chi.append(c) sim_u.append(nsum(u2) / N) chi_mu = nsum(sim_chi) / R_IMG chi_sd = sd(sim_chi) u_mu = nsum(sim_u) / R_IMG u_sd = sd(sim_u) z_chi = deep_z(chi16, chi_mu, chi_sd) p_chi = norm_two_sided_p(z_chi) z_u = deep_z(mean_u, u_mu, u_sd) p_u = norm_two_sided_p(z_u) # A modular square root of a residue is uniform on [0, p), so mean y/p # keeps its analytic U[0, 1) null. z_y = (mean_y - 0.5) / math.sqrt(1.0 / (12.0 * len(ys))) if ys else 0.0 p_y = norm_two_sided_p(z_y) alpha15 = 0.01 / 4 a15on = (p_leg < alpha15 or p_chi < alpha15 or p_u < alpha15 or (ys and p_y < alpha15)) out(f"keys outside [1, n-1] or zero : {bad_keys} of {N} " f"(group order n = 2^256 - 432420386565659656852420866394967595805)") out(f"quadratic residues of K^3+7 (mod p) : {resid} of {N} vs " f"{N / 2:.1f} expected -> exact two-sided binomial p " f"{f(p_leg, 4)}") out(f"K^3+7 == 0 (mod p) : {zero_a}; square roots " f"verified {root_ok}/{resid}") out(f"field image uniformity, 16 bins : chi2(15) {f(chi16, 2)} vs " f"length-matched null {f(chi_mu, 2)} +- {f(chi_sd, 2)} ({R_IMG} " f"draws) -> z {sg(z_chi, 3)}, p {f(p_chi, 4)} (alpha " f"{alpha15:.4f})") out(f"mean (K^3+7)/p : {f(mean_u, 4)} vs null " f"{f(u_mu, 4)} +- {f(u_sd, 4)} -> z {sg(z_u, 3)}, p {f(p_u, 4)} " f"(0.5 is not the null while 3n < 256: no reduction happens)") out(f"mean y/p over residues : {f(mean_y, 4)} vs 0.5 -> " f"z {sg(z_y, 3)}, p {f(p_y, 4)}") if a15on: flags.append("D15 secp256k1 residue") res["d15"] = {"badKeys": bad_keys, "resid": resid, "pLeg": p_leg, "zeroA": zero_a, "rootOk": root_ok, "chi": chi16, "chiMu": chi_mu, "chiSd": chi_sd, "zChi": z_chi, "pChi": p_chi, "bins": bin16, "meanU": mean_u, "uMu": u_mu, "uSd": u_sd, "zU": z_u, "pU": p_u, "meanY": mean_y, "zY": z_y, "pY": p_y, "nY": len(ys), "rImg": R_IMG, "alpha": alpha15, "flag": a15on} out("") out("=" * 74) out("DEEP VERDICT (family-wise alpha = 0.01 per check)") out("=" * 74) if flags: out(f"FLAGS: {len(flags)} of 5 checks flagged: {', '.join(flags)}") out("At least one deep statistic exceeded its null expectation;") out("inspect the flagged row before asserting a generator artifact.") else: out("FLAGS: 0 of 5 checks flagged.") out("No pairwise clustering, no compressible bitstream structure, no") out("non-linear segment dependence, no low-entropy chaotic") out("trajectory and no secp256k1 residue anomaly was uncovered: the") out("3403 pairwise Hamming distances are binomial, the runs and") out("LZ78 complexity match random bitstreams of the same length, the") out("bit-segment mutual informations sit inside their Bonferroni") out("band, the Collatz and Rule-30 statistics match their simulated") out("nulls, and K^3+7 mod p is a quadratic residue exactly as often") out("as chance predicts. The null hypothesis (true cryptographic") out("pseudo-randomness) holds across all five deep metrics — no") out("information-theoretic compression, pairwise clustering,") out("non-linear dependency or deterministic state machine is") out("detected.") res["flags"] = flags return res def main(argv=None): ap = argparse.ArgumentParser( description="Forensic math analysis of the Bitcoin Challenge dataset.") ap.add_argument("--data", help="path to BTC-Solved-Unsolved.txt or the JSON") ap.add_argument("--predict", type=int, default=71, metavar="N", help="unsolved puzzle to report a sub-space for (default 71)") args = ap.parse_args(argv) rows, source = load_dataset(args.data) print(f"dataset: {source}") solved, slope_p, intercept_p, resid_p = analyze(rows) target = next((r for r in rows if r["n"] == args.predict and not r["solved"]), None) if target is None: print(f"\n(puzzle #{args.predict} is not in the dataset as unsolved; " f"showing its range only)") predict(solved, args.predict, slope_p, intercept_p, resid_p) advanced(solved) lowlevel(solved) deep(solved) return 0 if __name__ == "__main__": sys.exit(main())