diff --git a/research/odd-budget-helper-capacity/PROOF.md b/research/odd-budget-helper-capacity/PROOF.md new file mode 100644 index 0000000000..72debc5fd2 --- /dev/null +++ b/research/odd-budget-helper-capacity/PROOF.md @@ -0,0 +1,129 @@ +# Exact helper capacity with prescribed odd denominators + +This scoped algebraic result was produced by Nick Petalas. Apache-2.0. + +## Contract and theorem + +Let D=Z[1/2]. Let K be any finite m-by-n matrix over D. Fix positive odd +integers a,b. A t-channel boundary is a pair of rational matrices + +```text +X in Mat(t,m,Q), Y in Mat(n,t,Q), X K Y = I_t, +a X over D, b Y over D. +``` + +These are denominator divisibility budgets. They do not bound coefficient +magnitudes, dyadic precision, operation count or native availability. +In particular, a bound on the numerical size of a denominator is a different +contract. The t=0 boundary is allowed. + +Choose a normalized dyadic Smith chart + +```text +U K V = diag(h_1,...,h_r,0), +U in GL_m(D), V in GL_n(D), +0 < h_1 | h_2 | ... | h_r, all h_i odd. +``` + +Such a chart exists by the classical integer Smith theorem after clearing a +power-of-two denominator and absorbing powers of two into U. The factors h_i +are the positive odd parts of the ordinary integer Smith factors. + +**The maximum t is the number of h_i dividing ab.** By divisibility order +these are an initial segment. Equivalently, for 1 <= t <= r a t-channel +boundary exists exactly when h_t divides ab. For t>r it is impossible. + +Consequently, for 1 <= t <= r, the least positive odd q permitting qX and +qY over D is + +```text +q_t = product over odd primes p of p^ceil(v_p(h_t)/2). +``` + +More strongly, every feasible common q is a multiple of q_t. For t=0 use +q_0=1; for t>r no q works. With a dyadic suffix only, b=1, the exact +one-sided denominator requirement is h_t | a. With a=b=1 the theorem +reduces to Ali Romman's #5615 criterion. That case and Smith reduction +are credited predecessors, not new results here. + +## Attainment + +Let h_t divide ab and take the first t rows and columns of the Smith chart. +Define + +```text +X[i,:] = (b/h_i) U[i,:], +Y[:,i] = V[:,i]/b, 1 <= i <= t. +``` + +Then XKY=I_t entry by entry. Also aX has rows (ab/h_i)U[i,:], hence is +dyadic, and bY consists of columns of V. This constructs all channels +simultaneously. It needs no unproved lifting of a residue-ring solution. +The symmetric denominator formula follows by applying this construction +with a=b=q_t. + +## Necessity, including rectangular maps and prime powers + +Write A=aX and B=bY. The endpoint implies AKB=ab I_t. Substitute the +Smith chart and absorb its dyadic inverse factors into A and B. It suffices +to consider + +```text +A' S B' = ab I_t, S=diag(h_1,...,h_r,0), +A',B' over D. +``` + +The rational rank of S gives t<=r. Suppose h_t does not divide ab. Choose +an odd prime p with v_p(h_t)>e=v_p(ab), and let k be the number of h_i +with v_p(h_i)<=e. Divisibility order gives k None: + if not condition: + raise ValueError(message) + + +def matrix(rows: object) -> Matrix: + require(isinstance(rows, (list, tuple)), "matrix must be a list or tuple") + require(1 <= len(rows) <= MAX_DIMENSION, "matrix row capacity") + require(all(isinstance(row, (list, tuple)) for row in rows), "matrix rows") + width = len(rows[0]) + require(1 <= width <= MAX_DIMENSION, "matrix column capacity") + require(all(len(row) == width for row in rows), "ragged matrix") + values = [] + for row in rows: + converted = [] + for value in row: + require(type(value) in (int, Fraction), "exact rational entry required") + value = Fraction(value) + require(value.denominator & (value.denominator - 1) == 0, + "entry must be dyadic") + require(max(abs(value.numerator).bit_length(), + value.denominator.bit_length()) <= MAX_BITS, "entry capacity") + converted.append(value) + values.append(tuple(converted)) + return tuple(values) + + +def identity(size: int) -> Matrix: + require(type(size) is int and 0 <= size <= MAX_DIMENSION, "identity capacity") + return tuple(tuple(Fraction(i == j) for j in range(size)) for i in range(size)) + + +def multiply(left: Matrix, right: Matrix) -> Matrix: + require(bool(left) and bool(right), "nonempty matrix product required") + require(len(left[0]) == len(right), "matrix product shape") + return tuple(tuple(sum((left[i][k] * right[k][j] for k in range(len(right))), + Fraction(0)) for j in range(len(right[0]))) + for i in range(len(left))) + + +def determinant(value: Matrix) -> Fraction: + require(len(value) == len(value[0]), "square determinant required") + work = [list(row) for row in value] + result = Fraction(1) + for column in range(len(work)): + pivot = next((i for i in range(column, len(work)) if work[i][column]), None) + if pivot is None: + return Fraction(0) + if pivot != column: + work[pivot], work[column] = work[column], work[pivot] + result = -result + leading = work[column][column] + result *= leading + for i in range(column + 1, len(work)): + ratio = work[i][column] / leading + for j in range(column + 1, len(work)): + work[i][j] -= ratio * work[column][j] + return result + + +def is_unit(value: Fraction) -> bool: + numerator = abs(value.numerator) + denominator = value.denominator + return (numerator > 0 and numerator & (numerator - 1) == 0 + and denominator & (denominator - 1) == 0) + + +def validate_chart(k: Matrix, u: Matrix, v: Matrix, smith: tuple[int, ...]) -> None: + for value in (k, u, v): + require(matrix(value) == value, "invalid dyadic matrix") + m, n = len(k), len(k[0]) + require(len(u) == m and len(u[0]) == m, "left chart shape") + require(len(v) == n and len(v[0]) == n, "right chart shape") + require(type(smith) is tuple and len(smith) <= min(m, n), "Smith shape") + for i, factor in enumerate(smith): + require(type(factor) is int and 0 < factor <= MAX_ODD and factor % 2 == 1, + "Smith factors must be bounded positive odd integers") + require(i == 0 or factor % smith[i - 1] == 0, "Smith divisibility") + require(is_unit(determinant(u)) and is_unit(determinant(v)), "nonunit chart") + expected = tuple(tuple(Fraction(smith[i] if i == j and i < len(smith) else 0) + for j in range(n)) for i in range(m)) + require(multiply(multiply(u, k), v) == expected, "false Smith identity") + + +def budget(value: int) -> None: + require(type(value) is int and 0 < value <= MAX_ODD and value % 2 == 1, + "budget must be a bounded positive odd integer") + + +def factors(value: int) -> tuple[tuple[int, int], ...]: + budget(value) + remaining = value + result = [] + for prime in range(3, isqrt(value) + 1, 2): + exponent = 0 + for _ in range(20): + if remaining % prime: + break + remaining //= prime + exponent += 1 + if exponent: + result.append((prime, exponent)) + if remaining > 1: + result.append((remaining, 1)) + return tuple(result) + + +def valuation(value: int, prime: int) -> int: + require(type(value) is int and 0 < value < 1 << 40, "valuation input") + require(type(prime) is int and 3 <= prime <= MAX_ODD, "valuation divisor") + for exponent in range(41): + if value % prime: + return exponent + value //= prime + raise ValueError("valuation capacity") + + +def construct(k: Matrix, u: Matrix, v: Matrix, smith: tuple[int, ...], + a: int, b: int) -> tuple[Matrix, Matrix, int, tuple[int, int] | None]: + validate_chart(k, u, v, smith) + budget(a) + budget(b) + count = sum(a * b % factor == 0 for factor in smith) + x = tuple(tuple(Fraction(b, smith[i]) * value for value in u[i]) + for i in range(count)) + y = tuple(tuple(v[j][i] / b for i in range(count)) for j in range(len(v))) + obstruction = None + if count < len(smith): + for prime, exponent in factors(smith[count]): + if exponent > valuation(a * b, prime): + obstruction = (prime, valuation(a * b, prime)) + break + require(obstruction is not None, "missing prime obstruction") + return x, y, count, obstruction + + +def common_denominator(smith: tuple[int, ...], count: int) -> int | None: + require(type(count) is int and 0 <= count <= MAX_DIMENSION + 1, + "channel count capacity") + require(type(smith) is tuple and len(smith) <= MAX_DIMENSION, "Smith capacity") + for i, factor in enumerate(smith): + budget(factor) + require(i == 0 or factor % smith[i - 1] == 0, "Smith divisibility") + if count > len(smith): + return None + result = 1 + if count: + for prime, exponent in factors(smith[count - 1]): + result *= prime ** ((exponent + 1) // 2) + return result diff --git a/research/odd-budget-helper-capacity/verify.py b/research/odd-budget-helper-capacity/verify.py new file mode 100644 index 0000000000..2af7785b44 --- /dev/null +++ b/research/odd-budget-helper-capacity/verify.py @@ -0,0 +1,269 @@ +"""Independent finite corroboration of the new odd-budget boundary contract. + +Produced by Nick Petalas. Apache-2.0. Writes only its computed stdout receipt. +""" + +import hashlib +import itertools +import json +from fractions import Fraction +from pathlib import Path + +import capacity as candidate + + +def check(condition: bool, message: str) -> None: + if not condition: + raise AssertionError(message) + + +def inverse(value: candidate.Matrix) -> candidate.Matrix: + size = len(value) + work = [list(row) + list(candidate.identity(size)[i]) for i, row in enumerate(value)] + for column in range(size): + pivot = next(i for i in range(column, size) if work[i][column]) + work[column], work[pivot] = work[pivot], work[column] + scale = work[column][column] + work[column] = [entry / scale for entry in work[column]] + for row in range(size): + if row != column: + scale = work[row][column] + work[row] = [a - scale * b for a, b in zip(work[row], work[column], strict=True)] + check(all(work[i][:size] == list(candidate.identity(size)[i]) for i in range(size)), + "inverse elimination failed") + return tuple(tuple(row[size:]) for row in work) + + +def chart(smith: tuple[int, ...], serial: int) -> tuple[candidate.Matrix, ...]: + m, n = len(smith) + 1, len(smith) + 2 + changes = [] + for size in (m, n): + work = [list(row) for row in candidate.identity(size)] + for index in range(size - 1): + coefficient = (serial + index) % 5 - 2 + work[index] = [a + coefficient * b + for a, b in zip(work[index], work[index + 1], strict=True)] + work[-1] = [entry * Fraction(1, 2) for entry in work[-1]] + changes.append(tuple(tuple(row) for row in work)) + u, v = changes + diagonal = tuple(tuple(Fraction(smith[i] if i == j and i < len(smith) else 0) + for j in range(n)) for i in range(m)) + k = candidate.multiply(candidate.multiply(inverse(u), diagonal), inverse(v)) + return k, u, v + + +def witness_ok(k: candidate.Matrix, x: candidate.Matrix, y: candidate.Matrix, + count: int, a: int, b: int) -> bool: + if len(x) != count or len(y) != len(k[0]): + return False + if any(len(row) != len(k) for row in x) or any(len(row) != count for row in y): + return False + for rows, scale in ((x, a), (y, b)): + for row in rows: + for entry in row: + denominator = (scale * entry).denominator + if denominator & (denominator - 1): + return False + for i in range(count): + for j in range(count): + entry = sum((x[i][s] * k[s][t] * y[t][j] + for s in range(len(k)) for t in range(len(k[0]))), Fraction(0)) + if entry != int(i == j): + return False + return True + + +def chart_cases() -> dict[str, int | list[int]]: + profiles = [()] + profiles.extend((a,) for a in (1, 3, 5)) + profiles.extend((a, a * b) for a in (1, 3, 5) for b in (1, 3, 5, 9)) + profiles.extend((a, a * b, a * b * c) for a in (1, 3, 5) + for b in (1, 3, 5, 9) for c in (1, 3, 5, 9)) + profiles.append((1, 27, 675)) + cases = 0 + obstructions = 0 + for serial, smith in enumerate(profiles): + k, u, v = chart(smith, serial) + for a, b in itertools.product((1, 3, 5, 9, 15, 25), repeat=2): + x, y, count, obstruction = candidate.construct(k, u, v, smith, a, b) + check(witness_ok(k, x, y, count, a, b), "non-diagonal boundary witness") + if obstruction: + prime, exponent = obstruction + local_count = sum(h % (prime ** (exponent + 1)) != 0 for h in smith) + check(local_count == count, "prime-local upper certificate") + obstructions += 1 + else: + check(count == len(smith), "rational-rank upper certificate") + cases += 1 + for count in range(len(smith) + 1): + q = candidate.common_denominator(smith, count) + x, y, achieved, _ = candidate.construct(k, u, v, smith, q, q) + check(achieved >= count and witness_ok(k, x, y, achieved, q, q), + "least common denominator attainer") + if count: + for prime, _ in candidate.factors(q): + _, _, fewer, _ = candidate.construct(k, u, v, smith, q // prime, q // prime) + check(fewer < count, "common denominator prime is indispensable") + check(candidate.common_denominator(smith, len(smith) + 1) is None, + "above-rank common denominator") + minima = [candidate.common_denominator((1, 27, 675), t) for t in (1, 2, 3)] + check(minima == [1, 9, 45], "prime-power example") + return {"profiles": len(profiles), "budget_cases": cases, + "prime_obstructions": obstructions, "example_minima": minima} + + +def factor_exists(k: tuple[tuple[int, int], tuple[int, int]], scalar: int, + width: int, modulus: int) -> tuple[bool, int, int]: + """Brute-force every right factor, then every possible left row modulo 9.""" + right_count = 0 + row_count = 0 + targets = [tuple(scalar * int(i == j) % modulus for j in range(width)) + for i in range(width)] + for flat in itertools.product(range(modulus), repeat=2 * width): + right_count += 1 + right = (flat[:width], flat[width:]) + middle = tuple(tuple(sum(k[i][s] * right[s][j] for s in range(2)) % modulus + for j in range(width)) for i in range(2)) + image = set() + for left in itertools.product(range(modulus), repeat=2): + image.add(tuple(sum(left[s] * middle[s][j] for s in range(2)) % modulus + for j in range(width))) + row_count += 1 + if all(target in image for target in targets): + return True, right_count, row_count + return False, right_count, row_count + + +def residue_cases() -> dict[str, int]: + systems = 0 + right_count = 0 + row_count = 0 + for d0, d1 in itertools.combinations_with_replacement((1, 3, 9), 2): + for changed in range(3): + k = ((d0, changed * d0), ((changed + 1) * d0, + (changed + 1) * changed * d0 + d1)) + for scalar in (1, 3): + expected = sum(scalar % d == 0 for d in (d0, d1)) + for width in (1, 2): + exists, rights, rows = factor_exists(k, scalar, width, 9) + check(exists == (width <= expected), "independent residue-ring existence") + systems += 1 + right_count += rights + row_count += rows + check(not factor_exists(((1, 0), (0, 9)), 3, 2, 9)[0], + "determinant-only false positive") + return {"systems": systems, "right_factors_enumerated": right_count, + "left_rows_enumerated": row_count, "determinant_control": 1} + + +def scalar_cases() -> dict[str, int]: + cases = 0 + pairs = 0 + for prime in (3, 5): + for power in range(5): + for exponent in range(3): + modulus = prime ** (exponent + 1) + scalar = prime ** exponent + found = False + for left, right in itertools.product(range(modulus), repeat=2): + pairs += 1 + if (left * prime ** power * right - scalar) % modulus == 0: + found = True + check(found == (power <= exponent), "complete scalar prime-power search") + cases += 1 + return {"systems": cases, "factor_pairs_enumerated": pairs} + + +def symmetric_cases() -> dict[str, int]: + pairs = 0 + for factor in range(1, 256, 2): + feasible = [] + for q in range(1, factor + 1, 2): + if q * q % factor == 0: + feasible.append(q) + pairs += 1 + predicted = candidate.common_denominator((factor,), 1) + check(predicted == min(feasible), "independent least common denominator") + check(all(q % predicted == 0 for q in feasible), "all feasible denominators divide") + largest = 1 << 255 + k = candidate.matrix([[largest]]) + u = candidate.matrix([[Fraction(1, largest)]]) + v = candidate.identity(1) + x, y, count, _ = candidate.construct(k, u, v, (1,), candidate.MAX_ODD, candidate.MAX_ODD) + check(witness_ok(k, x, y, count, candidate.MAX_ODD, candidate.MAX_ODD), + "inclusive bit and budget capacities") + k = tuple(tuple(Fraction(i == j and i < 1) for j in range(6)) for i in range(6)) + x, y, count, _ = candidate.construct(k, candidate.identity(6), candidate.identity(6), + (1,), 1, 1) + check(witness_ok(k, x, y, count, 1, 1), "inclusive dimension capacity") + return {"odd_factors": 128, "denominator_candidates": pairs, "capacity_positives": 2} + + +def expect_rejection(function, *arguments) -> None: + try: + function(*arguments) + except (ValueError, TypeError): + return + raise AssertionError("invalid contract accepted") + + +def controls() -> dict[str, int]: + k, u, v = chart((1, 27, 675), 2) + invalid_budgets = (0, -1, 2, True, (1 << 21) + 1) + for value in invalid_budgets: + expect_rejection(candidate.construct, k, u, v, (1, 27, 675), value, 1) + for smith in ((3, 1), (0,), (2,), ((1 << 21) + 1,)): + expect_rejection(candidate.construct, k, u, v, smith, 1, 1) + altered = tuple(tuple(value + int(i == 0 and j == 0) for j, value in enumerate(row)) + for i, row in enumerate(k)) + expect_rejection(candidate.construct, altered, u, v, (1, 27, 675), 1, 1) + bad_u = tuple(tuple(value * (3 if i == 0 else 1) for value in row) + for i, row in enumerate(u)) + expect_rejection(candidate.construct, k, bad_u, v, (1, 27, 675), 1, 1) + bad_matrices = ([[Fraction(1, 3)]], [[True]], [[1], [1, 2]], [], [[1 << 256]], + [[1] * 7], [[1]] * 7) + for value in bad_matrices: + expect_rejection(candidate.matrix, value) + expect_rejection(candidate.common_denominator, (1,) * 7, 1) + expect_rejection(candidate.common_denominator, (1,), -1) + x, y, count, _ = candidate.construct(k, u, v, (1, 27, 675), 9, 75) + check(witness_ok(k, x, y, count, 9, 75), "positive witness control") + bad_x = tuple(tuple(value + int(i == 0 and j == 0) for j, value in enumerate(row)) + for i, row in enumerate(x)) + check(not witness_ok(k, bad_x, y, count, 9, 75), "changed boundary coefficient") + check(not witness_ok(k, x, y[:-1], count, 9, 75), "missing suffix row") + check(not witness_ok(k, x, y, count, 1, 75), "undeclared left odd division") + return {"rejections": len(invalid_budgets) + 4 + 2 + len(bad_matrices) + 2 + 3, + "positive_witnesses": 1} + + +def authenticate() -> dict[str, str]: + package = Path(__file__).resolve().parent + root = package.parent.parent + sources = json.loads((package / "SOURCES.json").read_text()) + hashes = {} + for row in sources["files"]: + if row["usage"] == "input": + path = root / row["path"] + digest = hashlib.sha256(path.read_bytes()).hexdigest() + check(digest == row["sha256"], "live predecessor proof drift") + hashes[row["path"]] = digest + names = ("README.md", "PROOF.md", "SAFETY.md", "TASK.json", "SOURCES.json", + "VERIFY.json", "capacity.py", "verify.py") + for name in names: + hashes[name] = hashlib.sha256((package / name).read_bytes()).hexdigest() + return hashes + + +def main() -> None: + hashes = authenticate() + receipt = {"status": "PASS_ODD_BUDGET_HELPER_CAPACITY", "source_sha256": hashes, + "charts": chart_cases(), "residue_ring": residue_cases(), + "scalar_prime_powers": scalar_cases(), "common_denominators": symmetric_cases(), + "controls": controls(), + "scope": "Scoped algebraic proof with finite corroboration; no supplier or kappa."} + print(json.dumps(receipt, indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main()