"""Board stabrank: stabilizer rank of magic states (catalog W1-3, Q49; families of SPEC_v2_final s6 with review S1-S5). Instance families (n qubits, m = copies; bit q of x is qubit q, w = |x|) and their unnormalized target amplitudes: T n = m zeta^w, zeta = e^{i pi/4} basis over Q(i): 1, zeta CCZ n = 3m prod_b (-1)^{x_3b x_3b+1 x_3b+2} basis 1 face n = m (1 + sqrt3)^{m-w} (1 + i)^w basis 1, sqrt3 sqrtT n = m eta^w, eta = e^{i pi/8} basis 1, eta, eta^2, eta^3 CS n = 2m i^{sum_j x_2j x_2j+1} basis 1 cat n = m i^{w/2} for even w, 0 for odd w (one m-qubit state) basis 1 C3Z n = 4m (-1)^{sum_j x_4j x_4j+1 x_4j+2 x_4j+3} basis 1 A missing or unknown family, or a wrong copies/n relation, is an invalid instance (status error), never a fall-through. Artifact {"n": n, "states": [{"h": bits(n), "G": [bits(n)] * k, "Q": [[i, j], ...], "L": [l_0 .. l_{k-1}]}, ...]}. State vector: psi(h XOR yG) = i^{sum_i l_i y_i} * (-1)^{sum_{(i,j) in Q} y_i y_j} for y in {0,1}^k, 0 elsewhere. G must have full row rank (then psi is a stabilizer state; every stabilizer state has this form up to a scalar). Valid iff the target lies in the complex span of the states. Score r = number of states submitted (min). Exact decision (no floating point, no probabilistic acceptance): 1. M = amplitude matrix (2^n x r, entries 0, ±1, ±i). Target t = sum_c beta_c t_c with t_c in Z[i]^{2^n} and {beta_c} the family's basis above, which is linearly independent over Q(i). M has entries in Q(i), so t in span_C(M) iff every component t_c is in span_{Q(i)}(M). 2. Modulo a fixed prime p = 1 (mod 8) (i -> a square root of -1), Gauss-Jordan elimination picks pivot columns C and rows R with M[R, C] invertible mod p, hence invertible over Q(i) (a minor that is nonzero mod p is nonzero). 3. Exact rational solves (python-flint) on the rows R. Every non-pivot column is checked to lie in span(M_C) EXACTLY on all rows; if one does not, p was a bad prime and the next fixed prime is used. Then span(M) = span(M_C) and M_C has full column rank, so the unique candidate coefficients from rows R decide: each component t_c is in the span iff its candidates reproduce it EXACTLY on all rows. 4. "Exactly on all rows": both sides are integer vectors after clearing denominators, bounded by a computed B (from the actual component amplitudes); they are compared modulo distinct primes whose product exceeds 2B + 1, which proves equality over Z. If the fixed primes cannot cover the bound, or all three pivot primes are bad, the status is error. """ import math import numpy as np import flint from common import Reject, exact_keys, want_int, want_list, is_int, bitstring, load_instance, rank, short MAX_R = 256 WIDTH = {"T": 1, "CCZ": 3, "face": 1, "sqrtT": 1, "CS": 2, "cat": 1, "C3Z": 4} # qubits per copy class Exhausted(Exception): """The fixed CRT primes cannot cover the exactness bound (reported as status error, never a verdict).""" def _is_prime(n): if n < 2: return False for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37): if n % p == 0: return n == p d, s = n - 1, 0 while d % 2 == 0: d //= 2 s += 1 for a in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37): # deterministic for n < 3.3e24 x = pow(a, d, n) if x in (1, n - 1): continue for _ in range(s - 1): x = x * x % n if x == n - 1: break else: return False return True def _primes_1mod8(count, start=2 ** 31): out, p = [], start - (start % 8) + 1 while len(out) < count: p -= 8 if _is_prime(p): out.append(p) return out PRIMES = _primes_1mod8(200) # fixed and deterministic: 200 primes p = 1 (mod 8) just below 2^31 def _root8(p): """An element of multiplicative order 8 mod p.""" for g in range(2, 200): w = pow(g, (p - 1) // 8, p) if pow(w, 4, p) == p - 1: return w raise RuntimeError("no 8th root") def parse_state(s, n, idx): if not isinstance(s, dict): raise Reject(f"states[{idx}] must be an object") exact_keys(s, ["h", "G", "Q", "L"], where=f"states[{idx}]") h = bitstring(s["h"], n, f"states[{idx}].h") G = [bitstring(g, n, f"states[{idx}].G[{t}]") for t, g in enumerate(want_list(s["G"], f"states[{idx}].G", max_len=n))] k = len(G) if rank(G) != k: raise Reject(f"states[{idx}].G does not have full row rank over GF(2)") L = want_list(s["L"], f"states[{idx}].L", max_len=k, min_len=k) for t, l in enumerate(L): want_int(l, 0, 3, f"states[{idx}].L[{t}]") seen = set() Q = [] for t, pr in enumerate(want_list(s["Q"], f"states[{idx}].Q", max_len=k * (k + 1) // 2)): if not (isinstance(pr, list) and len(pr) == 2 and is_int(pr[0]) and is_int(pr[1]) and 0 <= pr[0] <= pr[1] < k): raise Reject(f"states[{idx}].Q[{t}] must be a pair [i, j] with 0 <= i <= j < k = {k}") if tuple(pr) in seen: raise Reject(f"states[{idx}].Q repeats the pair {short(pr)}") seen.add(tuple(pr)) Q.append(tuple(pr)) return h, G, Q, L def state_vector(h, G, Q, L, n): """Return (support indices, exponents e in Z_4) with amplitude i^e at each support index.""" k = len(G) y = np.arange(1 << k, dtype=np.int64) Y = [(y >> t) & 1 for t in range(k)] x = np.full(1 << k, h, dtype=np.int64) e = np.zeros(1 << k, dtype=np.int64) for t in range(k): x ^= Y[t] * G[t] e += L[t] * Y[t] for (a, b) in Q: e += 2 * (Y[a] & Y[b]) return x, e % 4 def _ipow(e, mask=1): """(re, im) int64 arrays of mask * i^e.""" e = e % 4 return (np.where(e == 0, 1, np.where(e == 2, -1, 0)) * mask, np.where(e == 1, 1, np.where(e == 3, -1, 0)) * mask) def target(fam, m, n): """Components [(re, im), ...] of the target over Q(i), one per basis element of the family (module docstring).""" x = np.arange(1 << n, dtype=np.int64) bit = [(x >> q) & 1 for q in range(n)] w = sum(bit) if fam == "T": # zeta^w = i^(w // 2) zeta^(w % 2) return [_ipow(w // 2, w % 2 == c) for c in range(2)] if fam == "sqrtT": # eta^w = i^(w // 4) eta^(w % 4) return [_ipow(w // 4, w % 4 == c) for c in range(4)] if fam == "cat": return [_ipow(w // 2, w % 2 == 0)] if fam == "CS": return [_ipow(sum(bit[2 * j] & bit[2 * j + 1] for j in range(m)))] if fam in ("CCZ", "C3Z"): d = WIDTH[fam] return [_ipow(2 * sum(np.bitwise_and.reduce(bit[d * j:d * j + d]) for j in range(m)))] if fam == "face": # (1 + sqrt3)^j = a_j + b_j sqrt3, (1 + i)^w = g_w; components a_{m-w} g_w and b_{m-w} g_w a, b, g = [1], [0], [(1, 0)] for _ in range(m): a, b, g = a + [a[-1] + 3 * b[-1]], b + [a[-1] + b[-1]], g + [(g[-1][0] - g[-1][1], g[-1][0] + g[-1][1])] A, B = np.array(a, dtype=np.int64)[m - w], np.array(b, dtype=np.int64)[m - w] gr, gi = np.array([v[0] for v in g], dtype=np.int64)[w], np.array([v[1] for v in g], dtype=np.int64)[w] return [(A * gr, A * gi), (B * gr, B * gi)] raise ValueError(f"unknown family {fam!r}") def instance(I): """(family, copies, n, r_max) of a parsed instance; ValueError if it is not a valid stabrank instance (review S1, S2).""" if not isinstance(I, dict): raise ValueError("not a JSON object") fam, m, n, r_max = I.get("family"), I.get("copies"), I.get("n"), I.get("r_max") if not (isinstance(fam, str) and fam in WIDTH): raise ValueError(f"unknown family {short(fam)}") if not (is_int(m) and m >= 1 and is_int(n) and n == WIDTH[fam] * m): raise ValueError(f"family {fam} needs n = {WIDTH[fam]} * copies") if not (is_int(r_max) and 1 <= r_max <= MAX_R): raise ValueError(f"r_max must be an integer in [1, {MAX_R}]") return fam, m, n, r_max def _modp_pivots(Mr, Mi, p): """Gauss-Jordan mod p on M = Mr + i Mi. Returns (rows R, cols C) with M[R, C] invertible mod p.""" im = pow(_root8(p), 2, p) B = (Mr % p + (Mi % p) * im) % p N, r = B.shape used = np.zeros(N, dtype=bool) R, C = [], [] for j in range(r): col = B[:, j].copy() col[used] = 0 nz = np.flatnonzero(col) if len(nz) == 0: continue i = nz[0] inv = pow(int(B[i, j]), p - 2, p) B[i] = (B[i] * inv) % p f = B[:, j].copy() f[i] = 0 rows = np.flatnonzero(f) if len(rows): B[rows] = (B[rows] - (f[rows, None] * B[i][None, :]) % p) % p used[i] = True R.append(int(i)) C.append(j) return R, C def _solve_exact(Ar, Ai, rhs): """Solve (Ar + i Ai) c = b over Q(i) for each complex RHS column b = (br, bi) in rhs. Returns list of (D, a, b): c = (a + i b) / D with integer vectors a, b and integer D > 0.""" rho = Ar.shape[0] big = [[0] * (2 * rho) for _ in range(2 * rho)] for u in range(rho): for v in range(rho): big[u][v] = int(Ar[u, v]) big[u][rho + v] = -int(Ai[u, v]) big[rho + u][v] = int(Ai[u, v]) big[rho + u][rho + v] = int(Ar[u, v]) A = flint.fmpq_mat(2 * rho, 2 * rho, [x for row in big for x in row]) Bm = flint.fmpq_mat(2 * rho, len(rhs), [int(rhs[c][0][u]) if u < rho else int(rhs[c][1][u - rho]) for u in range(2 * rho) for c in range(len(rhs))]) X = A.solve(Bm) out = [] for c in range(len(rhs)): vals = [X[u, c] for u in range(2 * rho)] D = 1 for q in vals: D = D * int(q.q) // math.gcd(D, int(q.q)) ints = [int(q.p) * (D // int(q.q)) for q in vals] out.append((D, ints[:rho], ints[rho:])) return out def _check_exact(MCr, MCi, D, a, b, br, bi): """Exactly test (MCr + i MCi)(a + i b) == D * (br + i bi) on all rows, via enough primes for a CRT proof.""" bound = (sum(abs(x) for x in a) + sum(abs(x) for x in b)) + D * (int(np.abs(br).max()) + int(np.abs(bi).max())) need, prod, k = 2 * bound + 1, 1, 0 while prod <= need: if k == len(PRIMES): raise Exhausted p = PRIMES[k] k += 1 prod *= p am = np.array([x % p for x in a], dtype=np.int64) bm = np.array([x % p for x in b], dtype=np.int64) lhs_r = (MCr @ am - MCi @ bm) % p lhs_i = (MCr @ bm + MCi @ am) % p Dm = D % p if np.any((lhs_r - (Dm * (br % p)) % p) % p) or np.any((lhs_i - (Dm * (bi % p)) % p) % p): return False return True def verify(inst, art): return check(load_instance("stabrank", inst), art) def check(I, art): """Verify an artifact against a parsed instance (verify() reads the instance file).""" try: fam, m, n, r_max = instance(I) except ValueError as e: return {"status": "error", "detail": {"reason": f"invalid instance file: {e}"}} exact_keys(art, ["n", "states"]) want_int(art["n"], n, n, "n") states = want_list(art["states"], "states", max_len=r_max, min_len=1) r = len(states) N = 1 << n Mr = np.zeros((N, r), dtype=np.int64) Mi = np.zeros((N, r), dtype=np.int64) for j, s in enumerate(states): h, G, Q, L = parse_state(s, n, j) x, e = state_vector(h, G, Q, L, n) Mr[x, j] = np.where(e == 0, 1, np.where(e == 2, -1, 0)) Mi[x, j] = np.where(e == 1, 1, np.where(e == 3, -1, 0)) comps = target(fam, m, n) try: for p in PRIMES[:3]: # pivots from the first rows that have any support, widening the block until it reaches full column rank R, C = None, None for block in (min(N, 4 * r + 64), N): rows = np.arange(block) Rb, Cb = _modp_pivots(Mr[rows], Mi[rows], p) if len(Cb) == r or block == N: R, C = Rb, Cb break rho = len(R) nonpiv = [j for j in range(r) if j not in C] MCr, MCi = Mr[:, C], Mi[:, C] Ar, Ai = MCr[R], MCi[R] rhs = [(Mr[R, j], Mi[R, j]) for j in nonpiv] + [(tr[R], ti[R]) for tr, ti in comps] sols = _solve_exact(Ar, Ai, rhs) if not all(_check_exact(MCr, MCi, *sol, Mr[:, j], Mi[:, j]) for sol, j in zip(sols, nonpiv)): continue # bad prime: a modular rank loss, not a verdict if not all(_check_exact(MCr, MCi, *sol, tr, ti) for sol, (tr, ti) in zip(sols[len(nonpiv):], comps)): raise Reject(f"the target is not in the span of the {r} states (exact rank of the states: {rho})") score = {"r": r} if fam == "T": score["log2r_over_m"] = round(math.log2(r) / n, 6) return {"status": "accepted", "score": score, "detail": {"n": n, "r": r, "exact_rank_of_states": rho, "note": ("target verified in the span exactly on all 2^n amplitudes" + ("" if rho == r else f"; {r - rho} states are redundant, an independent subset " f"of {rho} states also spans the target"))}} except Exhausted: return {"status": "error", "detail": {"reason": f"the exactness proof needs more than the {len(PRIMES)} fixed " "CRT primes; please report"}} return {"status": "error", "detail": {"reason": "three fixed primes were all bad for this matrix; please report"}}