#!/usr/bin/env python3 """Exact algebra and decoding experiments for the six OMEGA landmarks. The proofs are in OMEGA_Object_Exploration_02.md. Checks use an indeterminate seed s over Q(sqrt(5)), so the exact algebra is not a floating approximation to the canonical one-arcsecond construction. Decimal values are illustrative. """ import argparse import hashlib import itertools import json from pathlib import Path import sympy as sy from mpmath import mp seed = sy.Symbol("s") field = sy.QQ.algebraic_field(sy.sqrt(5)).frac_field(seed) F = field.from_sympy ZERO, ONE = field.zero, field.one phi = F((1 + sy.sqrt(5)) / 2) s = F(seed) A, B = s + phi**2, s + ONE nodes = (-A, -B, -s, s, B, A) checks = [] def equal(x, y): # Algebraic coefficient units need not have identical stored denominators. return x - y == ZERO def vequal(x, y): return len(x) == len(y) and all( equal(a, b) for a, b in zip(x, y) ) def record_equal(x, y): return len(x) == len(y) and all( vequal(a, b) for a, b in zip(x, y) ) def require(condition, name): if not condition: raise ArithmeticError(name) checks.append(name) print("PASS", name, flush=True) def vadd(x, y): if len(x) != len(y): raise ValueError("Coordinate dimensions differ") return tuple(a + b for a, b in zip(x, y)) def vsub(x, y): if len(x) != len(y): raise ValueError("Coordinate dimensions differ") return tuple(a - b for a, b in zip(x, y)) def smul(a, x): return tuple(a * b for b in x) def affine_record(m, q): return tuple(vadd(m, smul(t, q)) for t in nodes) def _require_six_slots(received): if len(received) != 6: raise ValueError("Decoders require exactly six labeled input slots") def _require_max_errors(max_errors): if type(max_errors) is not int or max_errors < 0: raise ValueError("max_errors must be a nonnegative built-in integer") def _erasure_count(received): return sum(value is None for value in received) def _require_budget(max_errors, received, bound): _require_six_slots(received) _require_max_errors(max_errors) r = _erasure_count(received) if 2 * max_errors + r > bound: raise ValueError("Requested recovery budget is outside the proved bound") def _affine_vector(value): if isinstance(value, (str, bytes)): raise ValueError("Scalar positions are not affine coordinate records") try: length = len(value) except TypeError as exc: raise ValueError("Scalar positions are not affine coordinate records") from exc if length == 0: raise ValueError("Empty vectors are not affine coordinate records") return tuple(value) def decode_affine(received, max_errors=2): _require_budget(max_errors, received, 4) retained = [i for i, value in enumerate(received) if value is not None] vectors = [_affine_vector(received[i]) for i in retained] if vectors: dim = len(vectors[0]) if any(len(item) != dim for item in vectors): raise ValueError("Inconsistent coordinate dimensions") accepted = [] for i, j in itertools.combinations(retained, 2): q = smul(ONE / (nodes[j] - nodes[i]), vsub(received[j], received[i])) if all(value == ZERO for value in q): continue m = vsub(received[i], smul(nodes[i], q)) predicted = affine_record(m, q) matches = sum(vequal(predicted[k], received[k]) for k in retained) if matches >= len(retained) - max_errors: if not any(record_equal(predicted, item[1]) for item in accepted): accepted.append(((m, q), predicted)) if len(accepted) != 1: raise ValueError(f"No unique admissible baseline: {len(accepted)} candidates") return accepted[0] def mobius_value(h, t): a, b, c, d = h denominator = c * t + d if denominator == ZERO: raise ValueError("Point at infinity") return (a * t + b) / denominator def fit_mobius(z, y): z0, z1, z2 = z y0, y1, y2 = y if any(equal(a, b) for a, b in itertools.combinations(y, 2)): raise ValueError("Repeated image coordinates") az, bz = z1 - z2, z1 - z0 ay, by = y1 - y2, y1 - y0 h = (by * y2 * az - ay * y0 * bz, -by * y2 * az * z0 + ay * y0 * bz * z2, by * az - ay * bz, -by * az * z0 + ay * bz * z2) if h[0] * h[3] - h[1] * h[2] == ZERO: raise ValueError("Degenerate projectivity") return h def decode_projective(received, max_errors=1): _require_budget(max_errors, received, 3) retained = [i for i, value in enumerate(received) if value is not None] candidates = [] for triple in itertools.combinations(retained, 3): try: h = fit_mobius(tuple(nodes[i] for i in triple), tuple(received[i] for i in triple)) predicted = tuple(mobius_value(h, t) for t in nodes) except ValueError: continue matches = sum(equal(predicted[i], received[i]) for i in retained) if matches >= len(retained) - max_errors: if not any(vequal(predicted, item) for item in candidates): candidates.append(predicted) if len(candidates) != 1: raise ValueError(f"No unique admissible projectivity: {len(candidates)} candidates") return candidates[0] def cross_ratio(x): x0, x1, x2, x3 = x return (x2 - x0) * (x3 - x1) / ((x2 - x1) * (x3 - x0)) def _raises(fn, *args, **kwargs): try: fn(*args, **kwargs) except ValueError: return True return False def main(): parser = argparse.ArgumentParser() parser.add_argument("--output", type=Path, required=True) args = parser.parse_args() ratio = cross_ratio((-A, -B, B, A)) require(equal(ratio, ONE + phi**2 / (4 * A * B)), "Outer cross-ratio formula") require(equal(ratio / (ratio - ONE), (2 * s / phi + F(sy.sqrt(5)))**2), "Seed inversion identity before positive square root") ratio_expr = field.to_sympy(ratio) expected_derivative = -phi**2 * (2 * s + phi**2 + ONE) / (4 * A**2 * B**2) require(equal(F(sy.diff(ratio_expr, seed)), expected_derivative), "Strictly negative cross-ratio derivative for s>0") h_true = tuple(F(sy.Integer(v)) for v in (2, 3, 1, 5)) outer_projected = tuple(mobius_value(h_true, t) for t in (-A, -B, B, A)) require(equal(cross_ratio(outer_projected), ratio), "Cross-ratio under an exact non-affine projectivity") # Coordinate values lie in the symbolic field Q(sqrt(5))(s). m = tuple(F(sy.Integer(v)) for v in (2, 3, 5)) q = tuple(F(sy.Integer(v)) for v in (3, -1, 2)) truth = affine_record(m, q) require(len(set(nodes)) == 6, "Six distinct symbolic evaluation nodes") recovered_pairs = 0 for bad in itertools.combinations(range(6), 2): received = list(truth) for j, index in enumerate(bad): offset = tuple(F(sy.Integer(v)) for v in (7 + j, -4 - j, 9 + 2*j)) received[index] = vadd(received[index], offset) (decoded_m, decoded_q), predicted = decode_affine(received) if not vequal(decoded_m, m) or not vequal(decoded_q, q) or not record_equal(predicted, truth): raise ArithmeticError(f"Two-error recovery failed at {bad}") recovered_pairs += 1 require(recovered_pairs == 15, "All 15 two-error location patterns decoded") recovered_erasures = 0 for good in itertools.combinations(range(6), 2): received = [truth[i] if i in good else None for i in range(6)] (decoded_m, decoded_q), predicted = decode_affine(received, max_errors=0) if not vequal(decoded_m, m) or not vequal(decoded_q, q) or not record_equal(predicted, truth): raise ArithmeticError(f"Four-erasure recovery failed at {good}") recovered_erasures += 1 require(recovered_erasures == 15, "All 15 four-erasure patterns reconstructed") # Sharp limitation: an observation at distance three from two valid records. other = affine_record(vadd(m, (ONE, ZERO, ZERO)), q) ambiguous = truth[:3] + other[3:] require(sum(not vequal(x, y) for x, y in zip(ambiguous, truth)) == 3 and sum(not vequal(x, y) for x, y in zip(ambiguous, other)) == 3, "Three-error ambiguity between two valid affine records") projected = tuple(mobius_value(h_true, t) for t in nodes) for bad in range(6): received = list(projected) received[bad] += ONE if not vequal(decode_projective(received), projected): raise ArithmeticError(f"Projective one-error recovery failed at {bad}") require(True, "All six one-error projective location patterns decoded") for good in itertools.combinations(range(6), 3): received = [projected[i] if i in good else None for i in range(6)] if not vequal(decode_projective(received, max_errors=0), projected): raise ArithmeticError(f"Projective three-erasure recovery failed at {good}") require(True, "All 20 three-erasure projective patterns reconstructed") r0, r1 = nodes[:2] kappa = F(sy.Rational(1, 10)) h_second = (ONE - kappa * (r0 + r1), kappa * r0 * r1, -kappa, ONE) second = tuple(mobius_value(h_second, t) for t in nodes) require(sum(not equal(x, y) for x, y in zip(nodes, second)) == 4, "Two distinct projective records can differ at exactly four nodes") projective_ambiguous = nodes[:4] + second[4:] require(sum(not equal(x, y) for x, y in zip(projective_ambiguous, nodes)) == 2 and sum(not equal(x, y) for x, y in zip(projective_ambiguous, second)) == 2, "Two-error ambiguity for projective coordinates") # Uniform probability integral on [-A,A]. def mean_monomial(k): return ZERO if k % 2 else A**k / (k + 1) D = 2 * A**2 - s**2 - B**2 w = 2 * A**2 / (3 * D) w_outer = ONE - 2 * w weights = (w_outer / 2, w / 2, w / 2, w / 2, w / 2, w_outer / 2) for degree in range(4): if not equal(sum((weight * t**degree for weight, t in zip(weights, nodes)), ZERO), mean_monomial(degree)): raise ArithmeticError(f"Positive rule moment failed at degree {degree}") require(True, "Six-node symmetric rule integrates degrees zero through three") obstruction_mean = mean_monomial(4) - (A**2 + B**2)*mean_monomial(2) + A**2*B**2 require(equal(obstruction_mean, 2*A**2*(5*B**2 - A**2)/15), "Degree-four obstruction integral") obstruction_values = tuple((t*t-B*B)*(t*t-A*A) for t in nodes) require(obstruction_values[0] == obstruction_values[1] == ZERO and obstruction_values[4] == obstruction_values[5] == ZERO and equal(obstruction_values[2], (s*s-B*B)*(s*s-A*A)) and equal(obstruction_values[3], obstruction_values[2]), "Obstruction vanishes at four nodes and is positive at the inner pair") threshold = (F(sy.sqrt(5)) - ONE) / 4 require(equal(phi**2 - F(sy.sqrt(5)), (F(sy.sqrt(5)) - ONE)*threshold), "Exact threshold s=(sqrt(5)-1)/4") require(equal(3*threshold - 4*threshold**3, ONE - 2*threshold**2), "Threshold sine satisfies sin(3 delta)=cos(2 delta), giving delta=18 degrees") lobatto_nodes = (-ONE, -ONE/F(sy.sqrt(5)), ONE/F(sy.sqrt(5)), ONE) lobatto_weights = tuple(F(sy.Rational(v, 12)) for v in (1, 5, 5, 1)) for degree in range(6): target = ZERO if degree % 2 else ONE/(degree + 1) actual = sum((weight*t**degree for weight,t in zip(lobatto_weights,lobatto_nodes)), ZERO) if not equal(actual, target): raise ArithmeticError(f"Boundary degree-five rule failed at {degree}") require(True, "At seed angle 36 degrees the four outer nodes give positive degree-five exactness") signed_inner_mass = obstruction_mean / ((s*s-B*B)*(s*s-A*A)) require(equal(signed_inner_mass, 2*A**2*(5*B**2-A**2)/(15*(s*s-B*B)*(s*s-A*A))), "Every symmetric degree-five rule has the displayed negative inner mass") truncated = list(truth) truncated[0] = truth[0][:2] empty = list(truth) empty[1] = () scalar = list(truth) scalar[2] = ONE mixed_dim = list(truth) mixed_dim[3] = truth[3] + (ONE,) require( (not vequal(truth[0][:2], truth[0])) and _raises(vadd, truth[0][:2], truth[0]) and _raises(vsub, truth[0], truth[0][:1]) and _raises(decode_affine, truncated) and _raises(decode_affine, empty) and _raises(decode_affine, scalar) and _raises(decode_affine, mixed_dim), "Malformed affine coordinate records are rejected", ) too_short = list(truth)[:5] too_long = list(truth) + [None] six_erasures = [None] * 6 three_none = [None, None, None, projected[3], projected[4], projected[5]] require( _raises(decode_affine, too_short) and _raises(decode_affine, too_long) and _raises(decode_affine, list(truth), max_errors=True) and _raises(decode_affine, list(truth), max_errors=-1) and _raises(decode_affine, list(truth), max_errors=1.5) and _raises(decode_affine, list(truth), max_errors=3) and _raises(decode_affine, six_erasures, max_errors=0) and _raises(decode_projective, list(projected), max_errors=2) and _raises(decode_projective, three_none, max_errors=1), "Invalid record counts and recovery budgets are rejected", ) mp.dps = 70 ph = (1 + mp.sqrt(5))/2 canonical_s = mp.sin(mp.pi / 1296000) outer, middle = canonical_s + ph**2, canonical_s + 1 canonical_ratio = (outer + middle)**2 / (4*outer*middle) canonical_w = 2*outer**2 / (3*(2*outer**2-canonical_s**2-middle**2)) numbers = { "s": mp.nstr(canonical_s, 45), "outer_cross_ratio": mp.nstr(canonical_ratio, 45), "seed_from_cross_ratio": mp.nstr(ph/2*(mp.sqrt(canonical_ratio/(canonical_ratio-1))-mp.sqrt(5)), 45), "epsilon_recovered_arcseconds": mp.nstr(2*mp.asin(ph/2*(mp.sqrt(canonical_ratio/(canonical_ratio-1))-mp.sqrt(5)))*648000/mp.pi, 40), "inner_pair_mass": mp.nstr(canonical_w, 35), "middle_pair_mass": mp.nstr(canonical_w, 35), "outer_pair_mass": mp.nstr(1-2*canonical_w, 35), "obstruction_uniform_mean": mp.nstr(2*outer**2*(5*middle**2-outer**2)/15, 35), } args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(json.dumps({ "symbolic_domain": "Q(sqrt(5))(s); s is an indeterminate", "proof_scope": "Exact statements under the hypotheses in the accompanying note", "passed_check_groups": len(checks), "checks": checks, "affine_error_location_patterns": 15, "affine_erasure_location_patterns": 15, "projective_error_location_patterns": 6, "projective_erasure_location_patterns": 20, "canonical_decimals_are_illustrations": numbers, "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "sympy_version": sy.__version__, "limitations": [ "Symbolic checks corroborate the written proofs.", "Exact sparse corruptions are not an ordinary measurement-noise model.", "Labels and the declared seed parameter must be trusted for decoding.", "No claim of prior-art priority or newly invented general methods.", ], }, indent=2) + "\n") print("COMPLETED", len(checks), "check groups", flush=True) if __name__ == "__main__": main()