"""Exact supplementary checks; Python 3 standard library only. The accompanying written proofs establish the parameter-wide claims. Acceptance uses explicit exceptions and is unchanged under python -O. """ from fractions import Fraction as F from itertools import combinations from math import isqrt, sin, sqrt, pi def require(condition, label): if not condition: raise ValueError(label) # Sparse polynomials in two formal variables, with rational coefficients. def add(*polys): out = {} for poly in polys: for monomial, coefficient in poly.items(): out[monomial] = out.get(monomial, F(0)) + coefficient return {m: c for m, c in out.items() if c} def scale(poly, scalar): return {m: c * scalar for m, c in poly.items() if c * scalar} def mul(*polys): out = {(0, 0): F(1)} for poly in polys: nxt = {} for (i, j), c in out.items(): for (k, l), d in poly.items(): m = (i + k, j + l) nxt[m] = nxt.get(m, F(0)) + c * d out = {m: c for m, c in nxt.items() if c} return out class K: """Exact a + b*sqrt(5), with rational a and b.""" def __init__(self, a=0, b=0): if isinstance(a, K): self.a, self.b = a.a, a.b else: self.a, self.b = F(a), F(b) def __add__(self, other): other = K(other) return K(self.a + other.a, self.b + other.b) __radd__ = __add__ def __neg__(self): return K(-self.a, -self.b) def __sub__(self, other): return self + (-K(other)) def __rsub__(self, other): return K(other) + (-self) def __mul__(self, other): other = K(other) return K(self.a * other.a + 5 * self.b * other.b, self.a * other.b + self.b * other.a) __rmul__ = __mul__ def __truediv__(self, other): other = K(other) den = other.a**2 - 5 * other.b**2 if not den: raise ZeroDivisionError return self * K(other.a / den, -other.b / den) def __rtruediv__(self, other): return K(other) / self def __pow__(self, n): require(isinstance(n, int) and n >= 0, 'exponent') out = K(1) for _ in range(n): out *= self return out def __eq__(self, other): other = K(other) return self.a == other.a and self.b == other.b def sign(self): if not self.b: return (self.a > 0) - (self.a < 0) if not self.a: return (self.b > 0) - (self.b < 0) if self.a * self.b > 0: return 1 if self.a > 0 else -1 difference = self.a**2 - 5 * self.b**2 if difference > 0: return 1 if self.a > 0 else -1 return 1 if self.b > 0 else -1 def masses(alpha, beta): return ( 2 * (5 * beta - 1) / (15 * (beta - alpha) * (1 - alpha)), 2 * (1 - 5 * alpha) / (15 * (beta - alpha) * (1 - beta)), (3 - 5 * (alpha + beta) + 15 * alpha * beta) / (15 * (1 - alpha) * (1 - beta)), ) def field_text(x): x = K(x) return f'{x.a.numerator}/{x.a.denominator};{x.b.numerator}/{x.b.denominator}' def emit(name, x): print('VALUE\t' + name + '\t' + field_text(x)) def floorq(x): return x.numerator // x.denominator def decompose(x, h): require(h > 0, 'positive scale') q = floorq(x / h) return q, x - h*q def reverse_step(Q, S, nxt, rho): require(Q > 0 and S > 0 and 0 <= rho < S, 'ledger range') numerator = S*nxt + rho require(numerator % Q == 0, 'ledger divisibility') return numerator // Q def solve(matrix, rhs): """Unique exact solution to a full-column-rank rectangular system, or None.""" require(bool(matrix) and len(matrix) == len(rhs), 'system rows') cols = len(matrix[0]) require(cols > 0 and all(len(row) == cols for row in matrix), 'system columns') a = [[K(x) for x in row] + [K(y)] for row, y in zip(matrix, rhs)] row = 0 for col in range(cols): pivot = next((j for j in range(row, len(a)) if a[j][col] != 0), None) if pivot is None: return None a[row], a[pivot] = a[pivot], a[row] factor = a[row][col] a[row] = [x/factor for x in a[row]] for j in range(len(a)): if j != row: factor = a[j][col] a[j] = [x-factor*y for x, y in zip(a[j], a[row])] row += 1 if any(all(x == 0 for x in line[:cols]) and line[cols] != 0 for line in a): return None return [a[j][cols] for j in range(cols)] def moments(degree): require(isinstance(degree, int) and 0 <= degree <= 6, 'degree 0..6') return [K(0 if j % 2 else F(1, j+1)) for j in range(degree+1)] def positive_rule(nodes, degree): """Enumerate independent supports; intended for the six-node instrument.""" require(len(nodes) == 6 and all(nodes[i] != nodes[j] for i in range(6) for j in range(i)), 'six distinct nodes') target = moments(degree) for size in range(1, min(6, degree+1)+1): for support in combinations(range(6), size): a = [[nodes[i]**j for i in support] for j in range(degree+1)] weights = solve(a, target) if weights is not None and all(x.sign() >= 0 for x in weights): result = [K() for _ in nodes] for i, x in zip(support, weights): result[i] = x return result return None def recover_affine(nodes, records, errors): require(len(nodes) == len(records) == 6, 'six labeled records') require(isinstance(errors, int) and not isinstance(errors, bool) and errors >= 0, 'error budget') require(all(nodes[i] != nodes[j] for i in range(6) for j in range(i)), 'distinct nodes') good = [i for i, record in enumerate(records) if record is not None] require(2*errors + 6-len(good) <= 4, 'recovery budget') dim = len(records[good[0]]) require(dim > 0 and all(len(records[i]) == dim for i in good), 'coordinate dimension') candidates = [] for i, j in combinations(good, 2): q = [(records[j][k]-records[i][k])/(nodes[j]-nodes[i]) for k in range(dim)] m = [records[i][k]-nodes[i]*q[k] for k in range(dim)] predicted = [[m[k]+nodes[t]*q[k] for k in range(dim)] for t in range(6)] if sum(predicted[t] != records[t] for t in good) <= errors and predicted not in candidates: candidates.append(predicted) require(len(candidates) == 1, 'unique affine recovery') return candidates[0] def atan_interval(denominator, n=24): partial = sum((F((-1)**k, (2*k+1)*denominator**(2*k+1)) for k in range(n)), F()) other = partial + F((-1)**n, (2*n+1)*denominator**(2*n+1)) return min(partial, other), max(partial, other) def rejected(operation): try: operation() except (ValueError, ZeroDivisionError): return True return False def additional_checks(labels): for c in [F(-3, 2), F(0), F(1)]: for r in [F(-2, 3), F(0), F(1, 3), F(2, 3)]: for x in [c-2, c, c+3]: z = x for n in range(19): require(z == c+r**n*(x-c), 'affine closed form') nxt = c+r*(z-c) if 0 < r < 1: require((nxt-z)*(c-z) >= 0, 'affine monotonicity') z = nxt labels.append('affine iteration and monotonicity') for h in [F(1, 3), F(1), F(7, 2)]: for x in [F(-10**50, 3), F(-3, 2), F(0), F(5, 4), F(2**160+1)]: q, rho = decompose(x, h) require(x == h*q+rho and 0 <= rho < h, 'floor decomposition') require(floorq(F(-3, 2)) == -2, 'floor versus truncation') labels.append('signed floor decomposition and large integers') for r in [F(1, 3), F(2, 3)]: for h in [F(1, 5), F(1)]: for x in [F(-17, 3), F(0), F(399, 7)]: z, history = x, [] for n in range(1, 24): nxt = h*floorq(r*z/h) history.append(r*z-nxt) z = nxt err = sum((r**(n-1-j)*v for j, v in enumerate(history)), F()) require(r**n*x == z+err, 'unrolled residual identity') require(0 <= err < h*(1-r**n)/(1-r) < h/(1-r), 'residual bound') labels.append('quantized recurrence and geometric error bounds') Q, S = 339949771344778, 10**15 for seed in [-32767, 0, 399, 32767, 2**160+1]: d, ledger = seed, [] for _ in range(23): nxt, rho = divmod(Q*d, S) ledger.append(rho) require(reverse_step(Q, S, nxt, rho) == d, 'integer reverse step') d = nxt for rho in reversed(ledger): d = reverse_step(Q, S, d, rho) require(d == seed, 'complete ledger reverse') require(rejected(lambda: reverse_step(Q, S, 0, S)), 'reject remainder bound') require(rejected(lambda: reverse_step(Q, S, 0, 1)), 'reject nondivisible reverse') labels.append('integer ledgers and invalid-witness rejection') nodes = [K(j) for j in [-3, -2, -1, 1, 2, 3]] for degree in range(6): coeff = [K(F((-1)**j*(j+1), j+2)) for j in range(degree+1)] matrix = [[t**j for j in range(degree+1)] for t in nodes] values = [sum((a*b for a, b in zip(row, coeff)), K()) for row in matrix] require(solve(matrix, values) == coeff, 'polynomial coefficient recovery') m, q = [K(2), K(-1), K(7)], [K(3), K(2), K(-2)] original = [[m[k]+t*q[k] for k in range(3)] for t in nodes] corrupt = [list(row) for row in original] corrupt[0][0] += 100 corrupt[5][2] -= 200 require(recover_affine(nodes, corrupt, 2) == original, 'two errors') erased = [original[i] if i in [1, 4] else None for i in range(6)] require(recover_affine(nodes, erased, 0) == original, 'four erasures') mixed = [list(row) if i not in [1, 4] else None for i, row in enumerate(original)] mixed[0][1] += 123 require(recover_affine(nodes, mixed, 1) == original, 'mixed errors and erasures') malformed = [list(row) for row in original] malformed[0] = malformed[0][:1] require(rejected(lambda: recover_affine(nodes, malformed, 0)), 'reject truncated coordinate') require(rejected(lambda: recover_affine(nodes, original, 3)), 'reject excess budget') labels.append('polynomial interpolation and affine record recovery') phi, root = K(F(1, 2), F(1, 2)), K(0, 1) star = (root-1)/4 for s in [K(F(1, 4)), star, K(F(1, 2))]: A, B = s+phi**2, s+1 u, v = s/A, B/A x = [K(-1), -v, -u, u, v, K(1)] g = masses(u**2, v**2) w = [g[2]/2, g[1]/2, g[0]/2, g[0]/2, g[1]/2, g[2]/2] matrix = [[t**j for t in x] for j in range(6)] require(solve(matrix, moments(5)) == w, 'independent moment solve') for degree in [3, 4, 5, 6]: result = positive_rule(x, degree) expected = degree == 3 or (degree <= 5 and (s-star).sign() >= 0) require((result is not None) == expected, 'positive support feasibility') if result is not None: require(all(sum((result[i]*x[i]**j for i in range(6)), K()) == moments(degree)[j] for j in range(degree+1)), 'returned positive-rule witness') for j in range(6): target = K(0) if j % 2 else A**j/(j+1) require(sum((w[i]*(A*x[i])**j for i in range(6)), K()) == target, 'interval transport') outer = [-A, -B, B, A] ratio = lambda y: (y[2]-y[0])*(y[3]-y[1])/((y[2]-y[1])*(y[3]-y[0])) image = [(2*t+1)/(t+10) for t in outer] R = ratio(outer) require(ratio(image) == R, 'projective invariance') positive_root = (A+B)/(A-B) require(positive_root.sign() > 0 and positive_root**2 == R/(R-1), 'positive root branch') require(phi/2*(positive_root-root) == s, 'seed inverse') labels.extend(['independent quadrature solve and positive support search', 'interval rescaling and projective parameter recovery']) aL, aU = atan_interval(5) bL, bU = atan_interval(239) piL, piU = 16*aL-4*bU, 16*aU-4*bL xL, xU = piL/1296000, piU/1296000 require(0 < xL < xU < 1, 'sine-bound domain') sL, sU = xL-xL**3/6, xU require(0 < sL < sU < F(1, 4), 'one-arcsecond seed enclosure') require((K(sU)-star).sign() < 0, 'canonical seed below threshold') scale = 10**24 integer_root = isqrt(5*scale*scale) require(integer_root**2 <= 5*scale*scale < (integer_root+1)**2, 'integer square-root witness') phiL, phiU = (1+F(integer_root, scale))/2, (1+F(integer_root+1, scale))/2 lambdaL, lambdaU = phiL/(phiL+piU), phiU/(phiU+piL) require(0 < lambdaL < lambdaU < 1 and lambdaU-lambdaL < F(1, 10**24), 'contraction interval') emit('canonical_seed_lower', sL) emit('canonical_seed_upper', sU) emit('contraction_lower', lambdaL) emit('contraction_upper', lambdaU) labels.append('certified transcendental seed and contraction enclosures') require(rejected(lambda: decompose(F(1), F(0))), 'reject zero grid') require(rejected(lambda: positive_rule([K(0)]*6, 3)), 'reject duplicate nodes') require(rejected(lambda: moments(-1)), 'reject invalid degree') require(rejected(lambda: K(1)/K(0)), 'reject zero field divisor') labels.append('input-domain rejection') p = (1+sqrt(5))/2 s = sin(pi/1296000) a, b = (s/(s+p*p))**2, ((s+1)/(s+p*p))**2 amplification = 1-4*(5*b-1)/(15*(b-a)*(1-a)) require(abs(amplification-1.4944217144773972) < 1e-12, 'numerical illustration only') print(f'APPROX\tcanonical_amplification\t{amplification:.8f}') def run(): labels = [] one = {(0, 0): F(1)} a, b = {(1, 0): F(1)}, {(0, 1): F(1)} da, db = add(one, scale(a, -1)), add(one, scale(b, -1)) gap = add(b, scale(a, -1)) den = scale(mul(gap, da, db), 15) na = scale(mul(add(scale(b, 5), scale(one, -1)), db), 2) nb = scale(mul(add(one, scale(a, -5)), da), 2) nc = mul(add(scale(one, 3), scale(add(a, b), -5), scale(mul(a, b), 15)), gap) for j in range(3): lhs = add(mul(na, {(j, 0): F(1)}), mul(nb, {(0, j): F(1)}), nc) require(lhs == scale(den, F(1, 2*j + 1)), 'formal moment identity') labels.append('formal degree-five moment identities') sixth = add(scale(one, F(1, 7)), scale(add(one, a, b), -F(1, 5)), scale(add(a, b, mul(a, b)), F(1, 3)), scale(mul(a, b), -1)) require(sixth == scale(add(scale(one, 3), scale(add(a, b), -7), scale(mul(a, b), 35)), -F(2, 105)), 'formal sixth-degree obstruction') labels.append('formal sixth-degree obstruction') fourth = add(scale(one, F(1, 5)), scale(add(one, b), -F(1, 3)), b) require(fourth == scale(add(scale(b, 5), scale(one, -1)), F(2, 15)), 'formal fourth-degree obstruction') labels.append('formal fourth-degree obstruction') phi, root = K(F(1, 2), F(1, 2)), K(0, 1) star = (root - 1) / 4 rstar = (5 + 3 * root) / 10 require(phi**2 == phi + 1, 'golden identity') require((1 + phi**2)**2 == 5 * phi**2, 'endpoint bound') require((3 - 7 / phi**2).sign() > 0, 'uniform obstruction bound') require((F(2, 3) - 1 / phi**2).sign() > 0, 'positive cubic rule bound') labels.append('exact field identities and endpoint bounds') cases = [K(F(1, 1000000)), K(F(1, 10)), K(F(1, 4)), star, K(F(1, 2)), K(F(3, 4)), K(F(999999, 1000000))] for case_number, s in enumerate(cases): A, B = s + phi**2, s + 1 u, v = s / A, B / A alpha, beta = u**2, v**2 g = masses(alpha, beta) require(g[0].sign() == (s - star).sign(), 'inner mass sign') require(g[1].sign() > 0 and g[2].sign() > 0, 'other mass signs') for j in range(3): require(g[0]*alpha**j + g[1]*beta**j + g[2] == F(1, 2*j+1), 'exact field moment') sixth_value = -F(2, 105)*(3 - 7*(alpha+beta) + 35*alpha*beta) require(sixth_value.sign() < 0, 'sixth moment obstruction') t = 2 / (3 * (2 - alpha - beta)) require(t.sign() > 0 and (1 - 2*t).sign() > 0, 'cubic positivity') require(t*alpha + t*beta + 1 - 2*t == F(1, 3), 'cubic moment') ratio = (1+v)**2 / (4*v) require((rstar - ratio).sign() == (s-star).sign(), 'cross-ratio threshold') total_variation = sum((x * x.sign() for x in g), K()) require(total_variation == 1 + (0 if g[0].sign() >= 0 else -2*g[0]), 'amplification formula and attaining error signs') for name, value in [('s', s), ('inner', g[0]), ('middle', g[1]), ('outer', g[2]), ('cross_ratio', ratio), ('amplification', total_variation)]: emit(f'case_{case_number}_{name}', value) if s == star: require(beta == F(1, 5), 'boundary location') require(g == (K(0), K(F(5, 6)), K(F(1, 6))), 'boundary masses') require(ratio == rstar, 'boundary cross-ratio') labels.extend(['field moments and mass signs on seven seeds', 'strictly positive cubic rules on seven seeds', 'degree-six failure on seven seeds', 'cross-ratio threshold on seven seeds', 'sharp amplification on seven seeds', 'exact boundary weights and cross-ratio']) limiting = masses(K(0), 1 / phi**4) require(limiting[0] == -2 / (5 * phi), 'limiting negative mass') require(1 - 2*limiting[0] == 1 + 4 / (5 * phi), 'limiting amplification') labels.append('one-sided small-seed amplification limit') c, h, r, x = F(1, 2), F(1), F(1, 2), F(3, 4) floor = lambda z: z.numerator // z.denominator T = c + r*(x-c) require(h*floor(T/h)-c == -F(1, 2), 'unaligned grid counterexample') require(h*floor(r*(x-c)/h) == 0, 'recentered comparison') require(c + h*floor((T-c)/h)-c == h*floor(r*(x-c)/h), 'aligned grid') labels.append('grid-origin compatibility counterexample') additional_checks(labels) for label in labels: print('PASS\t' + label) print(f'TOTAL\t{len(labels)}\tPASS') if __name__ == '__main__': run()