#!/usr/bin/env python3 """Check private exact KMeans primitives against CPython's Fraction oracle. Build tests/ml_exact_tests.cpp separately as a standalone C++20 executable and pass it with --binary. This tool never compiles or edits the runtime. Oracle values are exact rational interpretations of input binary64 bit patterns, and CPython's independently implemented Fraction-to-float conversion supplies the correctly rounded expected result. These are arithmetic tests, not benchmarks. """ from __future__ import annotations import argparse from collections import Counter from fractions import Fraction import hashlib import json from pathlib import Path import random import struct import subprocess import sys def fraction(bits: int) -> Fraction: exponent = (bits >> 52) & 0x7FF if exponent == 0x7FF: raise ValueError("oracle accepts finite inputs only") mantissa = bits & ((1 << 52) - 1) power = -1074 if exponent: mantissa |= 1 << 52 power = exponent - 1075 if bits >> 63: mantissa = -mantissa return Fraction(mantissa << power, 1) if power >= 0 else Fraction(mantissa, 1 << -power) def rounded_bits(value: Fraction) -> str: # Fraction.__float__ uses integer true division, independently of the C++ # accumulator and its explicit guard/sticky/remainder rounding code. try: result = float(value) except OverflowError: return "OVERFLOW" return struct.pack(">d", result).hex() def random_finite(rng: random.Random, max_exponent: int = 2046) -> int: boundaries = [0, 1, 2, 510, 511, 512, 513, 1021, 1022, 1023, 1024, 1533, 1534, 1535, 2044, 2045, 2046] candidates = [value for value in boundaries if value <= max_exponent] exponent = rng.choice(candidates) if rng.randrange(2) else rng.randrange(max_exponent + 1) fractions = [0, 1, (1 << 52) - 1, 1 << 51, (1 << 51) - 1] mantissa = rng.choice(fractions) if rng.randrange(2) else rng.getrandbits(52) return (rng.getrandbits(1) << 63) | (exponent << 52) | mantissa def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--binary", required=True, type=Path) parser.add_argument("--seed", type=int, default=20261030) parser.add_argument("--random-cases", type=int, default=2000) parser.add_argument("--output", type=Path) args = parser.parse_args() if args.random_cases < 0: parser.error("--random-cases must be nonnegative") rng = random.Random(args.seed) cases: list[tuple[str, str, str]] = [] def add_mean(name: str, values: list[int], count: int | None = None) -> None: divisor = len(values) if count is None else count expected = rounded_bits(sum((fraction(value) for value in values), Fraction()) / divisor) line = f"M {divisor} {len(values)} " + " ".join(f"{value:016x}" for value in values) cases.append((name, line, expected)) def add_squares(name: str, values: list[int], lanes: tuple[int, ...] = (1, 2, 4)) -> None: expected = rounded_bits(sum((fraction(value) ** 2 for value in values), Fraction())) for count in lanes: line = f"S {count} {len(values)} " + " ".join(f"{value:016x}" for value in values) cases.append((f"{name}/lanes={count}", line, expected)) def encode(value: float) -> int: return int.from_bytes(struct.pack(">d", value), "big") minimum_normal = 0x0010000000000000 maximum = 0x7FEFFFFFFFFFFFFF sign = 1 << 63 named_means = [ ("empty-positive-zero", [0, sign], 2), ("half-positive-subnormal", [1], 2), ("half-negative-subnormal", [sign | 1], 2), ("subnormal-even-up", [3], 2), ("subnormal-even-down", [5], 2), ("negative-subnormal-even-up", [sign | 3], 2), ("negative-subnormal-even-down", [sign | 5], 2), ("boundary-subnormal", [minimum_normal, sign | 1], 2), ("normal-even-down", [encode(1.0), encode(1.0) + 1], 2), ("normal-even-up", [encode(1.0) + 1, encode(1.0) + 2], 2), ("quotient-sticky", [encode(1.0), encode(1.0) + 1, encode(1.0) + 1], 3), ("divide-before-rounding", [maximum, maximum], 2), ("full-range-cancellation", [maximum, sign | maximum, 1], 3), ("full-range-negative-cancellation", [maximum, sign | maximum, sign | 1], 3), ("confirmed-1e16", [encode(1e16), encode(1.0), encode(-1e16)], 3), ("confirmed-1e150", [encode(1e150), encode(1.0), encode(-1e150)], 3), ("confirmed-scaled-underflow", [encode(-1e150), encode(1e150), encode(1e-200)], 3), ] for name, values, divisor in named_means: add_mean(name, values, divisor) add_mean(name + "/reverse", list(reversed(values)), divisor) add_squares("empty", [], (1, 2, 4, 7)) add_squares("confirmed-subnormal-SSE", [encode(1e-162)] * 4) add_squares("exact-minimum-subnormal", [encode(2.0 ** -537)]) add_squares("half-minimum-subnormal", [encode(2.0 ** -538)] * 2) add_squares("above-half-minimum-subnormal", [encode(2.0 ** -538)] * 3) add_squares("rounded-overflow", [maximum]) add_squares("small-after-large", [encode(1.0), encode(-1.0)] + [encode(2.0 ** -28)] * 8192) maximum_square_terms = [] for exponent in range(971, 1024): maximum_square_terms.extend([encode(2.0 ** (exponent // 2))] * (1 if exponent % 2 == 0 else 2)) quarter_ulp_terms = [encode(2.0 ** 484)] * 2 add_squares("maximum-exact-squares", maximum_square_terms, (1, 2, 4, 7)) add_squares("below-overflow-midpoint", maximum_square_terms + quarter_ulp_terms, (1, 2, 4, 7)) add_squares("overflow-midpoint-tie", maximum_square_terms + quarter_ulp_terms * 2, (1, 2, 4, 7)) # Every finite exponent field appears in both operations. Cancellation # fixtures additionally exercise sparse adds, long borrows, and signed zero # across the entire exponent range. for exponent in range(2047): bits = (exponent << 52) | rng.getrandbits(52) add_mean(f"exponent-{exponent}/single", [bits]) add_mean(f"exponent-{exponent}/cancel", [bits, sign | bits, 1], 3) add_squares(f"exponent-{exponent}/square", [bits], (1,)) for index in range(args.random_cases): count = rng.randrange(1, 65) values = [random_finite(rng) for _ in range(count)] if index % 3 == 0: base = random_finite(rng) & ~sign values = [base, sign | base] + values add_mean(f"random-{index}/mean", values) if index % 8 == 0: rng.shuffle(values) add_mean(f"random-{index}/mean-shuffled", values) # Most full-exponent square sums overflow; alternate a finite-range # distribution to exercise final rounding rather than just rejection. square_values = values if index % 2 else [random_finite(rng, 1534) for _ in range(count)] add_squares(f"random-{index}/squares", square_values) if index % 16 == 0: rng.shuffle(square_values) add_squares(f"random-{index}/squares-shuffled", square_values, (7,)) protocol = "\n".join(line for _, line, _ in cases) + "\n" run = subprocess.run([str(args.binary.resolve()), "--oracle"], input=protocol, capture_output=True, text=True, check=True, timeout=120) observed = run.stdout.splitlines() if len(observed) != len(cases): raise RuntimeError(f"oracle protocol returned {len(observed)} rows for {len(cases)} cases; stderr={run.stderr!r}") failures = [] for (name, line, expected), actual in zip(cases, observed): if actual != expected: failures.append({"name": name, "expected": expected, "actual": actual, "protocol": line}) root = Path(__file__).resolve().parent.parent hashes = {str(path.relative_to(root)): hashlib.sha256(path.read_bytes()).hexdigest() for path in [root / "src/ml_exact.hpp", root / "tests/ml_exact_tests.cpp", Path(__file__).resolve()]} report = { "schema": "aner-kmeans-exact-primitive-oracle-v1", "passed": not failures, "oracle": "Python fractions.Fraction exact rational arithmetic and CPython correctly rounded Fraction.__float__", "scope": "Mean exact sum/count and sum of squares of supplied binary64 distances; not kernel coordinate SSE or performance", "seed": args.seed, "random_cases": args.random_cases, "finite_exponent_fields_covered": [0, 2046], "case_count": len(cases), "expected_outcomes": dict(Counter("overflow" if expected == "OVERFLOW" else "finite" for _, _, expected in cases)), "binary": str(args.binary.resolve()), "binary_sha256": hashlib.sha256(args.binary.read_bytes()).hexdigest(), "source_sha256": hashes, "python": sys.version, "failure_count": len(failures), "failures": failures[:20], } if args.output: args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(json.dumps(report, indent=2) + "\n") print(json.dumps({key: report[key] for key in ["passed", "seed", "case_count", "expected_outcomes", "failure_count"]})) for failure in failures[:10]: print(json.dumps(failure), file=sys.stderr) return 1 if failures else 0 if __name__ == "__main__": raise SystemExit(main())