diff --git a/research/dyadic-affine-feasibility/README.md b/research/dyadic-affine-feasibility/README.md new file mode 100644 index 0000000000..1418d80b45 --- /dev/null +++ b/research/dyadic-affine-feasibility/README.md @@ -0,0 +1,215 @@ +# Dyadic affine feasibility beyond a finite precision cap + +This contribution was produced by Nick Petalas. Apache-2.0. + +Claim type: scoped algebraic proof with exact finite corroboration. No new kappa. +The result distinguishes a permanent dyadic obstruction from insufficient binary +precision for a rational linear parameter region with supplied affine equalities +and a strict rational witness. It also constructs a feasible point at an explicit +sufficient precision, without enumerating that grid. + +Ali Romman's [#5794](https://github.com/CrocSwap/integer-mult-bounds/pull/5794), +head `b67e65c5fcc313f3fcb1e92b199025a63cc1dec7`, constructs minimum-precision +points by bounded grid search. Its `infeasible_through_limit` deliberately leaves +finer grids open. The new step here is an all-precision decision for the supplied +strict affine model. The predecessor's search, receipts and example are cited, +not replayed or independently reviewed. + +## The exact criterion and sufficient precision + +Let E be an integer m by n matrix, b an integer column, A a rational k by n +matrix, and h a rational column. Consider + +```text +E x = b, A x < h. +``` + +All inequalities are strict. Supply a rational witness x* satisfying both, +and integer unimodular matrices U,V giving a Smith decomposition + +```text +U E V = D = diag(d_1,...,d_r,0,...), d_i > 0, +c = U b, y = V^-1 x. +``` + +The diagonal can be rectangular. The witness ensures c_i=0 for i>r. Rational +equalities can first be cleared of denominators row by row. A precision-p point +means x belongs to 2^(-p) Z^n, with one nonnegative common precision p. + +The region contains a dyadic point at some precision **if and only if** + +```text +oddpart(d_i) divides c_i, i=1,...,r. (1) +``` + +When (1) holds, put + +```text +p_eq = max(0, max_i [v2(d_i) - v2(c_i)]). (2) +``` + +A zero c_i contributes zero to this maximum. Equivalently, p_eq is the largest +binary denominator exponent of the reduced fractions c_i/d_i. It is the exact +minimum precision for the affine equalities alone, even after an arbitrary +unimodular change of coordinates. + +For each inequality row a_l define + +```text +delta_l = h_l - a_l x* > 0, +S_l = sum_(j=r+1..n) |(A V)_(l,j)|. +``` + +Choose any integer p>=p_eq such that + +```text +S_l < 2^(p+1) delta_l for every row l. (3) +``` + +This is an explicit finite bound. Let y*=V^-1 x*. Keep y_i=c_i/d_i for i<=r +and round every free coordinate y*_j to the nearest multiple of 2^(-p), with +halfway cases toward positive infinity. Then x=V y is a feasible precision-p +point. The checker computes the least p satisfying this sufficient bound by +integer bit lengths and one exact comparison per row. It does not claim that +this p minimizes precision for the inequalities. + +## Proof + +An integer unimodular matrix and its inverse both preserve Z^n. Therefore +x is in 2^(-p) Z^n exactly when y=V^-1 x is in that same grid. The transformed +equalities fix y_i=c_i/d_i for i<=r and leave the other coordinates free. +A fixed rational coordinate belongs to some dyadic grid exactly when its +reduced denominator is a power of two. This is precisely (1), and the largest +required exponent is (2). If (1) fails, the offending coordinate excludes every +precision, independently of the inequalities. If it holds, choosing all free +coordinates zero attains p_eq for the equalities alone. + +For the strict region, use the rounded free coordinates described above. Their +individual errors are at most 2^(-p-1); the fixed coordinates have zero error. +Thus, row by row, + +```text +|a_l (x-x*)| <= 2^(-p-1) S_l < delta_l. +``` + +Consequently a_l x 0 for x in diagonal[:rank]) and + all(x == 0 for x in diagonal[rank:]), "positive Smith prefix") + require(all(diagonal[i + 1] % diagonal[i] == 0 + for i in range(rank - 1)), "Smith divisibility") + return diagonal, apply(u, b), rank, inv_v + + +def dyadic_bits(x): + denominator = x.denominator + if denominator & (denominator - 1): + return None + return denominator.bit_length() - 1 + + +def strict_power_bits(ratio): + if ratio <= 0: + return 0 + bits = max(0, ratio.numerator.bit_length() - ratio.denominator.bit_length()) + return bits + int(F(2**bits) <= ratio) + + +def construct(e, b, u, v, witness, rows=(), bounds=()): + diagonal, c, rank, inv_v = smith_data(e, b, u, v) + n = len(v) + require(len(witness) == n and len(rows) == len(bounds) <= 32, "region shape") + require(all(len(row) == n for row in rows), "inequality shape") + values = (*witness, *bounds, *(x for row in rows for x in row)) + require(all(type(x) in (int, F) and abs(F(x).numerator).bit_length() <= 64 + and F(x).denominator.bit_length() <= 64 for x in values), + "rational entry bound") + require(apply(e, witness) == tuple(b), "false affine witness") + margins = tuple(bound - value for bound, value in + zip(bounds, apply(rows, witness), strict=True)) + require(all(x > 0 for x in margins), "witness must be strict") + fixed = tuple(F(c[i], diagonal[i]) for i in range(rank)) + for i, value in enumerate(fixed): + if dyadic_bits(value) is None: + return {"status": "obstructed", "index": i, + "forced_value": str(value)} + bits = max((dyadic_bits(x) for x in fixed), default=0) + equality_bits = bits + transformed_rows = multiply(rows, v) + for row, margin in zip(transformed_rows, margins, strict=True): + seminorm = sum(abs(x) for x in row[rank:]) + bits = max(bits, strict_power_bits(F(seminorm, 2 * margin))) + require(bits <= 4096, "executable precision bound exhausted") + y = apply(inv_v, witness) + scale = 2**bits + rounded = fixed + tuple(F((scale * z + F(1, 2)).__floor__(), scale) + for z in y[rank:]) + point = apply(v, rounded) + require(apply(e, point) == tuple(b), "constructed equality failed") + require(all(value < bound for value, bound in + zip(apply(rows, point), bounds, strict=True)), "constructed margin failed") + return {"status": "constructed", "equality_bits": equality_bits, + "bits": bits, "point": point} + + +def charts(n): + result = [identity(n)] + for i in range(n): + for j in range(n): + if i != j: + for sign in (-1, 1): + a = [list(row) for row in identity(n)] + a[i][j] = sign + result.append(tuple(map(tuple, a))) + a = list(identity(n)) + a[0], a[-1] = a[-1], a[0] + result.append(tuple(a)) + return tuple(dict.fromkeys(result)) + + +def integral(a): + require(all(F(x).denominator == 1 for row in a for x in row), "not integral") + return tuple(tuple(int(x) for x in row) for row in a) + + +def coupled_cases(): + counts = {"coupled_regions": 0, "constructed": 0, "obstructed": 0} + for n in (2, 3): + profiles = ((0,) * n, (1,) + (0,) * (n - 1), + (2,) + (0,) * (n - 1), (3,) + (0,) * (n - 1), + (1,) * (n - 1) + (6,), (2,) * n) + for u, v, ds, sign, denominator in itertools.product( + charts(n), charts(n), profiles, (-1, 1), (3, 17)): + inv_u, _ = inverse(u) + inv_v, _ = inverse(v) + d = tuple(tuple(ds[i] if i == j else 0 for j in range(n)) + for i in range(n)) + e = integral(multiply(multiply(inv_u, d), inv_v)) + c = tuple(sign if x else 0 for x in ds) + b = tuple(int(x) for x in apply(inv_u, c)) + y = tuple(F(sign, x) if x else F(sign, denominator) for x in ds) + witness = apply(v, y) + rows = tuple(identity(n)) + tuple(tuple(-x for x in row) + for row in identity(n)) + margins = (F(1, denominator),) * (2 * n) + bounds = tuple(x + margin for x, margin in + zip(apply(rows, witness), margins, strict=True)) + got = construct(e, b, u, v, witness, rows, bounds) + expected_obstruction = any(F(sign, x).denominator % 3 == 0 for x in ds if x) + require((got["status"] == "obstructed") == expected_obstruction, + "wrong coupled classification") + if got["status"] == "constructed": + require(all((x * 2**got["bits"]).denominator == 1 + for x in got["point"]), "point not on grid") + expected_bits = max((F(sign, x).denominator.bit_length() - 1 + for x in ds if x), default=0) + require(got["equality_bits"] == expected_bits, "wrong equality precision") + counts["constructed"] += 1 + else: + counts["obstructed"] += 1 + counts["coupled_regions"] += 1 + return counts + + +def exhaustive_grids(): + cases = points = 0 + for v, ds, sign in itertools.product(charts(2), ((0, 0), (1, 0), (2, 0), + (3, 0), (1, 2), (2, 4)), (-1, 1)): + inv_v, _ = inverse(v) + d = tuple(tuple(ds[i] if i == j else 0 for j in range(2)) for i in range(2)) + e = integral(multiply(d, inv_v)) + b = tuple(sign if x else 0 for x in ds) + for bits in range(3): + found = False + for numerator in itertools.product(range(-8, 9), repeat=2): + point = tuple(F(x, 2**bits) for x in numerator) + found |= apply(e, point) == b + points += 1 + expected = all((sign * 2**bits) % x == 0 for x in ds if x) + require(found == expected, "literal grid disagrees with equality threshold") + cases += 1 + return {"complete_grid_cases": cases, "grid_points": points} + + +def scalar_cases(): + count = 0 + for a, b in itertools.product(range(1, 33), range(-32, 33)): + got = construct(((a,),), (b,), ((1,),), ((1,),), (F(b, a),)) + denominator = F(b, a).denominator + odd = denominator + for _ in range(denominator.bit_length()): + if odd % 2 == 0: + odd //= 2 + require((got["status"] == "obstructed") == (odd != 1), "scalar obstruction") + if odd == 1: + require(got["bits"] == denominator.bit_length() - 1, "scalar minimum") + count += 1 + return {"scalar_equalities": count} + + +def controls(): + valid = [((2,),), (1,), ((1,),), ((1,),), (F(1, 2),)] + invalid = [] + for position, value in ((0, ((3,),)), (1, (2,)), (2, ((2,),)), + (3, ((2,),)), (4, (F(1, 3),))): + case = list(valid) + case[position] = value + invalid.append(case) + invalid.extend([ + [((1, 1), (0, 1)), (1, 1), identity(2), identity(2), (0, 1)], + [((2, 0), (0, 3)), (2, 3), identity(2), identity(2), (1, 1)], + [((0,),), (1,), ((1,),), ((1,),), (0,)], + [((True,),), (1,), ((1,),), ((1,),), (1,)], + ]) + rejected = 0 + for case in invalid: + try: + construct(*case) + except ValueError: + rejected += 1 + else: + raise ValueError("invalid input accepted") + for bound in (F(1, 2), F(0)): + try: + construct(*valid, rows=((1,),), bounds=(bound,)) + except ValueError: + rejected += 1 + else: + raise ValueError("nonstrict witness accepted") + edge = construct(((0,),), (0,), ((1,),), ((1,),), (F(1, 2),), ((1,),), (1,)) + require(edge["bits"] == 1 and edge["point"] == (F(1, 2),), "strict boundary") + require(F((F(1, 2) + F(1, 2)).__floor__()) == 1, "midpoint witness") + high = construct(((2**60,),), (1,), ((1,),), ((1,),), (F(1, 2**60),)) + require(high["equality_bits"] == 60 and high["bits"] == 60, "high precision") + coupled = construct(((3, 1),), (1,), ((1,),), ((0, 1), (1, -3)), + (F(1, 3), 0), ((1, 0), (-1, 0)), (F(1, 2), F(-1, 4))) + require(coupled["status"] == "constructed", "coupled cancellation") + return {"rejected_inputs": rejected, "semantic_controls": 3} + + +def main(): + require(len(sys.argv) == 1, "this standalone replay takes no arguments") + result = {"status": "PASS_DYADIC_AFFINE_FEASIBILITY"} + result.update(scalar_cases()) + result.update(coupled_cases()) + result.update(exhaustive_grids()) + result.update(controls()) + package = Path(__file__).resolve().parent + result["sha256"] = {name: hashlib.sha256((package / name).read_bytes()).hexdigest() + for name in ("README.md", "TASK.json", "SOURCES.json", + "VERIFY.json", "verify.py")} + print(json.dumps(result, sort_keys=True)) + + +if __name__ == "__main__": + main()