#!/usr/bin/env python3 """Independent D1 audit. Derived from definitions; no kx/scorer/audit imports. Nonnegative polynomial coefficients are packed as radix-2**(8*w) digits. GMP multiplies the packed integers exactly. Choosing a radix strictly above the coefficient bound prevents carries, so unpacking is exact convolution. """ import argparse import ctypes as C import ctypes.util from decimal import Decimal, localcontext from fractions import Fraction import hashlib import json import math from pathlib import Path import random import time class MPZ(C.Structure): _fields_ = [("alloc", C.c_int), ("size", C.c_int), ("limbs", C.POINTER(C.c_ulong))] GMP = C.CDLL(ctypes.util.find_library("gmp")) INIT = GMP.__gmpz_init CLEAR = GMP.__gmpz_clear IMPORT = GMP.__gmpz_import EXPORT = GMP.__gmpz_export MUL = GMP.__gmpz_mul INIT.argtypes = CLEAR.argtypes = [C.POINTER(MPZ)] IMPORT.argtypes = [C.POINTER(MPZ), C.c_size_t, C.c_int, C.c_size_t, C.c_int, C.c_size_t, C.c_void_p] EXPORT.argtypes = [C.c_void_p, C.POINTER(C.c_size_t), C.c_int, C.c_size_t, C.c_int, C.c_size_t, C.POINTER(MPZ)] MUL.argtypes = [C.POINTER(MPZ), C.POINTER(MPZ), C.POINTER(MPZ)] def exact_convolution(a, b): assert a and b and min(a) >= 0 and min(b) >= 0 bound = min(sum(a) * max(b), sum(b) * max(a)) width = max(1, (bound.bit_length() + 7) // 8) assert bound < 1 << (8 * width) packed_a = b"".join(x.to_bytes(width, "little") for x in a) packed_b = b"".join(x.to_bytes(width, "little") for x in b) x, y, z = MPZ(), MPZ(), MPZ() for v in (x, y, z): INIT(C.byref(v)) try: IMPORT(C.byref(x), len(packed_a), -1, 1, 0, 0, packed_a) IMPORT(C.byref(y), len(packed_b), -1, 1, 0, 0, packed_b) MUL(C.byref(z), C.byref(x), C.byref(y)) count = C.c_size_t() out = C.create_string_buffer(len(packed_a) + len(packed_b) + 1) EXPORT(out, C.byref(count), -1, 1, 0, 0, C.byref(z)) length = len(a) + len(b) - 1 assert count.value <= length * width raw = out.raw[:length * width] coeffs = [int.from_bytes(raw[i:i+width], "little") for i in range(0, length * width, width)] finally: for v in (x, y, z): CLEAR(C.byref(v)) assert max(coeffs) <= bound assert sum(coeffs) == sum(a) * sum(b) assert sum(i*x for i, x in enumerate(coeffs)) == ( sum(i*x for i, x in enumerate(a))*sum(b) + sum(a)*sum(i*x for i, x in enumerate(b))) assert sum(i*i*x for i, x in enumerate(coeffs)) == ( sum(i*i*x for i, x in enumerate(a))*sum(b) + 2*sum(i*x for i, x in enumerate(a))*sum(i*x for i, x in enumerate(b)) + sum(a)*sum(i*i*x for i, x in enumerate(b))) return coeffs, {"coefficient_bound": str(bound), "radix_bits": 8*width, "coefficient_count": len(coeffs), "moment_checks": "pass"} def direct_coeff(a, b, k): return sum(a[i]*b[k-i] for i in range(max(0, k-len(b)+1), min(len(a), k+1))) def decimal(q): with localcontext() as ctx: ctx.prec = 75 return str(Decimal(q.numerator) / Decimal(q.denominator)) def rational(q): return {"fraction": str(q), "decimal": decimal(q)} def parse(path): raw = path.read_bytes() values = [Fraction(t) for line in raw.decode().splitlines() for t in line.split("#", 1)[0].split()] assert values and min(values) >= 0 and sum(values) > 0 d = math.lcm(*(v.denominator for v in values)) a = [v.numerator*(d//v.denominator) for v in values] return raw, values, d, a def audit(path, problem): started = time.monotonic() raw, values, d, a = parse(path) n, s = len(a), sum(a) result = {"file": str(path), "problem": problem, "sha256": hashlib.sha256(raw).hexdigest(), "n": n, "decimal_scale": str(d), "integer_sum": str(s), "raw_sum": rational(Fraction(s, d)), "raw_min": rational(min(values)), "raw_max": rational(max(values)), "zero_cells": a.count(0), "nonnegative": True, "positive_integral": True} if problem == "ac1": b = a coeffs, method = exact_convolution(a, b) peak = max(coeffs) score = Fraction(2*n*peak, s*s) result["raw_integral"] = rational(Fraction(s, 2*n*d)) result["normalization"] = "scale-invariant ratio; mass can be normalized exactly" result["max_knots"] = [str(Fraction(k+1-n, 2*n)) for k,c in enumerate(coeffs) if c == peak] else: assert n*max(a) <= 2*s, "normalized height exceeds 1; no repair performed" b = [2*s - n*x for x in reversed(a)] coeffs, method = exact_convolution(a, b) peak = max(coeffs) score = Fraction(peak, 2*s*s) result.update({"normalized_integral": "1", "normalized_min": rational(Fraction(n*min(a),2*s)), "normalized_max": rational(Fraction(n*max(a),2*s)), "max_slack": rational(1-Fraction(n*max(a),2*s)), "raw_unit_bounds": max(values) <= 1, "repair_required": False, "max_shifts_for_h_x_times_complement_x_plus_t": [str(Fraction(2*(n-1-k), n)) for k,c in enumerate(coeffs) if c == peak]}) # Entire independent O(n^2) exact path, no packing/multiplication library. direct = [direct_coeff(a, b, k) for k in range(2*n-1)] assert direct == coeffs result["all_shifts_direct_integer_check"] = "pass" rng = random.Random(20261009) top = sorted(range(len(coeffs)), key=coeffs.__getitem__, reverse=True)[:8] indices = sorted(set([0, n-1, 2*n-2] + top + rng.sample(range(len(coeffs)), min(8,len(coeffs))))) for k in indices: assert direct_coeff(a, b, k) == coeffs[k] if problem != "ac1": h = [Fraction(n*x, 2*s) for x in a] k = top[0] shift = n-1-k direct_rational = Fraction(2,n)*sum( (h[i]*(1-h[i+shift]) for i in range(max(0,-shift), min(n,n-shift))), Fraction(0)) assert direct_rational == score result["peak_direct_fraction_check"] = "pass" result.update({"exact_score": rational(score), "fraction_sha256": hashlib.sha256(str(score).encode()).hexdigest(), "max_coefficient": str(peak), "max_indices": [k for k,c in enumerate(coeffs) if c == peak], "packing": method, "direct_checked_indices": indices, "direct_coefficient_checks": "pass", "elapsed_seconds": time.monotonic()-started}) return result def self_test(): rng = random.Random(816) cases = [([0], [0]), ([255], [1]), ([256], [1]), ([1,1], [1,1]), ([1,2,3], [4,0,5]), ([2**100,0,2**99], [0,2**80])] for _ in range(24): cases.append(([rng.randrange(10**10) for _ in range(rng.randrange(1,35))], [rng.randrange(10**10) for _ in range(rng.randrange(1,35))])) for a,b in cases: got,_ = exact_convolution(a,b) assert got == [direct_coeff(a,b,k) for k in range(len(a)+len(b)-1)] return {"status": "pass", "cases": len(cases), "reference": "schoolbook integer convolution"} if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument("--claims", type=Path, default=Path("claims")) parser.add_argument("--out", type=Path, required=True) args = parser.parse_args() results = {"derivation": "C1=2n max(a*a)/S^2; Erdos=max(a*reverse(2S-na))/(2S^2)", "method": "carry-free Kronecker substitution, GMP exact integer multiplication", "gmp_version": C.c_char_p.in_dll(GMP, "__gmp_version").value.decode(), "self_test": self_test(), "claims": []} for name, problem in [("ac1-SX-A43.txt","ac1"), ("erdos-SX-A51.txt","erdos"), ("ac1-C-A6.txt","ac1"), ("erdos-G-A2.txt","erdos"), ("erdos-SX-A16.txt","erdos")]: item = audit(args.claims/name, problem) results["claims"].append(item) print(name, item["exact_score"]["decimal"], "PASS", flush=True) args.out.write_text(json.dumps(results, indent=2)+"\n")