"""Homophonic / verbose-cipher solver (v2): 4-gram character LM with interpolated smoothing stored as a dense table, simulated annealing over many-to-one keys with incremental vectorised scoring (only the windows touched by the changed unit are rescored), frequency-constraint penalty, random restarts. Validated on planted keys (t_homophonic_control.py) before use on the manuscript. """ import math, random, collections import numpy as np class LM: """Interpolated 4-gram over characters incl. space; alphabet = letters with >= min_count.""" def __init__(self, words, keep_spaces=True, min_count=20, k=(0.5, 0.5, 0.5)): seq = [] for w in words: seq += list(w) + ([" "] if keep_spaces else []) cnt = collections.Counter(seq) self.alpha = [a for a, c in cnt.most_common() if c >= min_count and a != " "] if keep_spaces: self.alpha.append(" ") self.L = L = len(self.alpha); self.idx = {a: i for i, a in enumerate(self.alpha)} self.keep_spaces = keep_spaces s = np.array([self.idx[c] for c in seq if c in self.idx], dtype=np.int64) self.seq = s u = np.bincount(s, minlength=L).astype(float); P1 = (u + 1) / (u.sum() + L) c2 = np.zeros((L, L)); np.add.at(c2, (s[:-1], s[1:]), 1) P2 = (c2 + k[0] * P1[None, :]) / (c2.sum(1, keepdims=True) + k[0]) c3 = np.zeros((L, L, L)); np.add.at(c3, (s[:-2], s[1:-1], s[2:]), 1) P3 = (c3 + k[1] * P2[None, :, :]) / (c3.sum(2, keepdims=True) + k[1]) c4 = np.zeros((L, L, L, L)); np.add.at(c4, (s[:-3], s[1:-2], s[2:-1], s[3:]), 1) P4 = (c4 + k[2] * P3[None, :, :, :]) / (c4.sum(3, keepdims=True) + k[2]) self.logP4 = np.log(P4).reshape(-1) # index = ((a*L+b)*L+c)*L+d self.logP2 = np.log(P2) self.uni = u / u.sum() self.vocab = set(words) self.sp = self.idx.get(" ", None) def score_seq(self, s): """Mean 4-gram log-prob of an int sequence (fixed letters).""" if len(s) < 4: return 0.0 L = self.L; idx = ((s[:-3] * L + s[1:-2]) * L + s[2:-1]) * L + s[3:] return float(self.logP4[idx].mean()) def own_score(self): return self.score_seq(self.seq) class Cipher: def __init__(self, units, min_count=5, keep_spaces=True): cnt = collections.Counter(u for u in units if u != " ") self.alpha = [u for u, c in cnt.most_common() if c >= min_count] self.idx = {u: i for i, u in enumerate(self.alpha)} self.keep_spaces = keep_spaces seq = [] for u in units: if u == " ": if keep_spaces: seq.append(-1) elif u in self.idx: seq.append(self.idx[u]) # rare units are dropped (treated as absent) self.seq = np.array(seq, dtype=np.int64) self.nU = len(self.alpha) self.pos = [np.flatnonzero(self.seq == i) for i in range(self.nU)] self.freq = np.array([cnt[u] for u in self.alpha], float) def solve(ciph, lm, lam=0.3, restarts=4, iters=200000, T0=None, seed=0, verbose=False, polish_sweeps=6): rng = random.Random(seed) L = lm.L; sp = lm.sp letters = [i for i in range(L) if i != sp] uni = lm.uni[letters] / lm.uni[letters].sum() freq = ciph.freq / ciph.freq.sum() seq = ciph.seq; T = len(seq); nU = ciph.nU # window start indices affected by each unit (unique), for 4-grams starting at p-3..p starts = [] for i in range(nU): p = ciph.pos[i] w = np.concatenate([p - 3, p - 2, p - 1, p]); w = w[(w >= 0) & (w <= T - 4)] starts.append(np.unique(w)) logP4 = lm.logP4 def window_score(P, st): idx = ((P[st] * L + P[st + 1]) * L + P[st + 2]) * L + P[st + 3] return logP4[idx].sum() def kl_of(m): f = np.zeros(L); np.add.at(f, m, freq); f = f[letters] + 1e-9; f /= f.sum() return float((f * np.log(f / (uni + 1e-9))).sum()) T0 = T0 or 0.02 * T # temperature on the total-log-likelihood scale best = (-1e18, None, None, None) for r in range(restarts): order = np.argsort(-ciph.freq); lord = [letters[i] for i in np.argsort(-uni)] m = np.zeros(nU, int) for k, ui in enumerate(order): m[ui] = lord[k % len(lord)] if r > 0: for i in range(nU): if rng.random() < 0.5: m[i] = rng.choice(letters) P = np.where(seq >= 0, m[np.maximum(seq, 0)], sp if sp is not None else 0).astype(np.int64) all_idx = ((P[:-3] * L + P[1:-2]) * L + P[2:-1]) * L + P[3:] total = float(logP4[all_idx].sum()); kl = kl_of(m); score = total - lam * T * kl cur_T = T0 for it in range(iters): i = rng.randrange(nU) if rng.random() < 0.6: new = rng.choice(letters) if new == m[i]: continue st = starts[i]; old_s = window_score(P, st) P[ciph.pos[i]] = new; new_s = window_score(P, st) m2 = m.copy(); m2[i] = new; kl2 = kl_of(m2) score2 = score + (new_s - old_s) - lam * T * (kl2 - kl) if score2 >= score or rng.random() < math.exp((score2 - score) / cur_T): m, kl, score = m2, kl2, score2 else: P[ciph.pos[i]] = m[i] else: j = rng.randrange(nU) if m[i] == m[j]: continue st = np.union1d(starts[i], starts[j]); old_s = window_score(P, st) P[ciph.pos[i]] = m[j]; P[ciph.pos[j]] = m[i]; new_s = window_score(P, st) score2 = score + (new_s - old_s) # swap leaves the KL term unchanged only if freq equal; recompute m2 = m.copy(); m2[i], m2[j] = m[j], m[i]; kl2 = kl_of(m2); score2 -= lam * T * (kl2 - kl) if score2 >= score or rng.random() < math.exp((score2 - score) / cur_T): m, kl, score = m2, kl2, score2 else: P[ciph.pos[i]] = m[i]; P[ciph.pos[j]] = m[j] cur_T = T0 * (1 - it / iters) + 1e-9 # polish: coordinate ascent, each unit to its best letter given all others, until no change for sweep in range(polish_sweeps): changed = 0 for i in np.argsort(-ciph.freq): st = starts[i]; cur_letter = m[i]; base = window_score(P, st); best_l, best_gain = cur_letter, 0.0 for new in letters: if new == cur_letter: continue P[ciph.pos[i]] = new; gain = window_score(P, st) - base m2 = m.copy(); m2[i] = new; gain -= lam * T * (kl_of(m2) - kl) if gain > best_gain: best_l, best_gain = new, gain P[ciph.pos[i]] = best_l if best_l != cur_letter: m[i] = best_l; kl = kl_of(m); score += best_gain; changed += 1 if verbose: print(f" polish sweep {sweep}: {changed} changes, score/char {score / T:.4f}", flush=True) if not changed: break if verbose: print(f" restart {r}: score/char {score / T:.4f} kl {kl:.3f}", flush=True) if score > best[0]: best = (score, m.copy(), score / T, kl) return best def decode(ciph, lm, m, units): out = [] for u in units: if u == " ": out.append(" ") elif u in ciph.idx: out.append(lm.alpha[m[ciph.idx[u]]]) else: out.append("?") return "".join(out) def verify(ciph, lm, m, units_heldout): txt = decode(ciph, lm, m, units_heldout) s = np.array([lm.idx[c] for c in txt if c in lm.idx], dtype=np.int64) sc = lm.score_seq(s) words = [w for w in txt.split(" ") if w and "?" not in w] hits = sum(w in lm.vocab for w in words) / max(len(words), 1) long = [w for w in words if len(w) >= 4] hits4 = sum(w in lm.vocab for w in long) / max(len(long), 1) return dict(score4=sc, dict_hits=hits, dict_hits_len4=hits4, n_words=len(words), sample=txt[:200]) def random_key_verify(ciph, lm, units_heldout, n=3, seed=0): rng = random.Random(seed); letters = [i for i in range(lm.L) if i != lm.sp] out = [] for k in range(n): m = np.array([rng.choice(letters) for _ in range(ciph.nU)]) out.append(verify(ciph, lm, m, units_heldout)) return dict(score4=sum(o["score4"] for o in out) / n, dict_hits=sum(o["dict_hits"] for o in out) / n, dict_hits_len4=sum(o["dict_hits_len4"] for o in out) / n) def plant_homophonic(words, n_units=80, max_homo=4, seed=0, keep_spaces=True): rng = random.Random(seed) cnt = collections.Counter(ch for w in words for ch in w) total = sum(cnt.values()); homos = {}; uid = 0 for a, c in cnt.most_common(): k = max(1, min(max_homo, round(c / total * n_units))) homos[a] = [f"u{uid + j}" for j in range(k)]; uid += k units = [] for w in words: for ch in w: units.append(rng.choice(homos[ch])) if keep_spaces: units.append(" ") return units, {u: a for a, us in homos.items() for u in us} def key_accuracy(ciph, lm, m, key): ok = tot = 0 for u, i in ciph.idx.items(): if u in key: tot += 1; ok += (lm.alpha[m[i]] == key[u]) return ok / max(tot, 1)