diff --git a/benchmarks/quadratic_precision/README.md b/benchmarks/quadratic_precision/README.md new file mode 100644 index 0000000..6511433 --- /dev/null +++ b/benchmarks/quadratic_precision/README.md @@ -0,0 +1,232 @@ +# quadratic_precision — stable entry quadratic: precision vs. cost + +Self-contained, reproducible benchmark comparing four arithmetic realisations of +the **same** asymptotic-entry algebra. Nothing here modifies production code, +the `Makefile`, or git state. + +## Question + +The production camera pre-route solves the relative-distance quadratic + +``` +F(s) = |d + q s|^2 - (R0 - rr s)^2 = a s^2 + b s + c +d = x_cur - c_frame, q = w_frame + v_frame +``` + +in `long double`, with a stable root formula, an ordinary discriminant +`b*b - 4*a*c` and entry slope, a common power-of-two scaling, an +uncertainty band, geometric validation and a numerical fallback. Which matters +for accuracy and which for speed? + +Four variants share the *entire* rest of the module (scaling, uncertainty band, +validation, fallback, dispatch) and are generated from the current +`src/asymptotic.c` by rewriting exactly one lexical region: + +| variant | coefficients | type | discriminant / slope | +|---------------------|--------------|-------------|----------------------| +| `ld_plain` | as production| `long double` | plain `b*b - 4*a*c`, plain slope | +| `ld_fma` | as production| `long double` | compensated `fmal` (experimental arm) | +| `double_fma` | ordinary `+=`| `double` | compensated `fma` | +| `double_fma_coeff` | `fma` dots, `fma` products | `double` | compensated `fma` | + +`double_fma` is the control that isolates *type* precision from *FMA*; the two +double variants isolate *coefficient accumulation*. + +## Generation and build + +`build.py` locates the kernel region between two lexical anchors in the current +working-tree `src/asymptotic.c`: + +* start: the comment `/* Long-double coefficients of the relative-distance + quadratic` +* end: the forward declaration `static void minkowski_route_entry(...)` + +It asserts each anchor and the presence of `entry_quadratic_coeffs`, +`entry_solve`, `entry_discriminant` exactly once, then emits four full-module +copies under `/generated/`. Every other line (dispatch, fallback, +`asymptotic_route_camera`, Schwarzschild path, …) is copied verbatim. If a +future edit moves an anchor or changes the region, generation **fails** rather +than silently measuring the wrong code. + +`ld_plain` retains the production kernel verbatim; `ld_fma` rewrites only the +discriminant body and the entry slope. Double +variants replace `long double`→`double`, `fmal`/`fabsl`/`fmaxl`/`frexpl`/ +`scalbnl`/`sqrtl`/`copysignl`→their double forms, `LDBL_*`→`DBL_*`, and strip +`L` suffixes from literals. A guard rejects any remaining `long double` +promotion, `L` literal or `*l`/`fmal` call in the double kernel section, so no +double expression can be re-promoted to `long double`. `double_fma_coeff` +additionally changes the coefficient accumulation to `fma`. + +The probe `#include`s the generated module, so the private static +`entry_quadratic_coeffs`/`entry_solve` and the public `asymptotic_route_camera` +that are exercised are the **actual generated code**, not an independent copy. + +Compile flags, identical for all four variants: + +``` +-std=c11 -march=native -O2 -DNDEBUG -ffp-contract=off -fopenmp +``` + +`-ffp-contract=off` prevents implicit FMA contraction from silently changing +the `ld_plain` arm; explicit `fma`/`fmal` still lower as written. + +The probe links a minimal production subset (`geodesic.c`, +`asymptotic_entry.c`, `asymptotic_schwarzschild.c`, `spacetime_common.c`, +`spacetime_minkowski.c`) plus libm. No FFTW, PNG, catalog or observer-track +data are needed. + +## Inputs + +All fixtures are frozen as IEEE-754 hex literals (`cases.json` and the +generated C header are bit-identical), including the two production grazing +rows reproduced from the Alcubierre pre-route (`x=(0,-24,0)`, `centre=2t`, +`v=2`, `R=5`, canonical `w` hex values). Families: + +* **curated / adversarial**: head-on hit/miss, clear miss, grazing + `y=nextafter(R,0)`, exact tangent, near-boundary, `1e10 + 0.5, R=1` + cancellation, large-`t` `centre=0.1t`, exact linear `a==0` (inward and + initially outward), translated-origin cancellation. +* **fixed**: axes and three oblique frames, radii + `1e-100…1e100`, `D/R ∈ {1+ulp,2,10,100,512,1024,1e4,1e8,1e10}`, impact + ratios `{0,.5,.99,1-1e-6,nextafter(1,0),1,nextafter(1,∞),1+1e-6,1.1}`. + The `512`/`1024` ratios sample the crossover where long-double coefficient + rounding meets the `128 ε_D R²` geometry tolerance (`≈512`) and the double + crossover (`≈11`), separating fallback rates from the all-UNCERTAIN `1e10` + regime. +* **moving**: `v ∈ {0,.1,2,10}` along a transverse direction. +* **growing**: `rr = -1` and `nextafter(-1,∓∞)`, **inward and + initially-outward** photons (`w = ∓tow`), axis and oblique. +* **shrinking**: `rr>0`, labelled/domain-checked algebra (radius would go + negative on the open past; roots outside `R0 - rr·s > 0` are non-physical). +* **route**: `rr = 0` plus **`rr<0` growing inward/outward** families + (positive radius on the whole past, `valid_t_min=-1e300`) across all three + direction frames, driven through the public `asymptotic_route_camera`. The + oblique-2 growing-outward family reproduces the double kernel false-MISS + cases, where a kernel MISS bypasses the fallback and the route escapes + directly; the route reference uses the actual normalised canonical `w`, so + normalisation effects are explicit. + +Exact counts are emitted by the run (and in `summary.json`); they are not +hard-coded here. + +## Reference oracle + +Inputs are exact, so they are represented as `fractions.Fraction` (Python +stdlib). `d`, `q`, `a`, `b`, `c`, the discriminant and polynomial residuals +are **exact rationals**; only `sqrt` uses `decimal` (`--precision`, default +160). A custom hex parser handles both `%a` (double) and `%La` (long double) +without `float.fromhex`, which would silently drop a 64-bit long-double mantissa +to 53 bits. Reference classification mirrors production semantics: `c<0` is +INSIDE (not an entry failure), `c==0` uses the boundary slope, `a==0` is the +exact linear branch, a real double root/tangent is MISS. + +Route references are built from the **actual canonical state the probe printed** +(`canonx*`, `canonw*`) joined with the callback samples at `t0`, so the 1-ulp +normalisation of `w` and both callback models (`centre = c0 + v(t-t0)` and the +production `centre = v t`) are handled explicitly. + +## Metrics + +Kernel: + +* status counts (ENTRY/MISS/UNCERTAIN) vs the reference, **false MISS** + (reference ENTER, kernel MISS) and **false candidate** (reference MISS, + kernel ENTRY); +* root error in ULP of the returned double root and relative error, for the + common reference-ENTER set and the intersection where *all* variants returned + ENTRY (so more UNCERTAIN cannot look better by selection); +* polynomial residual at the candidate root using the actual printed + coefficients, the exact ideal polynomial, and the nearest-double-rounded + reference root as a baseline; +* `a` exact-zero vs kernel-zero, `a` sign mismatches, `a` collapse/spurious + non-zero, and conditioning labels (`near_linear`, `late`, `ill_conditioned`). + Outward/growing and near-linear cases are ill-conditioned; a rounded-zero `a` + can turn a very late true entry into a false MISS, so these are reported + separately and are not treated as a universal precision claim. + +Route: + +* kind/status vs reference, false MISS, false candidate, camera-INSIDE + mismatches, total `INVALID/ENTRY_UNCONFIRMED` outcomes independent of reference + classification (a conservative refusal to confirm is not a false escape); +* the **public kernel classification at the actual canonical inputs** + (`kern_status` in the route CSV), so a reference ENTER / kernel MISS / route + ESCAPED (`kernel_false_miss_escaped`) is visible and is not confused with + canonical normalisation; +* fast-path vs fallback (`entry_fallback_evaluations`), `F/tol` for accepted + entries, and Pi/`L_camera` preservation. + +The reference enforces `R0 - rr·s > 0` for quadratic roots; a positive root +outside the physical radius domain is reported as `domain_clipped` and is not +counted as a physical ENTER/false-MISS. A 160-vs-240-digit consistency check is +run for every curated and near-linear/late/domain-clipped id and must report the +same status, discriminant sign and nearest-double root. + +Timing: + +* kernel: coefficient assembly + `entry_solve`, serial, `noinline` + + `volatile` sink, and an inline-asm opaque loop index so GCC cannot hoist the + pure kernel out of the repetition loop (identical for all variants). A + linearity check runs the kernel at 1× and 2× calls and reports the ratio + (expect ≈2) to prove per-call execution. +* route: public pre-route, OpenMP `static` schedule + reductions, with + route-kind and fallback counts captured so a timing gain cannot come from + silently escaping more rays. `--threads` defaults to 4. + +Codegen evidence is collected with `objdump`/`nm`: for the double variants +`entry_solve` contains hardware `vfmadd` and zero x87; for `ld_fma` it contains +x87 plus `fmal` PLT relocations (extended-precision FMA is a libm software +routine on x86-64). + +## Run + +```sh +python3 benchmarks/quadratic_precision/run.py \ + --output-dir /tmp/opencode/quadratic-comparison +``` + +Useful overrides: `--rounds N`, `--threads T`, `--precision P`, +`--kernel-target-calls N`, `--route-target-calls N`, `--skip-build` (reuse the +last build), `--repo PATH`, `--cc CC`. + +Artifacts (all under `--output-dir`): + +``` +generated/{ld_plain,ld_fma,double_fma,double_fma_coeff}.c +cases/quadratic_cases.h, cases/cases.json +environment.json, build_manifest.json +raw/kernel_.csv, raw/route_.csv +raw/kernel_metrics_.csv, raw/route_metrics_.csv +raw/curated_kernel.csv, raw/curated_route.csv +raw/microbench__r*.json, raw/routebench__r*.json +raw/fma_verify.json, raw/{linearity,curated}* +summary.json +logs/{build_commands,accuracy_*,microbench_*,routebench_*}.log +``` + +`summary.json` is the machine-readable aggregate; the script also prints the +coverage counts, per-variant ULP distributions, the curated case table, codegen +evidence and median ns/call. + +## Interpretation caveats + +* FMA improves the rounding of one product only. It **cannot** restore the + information already lost in the double inputs and in the coefficient sums + (`d = x - c`, `q = w + v`, and the `qq`/`dq`/`dd` accumulations). The data + should be read stage by stage. +* A kernel UNCERTAIN is not a failure: the shared route fallback may confirm the + entry. A kernel **MISS bypasses the fallback entirely**, so growing/outward + false-MISS families are exercised through dedicated `rr<0` route fixtures and + reported both as the raw kernel decision (`kern_status`) and the route + outcome. This benchmark does **not** claim that the fallback rescues a raw + kernel MISS on the same case; it only observes that more fallback is used to + confirm inaccurate candidates. +* The double uncertainty band differs from the long-double one by the precision + term: production uses `64·(DBL_EPSILON + LDBL_EPSILON)·scale` while the + double variants use `64·(DBL_EPSILON + DBL_EPSILON)·scale`; since + `LDBL_EPSILON ≪ DBL_EPSILON` this is roughly a factor of two in the band + width, which only matters exactly at the UNCERTAIN/MISS border. +* Number/scale results are machine- and compiler-specific; `environment.json` + records the source hash, compiler, flags and platform. +* This benchmark does not establish a universal guarantee, only the sampled + behaviour on the frozen fixture set. diff --git a/benchmarks/quadratic_precision/build.py b/benchmarks/quadratic_precision/build.py new file mode 100644 index 0000000..0dc5c35 --- /dev/null +++ b/benchmarks/quadratic_precision/build.py @@ -0,0 +1,489 @@ +#!/usr/bin/env python3 +"""Generate the four precision variants of the asymptotic entry kernel and build +one self-contained probe executable per variant. + +This module never edits the production tree. It reads the *current working +tree* ``src/asymptotic.c``, extracts the "entry kernel" region delimited by two +lexical anchors, emits four generated full-module copies under the scratch +directory, and compiles a probe that textually includes exactly one copy so the +private static ``entry_solve``/``entry_quadratic_coeffs`` are exercised as the +real generated code (not an independent toy copy). + +Variants (all share the same stable-entry algebra, scaling, uncertainty policy, +validation and fallback; only the arithmetic under test changes): + + ld_plain long double coefficients/solver, ordinary b*b-4ac and slope + ld_fma reconstructed experimental compensated fma disc/slope + double_fma double type/solver with compensated fma disc/slope, + ordinary coefficient accumulation (type-precision control) + double_fma_coeff double like double_fma plus fma coefficient accumulation + +Only the kernel region is rewritten. Every other production dispatch/fallback +path is copied verbatim. +""" +from __future__ import annotations + +import hashlib +import json +import os +import re +import shutil +import subprocess +import sys +from pathlib import Path + +# --------------------------------------------------------------------------- +# Fixed compiler configuration. Identical for all four variants so the +# comparison is not confounded by flags. ``-ffp-contract=off`` disables the +# implicit FMA contraction GCC enables at -O2 for this target, so the ld_plain +# arm cannot silently gain extra precision; explicit fma()/fmal() calls still +# lower to hardware FMA. +# --------------------------------------------------------------------------- +BASE_FLAGS = [ + "-std=c11", + "-march=native", + "-O2", + "-DNDEBUG", + "-ffp-contract=off", + "-fopenmp", +] + +VARIANTS = ("ld_plain", "ld_fma", "double_fma", "double_fma_coeff") + +# Common production translation units needed by the probe. Chosen as the +# minimal dependency closure of the entry kernel plus the public pre-route: +# geodesic initialisation, the numerical entry localizer, the Schwarzschild +# exterior referenced by the full module, the dispatch wrapper and the +# Minkowski provider. main/asymptotic/frame/catalog/PSF/FFTW are excluded. +COMMON_SOURCES = ( + "src/geodesic.c", + "src/asymptotic_entry.c", + "src/asymptotic_schwarzschild.c", + "src/spacetime_common.c", + "src/spacetime_minkowski.c", +) + +START_ANCHOR = "/* Long-double coefficients of the relative-distance quadratic" +FWD_PREFIX = "static void minkowski_route_entry(" +END_ANCHOR = ( + "static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, " + "double t0," +) + + +class BuildError(RuntimeError): + pass + + +def _replace_once(text: str, old: str, new: str, what: str) -> str: + count = text.count(old) + if count != 1: + raise BuildError(f"expected exactly one occurrence of {what}, found {count}") + return text.replace(old, new, 1) + + +def locate_kernel(lines: list[str]) -> tuple[int, int]: + """Return [start, end) line indices of the rewritten kernel region. + + start is the 'Long-double coefficients' comment; end is the forward + declaration of ``minkowski_route_entry`` (preserved verbatim). + """ + starts = [i for i, l in enumerate(lines) if l.startswith(START_ANCHOR)] + if len(starts) != 1: + raise BuildError( + f"kernel start anchor must appear exactly once, found {len(starts)}" + ) + start = starts[0] + ends = [ + i + for i, l in enumerate(lines) + if i > start and l.startswith(FWD_PREFIX) + ] + if not ends: + raise BuildError("kernel end anchor (minkowski_route_entry) not found") + end = ends[0] + if not lines[end].startswith(END_ANCHOR): + raise BuildError( + "kernel end anchor does not match the expected signature; " + "production source changed" + ) + region = "".join(lines[start:end]) + for needle in ( + "typedef struct {", + "static EntryQuadratic entry_quadratic_coeffs(", + "static EntrySolveResult entry_solve(", + "static long double entry_discriminant(", + "EntryQuadratic;", + ): + if region.count(needle) != 1: + raise BuildError( + f"kernel region must contain exactly one '{needle}', " + f"found {region.count(needle)}" + ) + # The region must not contain the forward declaration or any later route + # code: those must be byte-for-byte preserved. + if "minkowski_preroute" in region or "AsymptoticStatus asymptotic_route" in region: + raise BuildError("kernel region overran into preserved route code") + return start, end + + +def transform_ld_plain(kernel: str) -> str: + """Convert a compensated experimental kernel back to plain long double.""" + old_disc = ( + " const long double four_a = 4.0L * k->a;\n" + " const long double ac = four_a * k->c;\n" + " const long double ac_error = fmal(four_a, k->c, -ac);\n" + " *scale = k->b * k->b + fabsl(ac);\n" + " return fmal(k->b, k->b, -ac) - ac_error;\n" + ) + new_disc = ( + " const long double four_a = 4.0L * k->a;\n" + " const long double ac = four_a * k->c;\n" + " *scale = k->b * k->b + fabsl(ac);\n" + " return k->b * k->b - ac;\n" + ) + out = _replace_once(kernel, old_disc, new_disc, "ld_plain discriminant body") + old_slope = "const long double slope = fmal(2.0L * k->a, s, k->b);" + new_slope = "const long double slope = 2.0L * k->a * s + k->b;" + out = _replace_once(out, old_slope, new_slope, "ld_plain slope") + return out + + +def _assert_double_clean(t: str) -> None: + """The rewritten double kernel must contain no long-double promotion.""" + forbidden = ("long double", "fmal(", "fabsl(", "fmaxl(", "frexpl(", + "scalbnl(", "sqrtl(", "copysignl(", "LDBL_") + for tok in forbidden: + if tok in t: + raise BuildError(f"double kernel still contains '{tok}'") + if re.search(r"(?<=[0-9])L(?![A-Za-z0-9_])", t): + raise BuildError("double kernel still contains an L-suffixed literal") + + +def transform_to_double(kernel: str) -> str: + """Convert the kernel from long double to double arithmetic. + + Explicit fmal->fma, *l->*, LDBL_*->DBL_*, long double->double and strips the + L suffix from the (few) decimal literals. Explicit fma calls are retained + (double_fma keeps compensated discriminant/slope). + """ + t = kernel + for old, new in ( + ("fmal(", "fma("), + ("fabsl(", "fabs("), + ("fmaxl(", "fmax("), + ("frexpl(", "frexp("), + ("scalbnl(", "scalbn("), + ("sqrtl(", "sqrt("), + ("copysignl(", "copysign("), + ): + t = t.replace(old, new) + t = t.replace("LDBL_MIN", "DBL_MIN") + t = t.replace("LDBL_EPSILON", "DBL_EPSILON") + t = t.replace("long double", "double") + # Strip L suffixes from decimal literals; otherwise 4.0L/2.0L would + # re-promote the double expression and confound type accuracy/timing. + t = re.sub(r"(?<=[0-9])L(?![A-Za-z0-9_])", "", t) + _assert_double_clean(t) + return t + + +def transform_double_fma_coeff(kernel: str) -> str: + """double variants plus fma coefficient accumulation and products.""" + dbl = transform_to_double(kernel) + old_block = ( + " double qq = 0.0, dot_dq = 0.0, dot_dd = 0.0;\n" + " for (int i = 0; i < 3; ++i) {\n" + " qq += q[i] * q[i];\n" + " dot_dq += d[i] * q[i];\n" + " dot_dd += d[i] * d[i];\n" + " }\n" + " const double R0_ld = (double)R0;\n" + " const double rr_ld = (double)rr;\n" + " EntryQuadratic k;\n" + " k.a = qq - rr_ld * rr_ld;\n" + " k.b = 2.0 * (dot_dq + R0_ld * rr_ld);\n" + " k.c = dot_dd - R0_ld * R0_ld;\n" + ) + new_block = ( + " double qq = 0.0, dot_dq = 0.0, dot_dd = 0.0;\n" + " for (int i = 0; i < 3; ++i) {\n" + " qq = fma(q[i], q[i], qq);\n" + " dot_dq = fma(d[i], q[i], dot_dq);\n" + " dot_dd = fma(d[i], d[i], dot_dd);\n" + " }\n" + " const double R0_ld = (double)R0;\n" + " const double rr_ld = (double)rr;\n" + " EntryQuadratic k;\n" + " k.a = fma(-rr_ld, rr_ld, qq);\n" + " k.b = 2.0 * fma(R0_ld, rr_ld, dot_dq);\n" + " k.c = fma(-R0_ld, R0_ld, dot_dd);\n" + ) + return _replace_once(dbl, old_block, new_block, "double_fma_coeff block") + + +def generate_variants(repo: Path, out_dir: Path) -> dict: + """Read current src/asymptotic.c and emit the four variant modules.""" + src_path = repo / "src" / "asymptotic.c" + text = src_path.read_text() + lines = text.splitlines(keepends=True) + start, end = locate_kernel(lines) + kernel = "".join(lines[start:end]) + prefix = "".join(lines[:start]) + suffix = "".join(lines[end:]) + sha = hashlib.sha256(text.encode()).hexdigest() + + # Production now uses plain long-double arithmetic. Reconstruct the former + # compensated arm so the original four-way comparison remains reproducible. + plain_disc = ( + " const long double four_a = 4.0L * k->a;\n" + " const long double ac = four_a * k->c;\n" + " *scale = k->b * k->b + fabsl(ac);\n" + " return k->b * k->b - ac;\n" + ) + fused_disc = plain_disc.replace( + " *scale =", " const long double ac_error = fmal(four_a, k->c, -ac);\n *scale =" + ).replace("return k->b * k->b - ac;", "return fmal(k->b, k->b, -ac) - ac_error;") + fused = _replace_once(kernel, plain_disc, fused_disc, "ld_fma discriminant body") + fused = _replace_once( + fused, "const long double slope = 2.0L * k->a * s + k->b;", + "const long double slope = fmal(2.0L * k->a, s, k->b);", "ld_fma slope" + ) + kernels = { + "ld_plain": kernel, + "ld_fma": fused, + "double_fma": transform_to_double(fused), + "double_fma_coeff": transform_double_fma_coeff(fused), + } + + gen_dir = out_dir / "generated" + gen_dir.mkdir(parents=True, exist_ok=True) + paths = {} + metas = {} + for name in VARIANTS: + body = kernels[name] + banner = ( + f"/* GENERATED FILE - do not edit.\n" + f" * variant: {name}\n" + f" * source: src/asymptotic.c sha256={sha}\n" + f" * region lines [{start + 1}, {end}] rewritten; all other code verbatim.\n" + f" */\n" + ) + out = banner + prefix + body + suffix + path = gen_dir / f"{name}.c" + path.write_text(out) + paths[name] = path + metas[name] = { + "path": str(path), + "sha256": hashlib.sha256(out.encode()).hexdigest(), + "kernel_sha256": hashlib.sha256(body.encode()).hexdigest(), + } + return { + "source": str(src_path), + "source_sha256": sha, + "kernel_line_range": [start + 1, end], + "variants": metas, + } + + +def cc_version(cc: str) -> str: + try: + out = subprocess.run( + [cc, "--version"], capture_output=True, text=True, check=False + ) + return (out.stdout or out.stderr).splitlines()[0] if out.stdout or out.stderr else "" + except OSError: + return "" + + +def environment(repo: Path, cc: str, threads: int) -> dict: + src_path = repo / "src" / "asymptotic.c" + env = { + "repo": str(repo), + "cc": cc, + "cc_version": cc_version(cc), + "base_flags": list(BASE_FLAGS), + "threads": threads, + "source_sha256": hashlib.sha256(src_path.read_bytes()).hexdigest(), + "git_head": _git(repo, "rev-parse", "HEAD"), + "git_status_short": _git(repo, "status", "--short"), + "uname": os.uname().sysname + " " + os.uname().release + " " + os.uname().machine, + "cpu_model": _cpu_model(), + } + return env + + +def _git(repo: Path, *args: str) -> str: + try: + out = subprocess.run( + ["git", "-C", str(repo), *args], capture_output=True, text=True, check=False + ) + return out.stdout.strip() + except OSError: + return "" + + +def _cpu_model() -> str: + try: + for line in Path("/proc/cpuinfo").read_text().splitlines(): + if line.startswith("model name"): + return line.split(":", 1)[1].strip() + except OSError: + pass + return "unknown" + + +def compile_common(repo: Path, out_dir: Path, cc: str, flags: list[str]) -> list[Path]: + obj_dir = out_dir / "build" / "common" + obj_dir.mkdir(parents=True, exist_ok=True) + objs = [] + for rel in COMMON_SOURCES: + src = repo / rel + obj = obj_dir / (Path(rel).stem + ".o") + cmd = [cc, *flags, f"-I{repo / 'src'}", "-c", str(src), "-o", str(obj)] + _run(cmd, out_dir, f"compile common {rel}") + objs.append(obj) + return objs + + +def compile_probe( + repo: Path, + out_dir: Path, + cc: str, + flags: list[str], + variant: str, + probe_src: Path, + cases_dir: Path, + common_objs: list[Path], +) -> Path: + """Compile the probe (which #includes the generated variant) and link.""" + exe_dir = out_dir / "build" / "bin" + exe_dir.mkdir(parents=True, exist_ok=True) + obj = out_dir / "build" / f"probe_{variant}.o" + variant_src = out_dir / "generated" / f"{variant}.c" + cmd = [ + cc, + *flags, + f"-I{repo / 'src'}", + f"-I{cases_dir}", + f'-DPROBE_VARIANT_SOURCE="{variant_src}"', + f'-DPROBE_VARIANT_NAME="{variant}"', + "-c", + str(probe_src), + "-o", + str(obj), + ] + _run(cmd, out_dir, f"compile probe {variant}") + exe = exe_dir / f"probe_{variant}" + link = [cc, *flags, str(obj), *[str(o) for o in common_objs], "-lm", "-o", str(exe)] + _run(link, out_dir, f"link probe {variant}") + return exe + + +def _run(cmd: list[str], out_dir: Path, what: str) -> None: + log_dir = out_dir / "logs" + log_dir.mkdir(parents=True, exist_ok=True) + with open(log_dir / "build_commands.log", "a") as fh: + fh.write(what + "\n$ " + " ".join(cmd) + "\n") + try: + proc = subprocess.run(cmd, capture_output=True, text=True, check=False) + except OSError as exc: + raise BuildError(f"{what}: {exc}") from exc + with open(log_dir / "build_commands.log", "a") as fh: + fh.write(f"exit={proc.returncode}\n") + if proc.stdout: + fh.write(proc.stdout) + if proc.stderr: + fh.write(proc.stderr) + if proc.returncode != 0: + raise BuildError( + f"{what} failed (exit {proc.returncode}); see logs/build_commands.log\n" + + (proc.stderr or "")[-4000:] + ) + + +def build_all(repo: Path, out_dir: Path, cc: str, flags: list[str], threads: int) -> dict: + out_dir.mkdir(parents=True, exist_ok=True) + (out_dir / "logs").mkdir(parents=True, exist_ok=True) + variant_info = generate_variants(repo, out_dir) + common = compile_common(repo, out_dir, cc, flags) + probe_src = repo / "benchmarks" / "quadratic_precision" / "fixtures" / "probe.c" + cases_dir = out_dir / "cases" + exes = {} + for name in VARIANTS: + exes[name] = compile_probe( + repo, out_dir, cc, flags, name, probe_src, cases_dir, common + ) + result = {"variants": variant_info["variants"], + "source_sha256": variant_info["source_sha256"], + "base_flags": list(flags), + "dependencies": dependency_hashes(repo), + "compiler": {"command": cc, "version": cc_version(cc)}, + "probe_sha256": _file_sha(probe_src), + "cases_header_sha256": _file_sha(cases_dir / "quadratic_cases.h")} + result["executables"] = {k: str(v) for k, v in exes.items()} + result["common_objects"] = [str(o) for o in common] + (out_dir / "build_manifest.json").write_text(json.dumps(result, indent=2) + "\n") + return result + + +def _file_sha(path: Path): + if not Path(path).exists(): + return None + return hashlib.sha256(Path(path).read_bytes()).hexdigest() + + +def dependency_hashes(repo: Path): + # Include all project headers conservatively, including transitive includes. + paths = {repo / rel for rel in COMMON_SOURCES} + paths.update((repo / "src").rglob("*.h")) + return {str(p.relative_to(repo)): _file_sha(p) for p in sorted(paths)} + + +def verify_manifest(repo: Path, out_dir: Path, flags: list[str], cc: str): + """Check that an existing build matches the current source, flags and + fixtures. Returns (manifest, list_of_mismatch_reasons). Regenerates the + variant sources as a side effect so their hashes can be compared.""" + manifest = json.loads((out_dir / "build_manifest.json").read_text()) + reasons = [] + if manifest.get("dependencies") != dependency_hashes(repo): + reasons.append("dependencies") + if manifest.get("compiler") != {"command": cc, "version": cc_version(cc)}: + reasons.append("compiler") + current = generate_variants(repo, out_dir) + if manifest.get("source_sha256") != current["source_sha256"]: + reasons.append("source_sha256") + if list(manifest.get("base_flags", [])) != list(flags): + reasons.append("base_flags") + for name, meta in current["variants"].items(): + if manifest.get("variants", {}).get(name, {}).get("sha256") != meta["sha256"]: + reasons.append(f"variant:{name}") + probe_src = repo / "benchmarks" / "quadratic_precision" / "fixtures" / "probe.c" + if manifest.get("probe_sha256") != _file_sha(probe_src): + reasons.append("probe_sha256") + cases_header = out_dir / "cases" / "quadratic_cases.h" + if cases_header.exists() and manifest.get("cases_header_sha256") != _file_sha(cases_header): + reasons.append("cases_header_sha256") + for name, exe in manifest.get("executables", {}).items(): + if not Path(exe).exists(): + reasons.append(f"missing_exe:{name}") + return manifest, reasons + + +if __name__ == "__main__": + import argparse + + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--repo", default=str(Path(__file__).resolve().parents[2])) + ap.add_argument("--output-dir", default="/tmp/opencode/quadratic-comparison") + ap.add_argument("--cc", default=os.environ.get("CC", "cc")) + ap.add_argument("--threads", type=int, default=4) + args = ap.parse_args() + repo = Path(args.repo).resolve() + out = Path(args.output_dir).resolve() + try: + info = build_all(repo, out, args.cc, BASE_FLAGS, args.threads) + except BuildError as exc: + print(f"build failed: {exc}", file=sys.stderr) + sys.exit(1) + print(json.dumps(info, indent=2)) diff --git a/benchmarks/quadratic_precision/cases.py b/benchmarks/quadratic_precision/cases.py new file mode 100644 index 0000000..9ef970e --- /dev/null +++ b/benchmarks/quadratic_precision/cases.py @@ -0,0 +1,441 @@ +#!/usr/bin/env python3 +"""Deterministic fixture generation for the quadratic precision benchmark. + +All values are frozen as IEEE-754 hex literals in the generated C header and as +hex strings in cases.json, so the C probe and the Python oracle see bit-for-bit +identical inputs. No external catalog/observer/slab data is required. + +Case families (kernel): + curated explicit adversarial / production fixtures + fixed fixed sphere, axis and oblique directions + moving constant sphere velocity (radial/transverse), inward photons + growing rr < 0 (radius grows along the past parameter), both inward + (toward-center) and initially-outward photons, axis + oblique + shrinking rr > 0, kernel only (radius would go negative on the far past) + cancellation translated/far-origin style input cancellation + +Route cases use only rr == 0 (positive radius on the whole open past segment) +plus the two production grazing rows. +""" +from __future__ import annotations + +import json +import math +from pathlib import Path + +# Radii spanning subnormal-adjacent to huge scales. +RADII = [1e-100, 1e-10, 1.0, 1e10, 1e100] +DIST = [ + math.nextafter(1.0, math.inf), + 2.0, + 10.0, + 100.0, + 512.0, + 1024.0, + 1e4, + 1e8, + 1e10, +] +IMPACTS = [ + 0.0, + 0.5, + 0.99, + 1.0 - 1e-6, + math.nextafter(1.0, 0.0), + 1.0, + math.nextafter(1.0, math.inf), + 1.0 + 1e-6, + 1.1, +] +VELS = [0.0, 0.1, 2.0, 10.0] +RR_GROW = [ + -1.0, + math.nextafter(-1.0, -math.inf), + math.nextafter(-1.0, math.inf), + -0.999999, +] + + +def norm3(v): + n = math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) + return [v[0] / n, v[1] / n, v[2] / n] + + +def cross(a, b): + return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0]] + + +# Fixed direction frames: (tow, perp). `tow` points from the sphere centre to +# the camera; `w = -tow` is the inward (toward-centre) past direction. +DIRS_AXIS = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +DIRS_OBLIQUE1 = ( + norm3([1.0, 2.0, 3.0]), + norm3(cross(norm3([1.0, 2.0, 3.0]), [0.0, 0.0, 1.0])), +) +DIRS_OBLIQUE2 = ( + norm3([-2.0, 1.0, 0.7]), + norm3(cross(norm3([-2.0, 1.0, 0.7]), [0.0, 1.0, 0.0])), +) +DIRS = [DIRS_AXIS, DIRS_OBLIQUE1, DIRS_OBLIQUE2] + + +def _case(cid, category, x, c, w, v, R0, rr): + return { + "id": cid, + "category": category, + "x": list(x), + "c": list(c), + "w": list(w), + "v": list(v), + "R0": R0, + "rr": rr, + } + + +def curated_kernel_cases(): + """Explicit adversarial and production-reproduction fixtures.""" + cases = [] + t63 = float.fromhex("0x1.fa8f5c28f5c29p+3") + t64 = float.fromhex("0x1.fb17e4b17e4b1p+3") + w63 = [-float.fromhex("0x1.d4afba4704cap-2"), + float.fromhex("0x1.c7378f8e872d1p-1"), + float.fromhex("0x1.15bad4e30e8ddp-8")] + w64 = [-float.fromhex("0x1.f98ae1a782104p-2"), + float.fromhex("0x1.bd3bb364ac492p-1"), + float.fromhex("0x1.102d2a1c6ac74p-7")] + # Production Alcubierre grazing rows: x=(0,-24,0), centre=2t, v=2, R=5. + cases.append(_case(0, "real63", [0.0, -24.0, 0.0], + [2.0 * t63, 0.0, 0.0], w63, [2.0, 0.0, 0.0], 5.0, 0.0)) + cases.append(_case(1, "real64", [0.0, -24.0, 0.0], + [2.0 * t64, 0.0, 0.0], w64, [2.0, 0.0, 0.0], 5.0, 0.0)) + + D = 10.0 + R = 5.0 + cases += [ + _case(2, "headon_hit", [D, 0.0, 0.0], [0.0] * 3, [-1.0, 0.0, 0.0], + [0.0] * 3, R, 0.0), + _case(3, "headon_miss", [D, 0.0, 0.0], [0.0] * 3, [1.0, 0.0, 0.0], + [0.0] * 3, R, 0.0), + _case(4, "clear_miss", [D, 6.0, 0.0], [0.0] * 3, [-1.0, 0.0, 0.0], + [0.0] * 3, R, 0.0), + _case(5, "grazing_in", [D, math.nextafter(R, 0.0), 0.0], [0.0] * 3, + [-1.0, 0.0, 0.0], [0.0] * 3, R, 0.0), + _case(6, "exact_tangent", [D, R, 0.0], [0.0] * 3, [-1.0, 0.0, 0.0], + [0.0] * 3, R, 0.0), + _case(7, "grazing_out", [D, math.nextafter(R, math.inf), 0.0], + [0.0] * 3, [-1.0, 0.0, 0.0], [0.0] * 3, R, 0.0), + _case(8, "near_boundary_hit", + [math.nextafter(R, math.inf), 0.0, 0.0], [0.0] * 3, + [-1.0, 0.0, 0.0], [0.0] * 3, R, 0.0), + # 1e10 + 0.5, R=1: c loses the transverse term, b*b-4ac rounds to 0. + _case(9, "cancel_1e10_05", [1e10, 0.5, 0.0], [0.0] * 3, + [-1.0, 0.0, 0.0], [0.0] * 3, 1.0, 0.0), + _case(10, "cancel_1e10_1", [1e10, 1.0, 0.0], [0.0] * 3, + [-1.0, 0.0, 0.0], [0.0] * 3, 1.0, 0.0), + ] + # Large-t positive reconstructed minimum (tests/test_asymptotic.c). + camera_x = float.fromhex("0x1.6bcc41e901908p+46") + t0 = 1e15 + vx = 0.1 + cases.append(_case(11, "large_t_01t", [camera_x, 0.0, 0.0], + [vx * t0, 0.0, 0.0], [-1.0, 0.0, 0.0], [vx, 0.0, 0.0], + 10.0, 0.0)) + # Exact linear growing sphere (a == 0), inward. + cases.append(_case(12, "linear_grow_in", [100.0, 0.0, 0.0], [0.0] * 3, + [-1.0, 0.0, 0.0], [0.0] * 3, 10.0, -1.0)) + # Exact linear, initially outward: constant gap, honest miss. + cases.append(_case(13, "linear_grow_out", [100.0, 0.0, 0.0], [0.0] * 3, + [1.0, 0.0, 0.0], [0.0] * 3, 10.0, -1.0)) + # Near-linear oblique, inward vs initially outward, rr just around -1. + tow = DIRS_OBLIQUE1[0] + perp = DIRS_OBLIQUE1[1] + for rr in RR_GROW: + for sign, tag in ((-1.0, "in"), (1.0, "out")): + w = [sign * tow[i] for i in range(3)] + x = [50.0 * tow[i] + 0.25 * perp[i] for i in range(3)] + cases.append(_case(len(cases), f"nearlin_{tag}", x, [0.0] * 3, w, + [0.0] * 3, 4.0, rr)) + # Translated-origin style cancellation: d = x - c with x = c + small. + base = 1e10 + cbase = [base, -base, base * 0.5] + xb = [cbase[0] + 10.0, cbase[1] + 0.5, cbase[2] + 0.0] + cases.append(_case(len(cases), "translated_origin", xb, cbase, + [-1.0, 0.0, 0.0], [0.0] * 3, 1.0, 0.0)) + return cases + + +def _broad_fixed(cid): + for di, (tow, perp) in enumerate(DIRS): + for R in RADII: + for dr in DIST: + for ir in IMPACTS: + D = dr * R + b = ir * R + x = [D * tow[i] + b * perp[i] for i in range(3)] + w = [-tow[i] for i in range(3)] + yield _case(cid, f"fixed_d{di}", x, [0.0] * 3, w, [0.0] * 3, + R, 0.0) + cid += 1 + return cid + + +def _broad_moving(cid): + for (tow, perp) in DIRS: + for R in (1e-10, 1.0, 1e10): + for dr in (2.0, 100.0, 1e4): + for ir in (0.0, 0.99, 1.0, 1.1): + for vel in VELS: + D = dr * R + b = ir * R + x = [D * tow[i] + b * perp[i] for i in range(3)] + w = [-tow[i] for i in range(3)] + v = [vel * perp[i] for i in range(3)] + yield _case(cid, "moving", x, [0.0] * 3, w, v, R, 0.0) + cid += 1 + return cid + + +def _broad_growing(cid): + for di, (tow, perp) in enumerate(DIRS): + for sign, tag in ((-1.0, "in"), (1.0, "out")): + for rr in RR_GROW: + for dr in (2.0, 10.0, 100.0, 1e4, 1e8): + for ir in (0.0, 0.5, 0.99, 1.0, 1.1): + D = dr * 4.0 + b = ir * 4.0 + x = [D * tow[i] + b * perp[i] for i in range(3)] + w = [sign * tow[i] for i in range(3)] + yield _case(cid, f"growing_{tag}", x, [0.0] * 3, w, + [0.0] * 3, 4.0, rr) + cid += 1 + return cid + + +def _broad_shrinking(cid): + for (tow, perp) in (DIRS_AXIS, DIRS_OBLIQUE1): + for dr in (2.0, 10.0, 100.0): + for ir in (0.0, 0.99, 1.0): + D = dr * 10.0 + b = ir * 10.0 + x = [D * tow[i] + b * perp[i] for i in range(3)] + w = [-tow[i] for i in range(3)] + yield _case(cid, "shrinking", x, [0.0] * 3, w, [0.0] * 3, 10.0, + 0.1) + cid += 1 + return cid + + +def _materialize(gen, cases): + for c in gen: + cases.append(c) + + +def _build_all(): + cases = curated_kernel_cases() + cid = 1000 + for gen in (_broad_fixed, _broad_moving, _broad_growing, _broad_shrinking): + gen_cases = [] + _materialize(gen(cid), gen_cases) + if gen_cases: + cid = gen_cases[-1]["id"] + 1 + cases.extend(gen_cases) + return cases + + +# --------------------------------------------------------------------------- +def _route(cid, category, t0, obs, direction, c0, v, R0, rr, model=0, + valid_t_min=-1e300): + return { + "id": cid, + "category": category, + "t0": t0, + "obs": list(obs), + "dir": list(direction), + "c0": list(c0), + "v": list(v), + "R0": R0, + "rr": rr, + "valid_t_min": valid_t_min, + "model": model, + } + + +def route_cases(): + cases = [] + t63 = float.fromhex("0x1.fa8f5c28f5c29p+3") + t64 = float.fromhex("0x1.fb17e4b17e4b1p+3") + w63 = [-float.fromhex("0x1.d4afba4704cap-2"), + float.fromhex("0x1.c7378f8e872d1p-1"), + float.fromhex("0x1.15bad4e30e8ddp-8")] + w64 = [-float.fromhex("0x1.f98ae1a782104p-2"), + float.fromhex("0x1.bd3bb364ac492p-1"), + float.fromhex("0x1.102d2a1c6ac74p-7")] + # Production rows reproduced through the Alcubierre-style callback (model 1). + cases.append(_route(0, "real63", t63, [0.0, -24.0, 0.0], w63, + [2.0 * t63, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0, + model=1)) + cases.append(_route(1, "real64", t64, [0.0, -24.0, 0.0], w64, + [2.0 * t64, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0, + model=1)) + # Same rows through the input-stable callback (model 0). + cases.append(_route(2, "real63_stable", t63, [0.0, -24.0, 0.0], w63, + [2.0 * t63, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0, + model=0)) + cases.append(_route(3, "real64_stable", t64, [0.0, -24.0, 0.0], w64, + [2.0 * t64, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0, + model=0)) + R = 5.0 + cases += [ + _route(10, "headon_hit", 0.0, [10.0, 0.0, 0.0], [-1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + _route(11, "headon_miss", 0.0, [10.0, 0.0, 0.0], [1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + _route(12, "clear_miss", 0.0, [10.0, 6.0, 0.0], [-1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + _route(13, "grazing_in", 0.0, [10.0, math.nextafter(R, 0.0), 0.0], + [-1.0, 0.0, 0.0], [0.0] * 3, [0.0] * 3, R, 0.0), + _route(14, "exact_tangent", 0.0, [10.0, R, 0.0], [-1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + _route(15, "grazing_out", 0.0, + [10.0, math.nextafter(R, math.inf), 0.0], [-1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + _route(16, "camera_inside", 0.0, [2.0, 0.0, 0.0], [1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + _route(17, "near_boundary_hit", 0.0, + [math.nextafter(R, math.inf), 0.0, 0.0], [-1.0, 0.0, 0.0], + [0.0] * 3, [0.0] * 3, R, 0.0), + ] + # Broad fixed-sphere subset through the public route. + cid = 100 + for (tow, perp) in (DIRS_AXIS, DIRS_OBLIQUE1): + for R0 in (1e-10, 1.0, 1e10): + for dr in (2.0, 10.0, 100.0, 512.0, 1024.0, 1e4): + for ir in (0.0, 0.5, 0.99, 1.0, 1.1): + D = dr * R0 + b = ir * R0 + obs = [D * tow[i] + b * perp[i] for i in range(3)] + direction = [-tow[i] for i in range(3)] + cases.append(_route(cid, f"route_fixed_d{dr:g}", 0.0, obs, + direction, [0.0] * 3, [0.0] * 3, R0, + 0.0)) + cid += 1 + # Moving spheres (v = 2 tow), model 0. + for ir in (0.0, 0.99, 1.1): + cases.append(_route(cid, "route_moving", 0.0, [100.0, 0.0, 0.0], + [-1.0, 0.0, 0.0], [0.0] * 3, [2.0, 0.0, 0.0], 5.0, + 0.0)) + cid += 1 + # Growing worldtubes (rr <= 0, positive radius on the whole past), inward + # and initially-outward photons, axis + oblique. These mirror the + # growing_out kernel false-MISS family and exercise the public route where + # a kernel MISS bypasses the fallback entirely. valid_t_min is far below + # any sampled time so no artificial history clip is introduced. + R0 = 4.0 + for (tow, perp) in DIRS: + for sign, tag in ((-1.0, "in"), (1.0, "out")): + for rr in RR_GROW: + for dr in (2.0, 10.0, 100.0, 512.0, 1024.0, 1e4): + for ir in (0.0, 1.0, 1.1): + D = dr * R0 + b = ir * R0 + obs = [D * tow[i] + b * perp[i] for i in range(3)] + direction = [sign * tow[i] for i in range(3)] + cases.append(_route(cid, f"route_grow_{tag}_d{dr:g}", + 0.0, obs, direction, [0.0] * 3, + [0.0] * 3, R0, rr)) + cid += 1 + return cases + + +# --------------------------------------------------------------------------- +def _hex(x): + return float(x).hex() + + +def _fmt(v): + return _hex(v) + + +def write_header(kernel, route, out_path: Path): + lines = [] + lines.append("/* Generated by benchmarks/quadratic_precision/cases.py. */\n") + lines.append("#ifndef QUADRATIC_CASES_H\n#define QUADRATIC_CASES_H\n") + lines.append("typedef struct {\n int id;\n const char *category;\n" + " double x[3], c[3], w[3], v[3];\n double R0, rr;\n" + "} QuadKernelCase;\n\n") + lines.append("typedef struct {\n int id;\n const char *category;\n" + " double t0;\n double obs[3];\n double dir[3];\n" + " double c0[3];\n double v[3];\n double R0, rr;\n" + " double valid_t_min;\n int model;\n} QuadRouteCase;\n\n") + lines.append("static const QuadKernelCase quad_kernel_cases[] = {\n") + for c in kernel: + lines.append( + " {.id=%d,.category=\"%s\"," + ".x={%s,%s,%s},.c={%s,%s,%s},.w={%s,%s,%s},.v={%s,%s,%s}," + ".R0=%s,.rr=%s},\n" + % (c["id"], c["category"], *[_fmt(z) for z in c["x"]], + *[_fmt(z) for z in c["c"]], *[_fmt(z) for z in c["w"]], + *[_fmt(z) for z in c["v"]], _fmt(c["R0"]), _fmt(c["rr"]))) + lines.append("};\n") + lines.append("static const int quad_kernel_case_count = %d;\n\n" + % len(kernel)) + lines.append("static const QuadRouteCase quad_route_cases[] = {\n") + for c in route: + lines.append( + " {.id=%d,.category=\"%s\",.t0=%s," + ".obs={%s,%s,%s},.dir={%s,%s,%s},.c0={%s,%s,%s},.v={%s,%s,%s}," + ".R0=%s,.rr=%s,.valid_t_min=%s,.model=%d},\n" + % (c["id"], c["category"], _fmt(c["t0"]), + *[_fmt(z) for z in c["obs"]], *[_fmt(z) for z in c["dir"]], + *[_fmt(z) for z in c["c0"]], *[_fmt(z) for z in c["v"]], + _fmt(c["R0"]), _fmt(c["rr"]), _fmt(c["valid_t_min"]), + c["model"])) + lines.append("};\n") + lines.append("static const int quad_route_case_count = %d;\n" + % len(route)) + lines.append("#endif\n") + out_path.write_text("".join(lines)) + + +def write_json(kernel, route, out_path: Path): + def enc(c): + d = {} + for k, v in c.items(): + if isinstance(v, float): + d[k] = v.hex() + elif isinstance(v, list): + d[k] = [z.hex() for z in v] + else: + d[k] = v + return d + + out_path.write_text( + json.dumps( + {"deterministic": True, "kernel": [enc(c) for c in kernel], + "route": [enc(c) for c in route]}, + indent=1, + ) + + "\n") + + +def build(): + kernel = _build_all() + route = route_cases() + return kernel, route + + +if __name__ == "__main__": + k, r = build() + print(f"kernel cases: {len(k)}, route cases: {len(r)}") + cats = {} + for c in k: + cats[c["category"]] = cats.get(c["category"], 0) + 1 + for key in sorted(cats): + print(f" {key}: {cats[key]}") + rcats = {} + for c in r: + rcats[c["category"]] = rcats.get(c["category"], 0) + 1 + print("route categories:") + for key in sorted(rcats): + print(f" {key}: {rcats[key]}") diff --git a/benchmarks/quadratic_precision/fixtures/probe.c b/benchmarks/quadratic_precision/fixtures/probe.c new file mode 100644 index 0000000..bd3f0cf --- /dev/null +++ b/benchmarks/quadratic_precision/fixtures/probe.c @@ -0,0 +1,614 @@ +/* Quadratic precision benchmark probe. + * + * Compiled once per generated variant, with + * -DPROBE_VARIANT_SOURCE=".../ld_fma.c" -DPROBE_VARIANT_NAME="ld_fma" + * The probe #includes the generated full module, so the private static + * entry_quadratic_coeffs / entry_solve and the public asymptotic_route_camera + * are the *actual generated* code. No production file is modified. + * + * Modes: + * accuracy -- run every kernel case through the generated kernel and every + * route case through the generated public pre-route; emit CSVs. + * microbench -- timed repeated kernel assembly + entry_solve (serial). + * routebench -- timed repeated public pre-route (OpenMP, static schedule). + */ +#include PROBE_VARIANT_SOURCE + +#include +#include +#include +#include +#include +#include +#include + +#include "quadratic_cases.h" + +#ifndef PROBE_VARIANT_NAME +#define PROBE_VARIANT_NAME "unknown" +#endif + +/* ------------------------------------------------------------------ */ +/* Exact-value coefficient printing: long double uses %La (full mantissa), + * double uses %a. _Generic selects the right printer without knowing the + * generated EntryQuadratic field type. */ +static void coeff_ld(FILE *f, long double v) { fprintf(f, "%La", v); } +static void coeff_d(FILE *f, double v) { fprintf(f, "%a", v); } +#define PRINT_COEFF(f, v) \ + _Generic((v), long double : coeff_ld, double : coeff_d)((f), (v)) + +/* ------------------------------------------------------------------ */ +/* Fixture source: flat identity metric, one Minkowski end whose worldtube is + * a constant-velocity (optionally linearly growing/shrinking) sphere. + * + * model 0 (input-stable): center = c0 + v*(t - t0), R = R0 + rr*(t - t0) + * model 1 (production-style): center = v*t, R = R0 + rr*t + * + * model 0 evaluates exactly at the camera: center(t0) == c0, R(t0) == R0, so + * the kernel inputs are the frozen case values. model 1 mirrors the + * Alcubierre-style callback (used only with rr == 0). */ +typedef struct { + double c0[3]; + double v[3]; + double R0; + double rr; + double t0; + double valid_t_min; + int model; +} FixtureContext; + +static double fixture_center(const FixtureContext *ctx, int i, double t) { + if (ctx->model == 0) + return ctx->c0[i] + ctx->v[i] * (t - ctx->t0); + return ctx->v[i] * t; +} + +static double fixture_radius(const FixtureContext *ctx, double t) { + if (ctx->model == 0) + return ctx->R0 + ctx->rr * (t - ctx->t0); + return ctx->R0 + ctx->rr * t; +} + +static SpacetimePointStatus fixture_eval(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { + (void)source; + (void)t; + (void)x; + *metric = (MetricData){.alpha = 1.0, + .gamma = {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}}; + return SPACETIME_POINT_OK; +} + +static SpacetimeRayStatus fixture_classify(const SpacetimeSource *source, + double t, const double x[3]) { + (void)source; + (void)t; + (void)x; + return SPACETIME_RAY_ACTIVE; +} + +static size_t fixture_end_count(const SpacetimeSource *source) { + (void)source; + return 1; +} + +static int fixture_end(const SpacetimeSource *source, size_t index, + SpacetimeAsymptoticEnd *out) { + const FixtureContext *ctx = source->context; + if (index != 0) + return -1; + *out = (SpacetimeAsymptoticEnd){ + .end_id = 0, + .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, + .mass = 0.0, + .frame_origin = {0.0, 0.0, 0.0}, + .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + return 0; +} + +static int fixture_worldtube(const SpacetimeSource *source, + SpacetimeEndId end_id, double t, + SpacetimeEscapeWorldtubeSample *out) { + const FixtureContext *ctx = source->context; + if (end_id != 0) + return -1; + if (!isfinite(t) || t < ctx->valid_t_min) { + *out = (SpacetimeEscapeWorldtubeSample){.valid = 0}; + return 0; + } + *out = (SpacetimeEscapeWorldtubeSample){ + .center = {fixture_center(ctx, 0, t), fixture_center(ctx, 1, t), + fixture_center(ctx, 2, t)}, + .velocity = {ctx->v[0], ctx->v[1], ctx->v[2]}, + .radius = fixture_radius(ctx, t), + .radius_rate = ctx->rr, + .velocity_constant = 1, + .valid = 1}; + return 0; +} + +static double fixture_next_segment(const SpacetimeSource *source, + SpacetimeEndId end_id, double t) { + (void)source; + (void)end_id; + (void)t; + return NAN; +} + +static void fixture_destroy(SpacetimeSource *source) { + source->context = NULL; + source->ops = NULL; +} + +static const SpacetimeOps fixture_ops = { + .eval = fixture_eval, + .classify = fixture_classify, + .asymptotic_end_count = fixture_end_count, + .asymptotic_end = fixture_end, + .escape_worldtube_sample = fixture_worldtube, + .escape_worldtube_next_segment = fixture_next_segment, + .destroy = fixture_destroy, +}; + +/* ------------------------------------------------------------------ */ +static void print_environment(void) { + fprintf(stdout, + "{\"variant\":\"%s\",\"sizeof_long_double\":%zu,\"LDBL_MANT_DIG\":%d," + "\"DBL_MANT_DIG\":%d,\"LDBL_MAX_EXP\":%d,\"hardware_threads\":%d," + "\"omp_max_threads\":%d}\n", + PROBE_VARIANT_NAME, sizeof(long double), LDBL_MANT_DIG, DBL_MANT_DIG, + LDBL_MAX_EXP, omp_get_num_procs(), omp_get_max_threads()); +} + +/* ------------------------------------------------------------------ */ +/* accuracy mode */ +static void mode_accuracy(const char *kernel_out, const char *route_out) { + FILE *kf = fopen(kernel_out, "w"); + if (kf == NULL) { + fprintf(stderr, "cannot open %s\n", kernel_out); + exit(2); + } + fprintf(kf, + "id,category,status,sigma,x0,x1,x2,c0,c1,c2,w0,w1,w2,v0,v1,v2,R0,rr," + "a,b,c\n"); + for (int i = 0; i < quad_kernel_case_count; ++i) { + const QuadKernelCase *c = &quad_kernel_cases[i]; + EntryQuadratic k = entry_quadratic_coeffs(c->x, c->c, c->w, c->v, c->R0, + c->rr); + double s = -1.0; + EntrySolveResult r = entry_solve(&k, &s); + fprintf(kf, "%d,%s,%d,%a", c->id, c->category, (int)r, s); + for (int j = 0; j < 3; ++j) + fprintf(kf, ",%a", c->x[j]); + for (int j = 0; j < 3; ++j) + fprintf(kf, ",%a", c->c[j]); + for (int j = 0; j < 3; ++j) + fprintf(kf, ",%a", c->w[j]); + for (int j = 0; j < 3; ++j) + fprintf(kf, ",%a", c->v[j]); + fprintf(kf, ",%a,%a,", c->R0, c->rr); + PRINT_COEFF(kf, k.a); + fputc(',', kf); + PRINT_COEFF(kf, k.b); + fputc(',', kf); + PRINT_COEFF(kf, k.c); + fputc('\n', kf); + } + fclose(kf); + + FILE *rf = fopen(route_out, "w"); + if (rf == NULL) { + fprintf(stderr, "cannot open %s\n", route_out); + exit(2); + } + fprintf(rf, + "id,category,status,kind,failure_reason,fallback_evals,pi_match," + "lcam_match,F,tol,why,x0,x1,x2,Pi0,Pi1,Pi2,ninf0,ninf1,ninf2," + "canon_t,canonx0,canonx1,canonx2,canonw0,canonw1,canonw2," + "cbx0,cbx1,cbx2,cbr,cbok,ct00,ct01,ct02,rt0,cb0ok," + "kern_ok,kern_status,kern_sigma,kern_a,kern_b,kern_c," + "logcamera,logentry," + "t0,oc0,oc1,oc2,d0,d1,d2,fc0,fc1,fc2,v0,v1,v2," + "R0,rr,model,reason_name\n"); + for (int i = 0; i < quad_route_case_count; ++i) { + const QuadRouteCase *c = &quad_route_cases[i]; + FixtureContext ctx = {.c0 = {c->c0[0], c->c0[1], c->c0[2]}, + .v = {c->v[0], c->v[1], c->v[2]}, + .R0 = c->R0, + .rr = c->rr, + .t0 = c->t0, + .valid_t_min = c->valid_t_min, + .model = c->model}; + SpacetimeSource source = {.ops = &fixture_ops, .context = &ctx}; + ObserverState obs = {0}; + obs.coordinate_time = c->t0; + for (int j = 0; j < 3; ++j) + obs.coordinate_position[j] = c->obs[j]; + obs.tetrad[0][0] = 1.0; + obs.tetrad[1][1] = 1.0; + obs.tetrad[2][2] = 1.0; + obs.tetrad[3][3] = 1.0; + + int canon_ok = 0; + AsymptoticPhotonState canon = {0}; + GeodesicRayState st = {0}; + { + MetricData metric; + if (spacetime_eval(&source, c->t0, obs.coordinate_position, &metric) == + SPACETIME_POINT_OK && + geodesic_initialize_past_ray_metric(&metric, &obs, c->dir, &st) == 0 && + asymptotic_canonical_from_backend(&source, 0, &metric, c->t0, st.x, + st.Pi, st.log_alpha_p0, + &canon) == 0) + canon_ok = 1; + } + + AsymptoticRoute route; + AsymptoticStatus status = + asymptotic_route_camera(&source, &obs, c->dir, &route); + + int pi_match = 1, lcam_match = 1; + if (canon_ok && status == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY) { + for (int j = 0; j < 3; ++j) + if (route.Pi[j] != st.Pi[j]) + pi_match = 0; + if (route.log_alpha_p0_camera != st.log_alpha_p0) + lcam_match = 0; + } + + double F = NAN, tol = NAN; + RayReason why = RAY_REASON_NONE; + if (status == ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY) + asymptotic_entry_geometry(&source, route.end_id, route.activate_t, + route.x, &F, &tol, &why); + + SpacetimeEscapeWorldtubeSample cb = {0}; + int cbok = 0; + { + const double tt = (status == ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY) + ? route.activate_t + : c->t0; + if (spacetime_escape_worldtube_sample(&source, route.end_id, tt, &cb) == + 0 && + cb.valid) + cbok = 1; + } + /* Worldtube sample at the segment start used by the quadratic kernel. */ + SpacetimeEscapeWorldtubeSample cb0 = {0}; + int cb0ok = 0; + if (spacetime_escape_worldtube_sample(&source, route.end_id, c->t0, + &cb0) == 0 && + cb0.valid) + cb0ok = 1; + + /* Reproduce the first-segment public kernel classification at the actual + * canonical inputs, so a reference ENTER / kernel MISS / route ESCAPED is + * directly visible and not confused with canonical normalisation. */ + int kern_ok = 0; + int kern_status = -1; + double kern_sigma = -1.0; + EntryQuadratic kk = {0}; + SpacetimeAsymptoticEnd end_desc; + if (canon_ok && cb0ok && + spacetime_asymptotic_end(&source, 0, &end_desc) == 0) { + double c_frame[3], v_frame[3]; + backend_position_to_frame(&end_desc, cb0.center, c_frame); + backend_vector_to_frame(&end_desc, cb0.velocity, v_frame); + kk = entry_quadratic_coeffs(canon.x, c_frame, canon.w, v_frame, + cb0.radius, cb0.radius_rate); + kern_status = (int)entry_solve(&kk, &kern_sigma); + kern_ok = 1; + } + + fprintf(rf, "%d,%s,%d,%d,%d,%u,%d,%d,%a,%a,%d", c->id, c->category, + (int)status, (int)route.kind, (int)route.failure_reason, + route.entry_fallback_evaluations, pi_match, lcam_match, F, tol, + (int)why); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", route.x[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", route.Pi[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", route.n_infinity[j]); + fprintf(rf, ",%a", canon.t); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", canon.x[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", canon.w[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", cb.center[j]); + fprintf(rf, ",%a,%d", cb.radius, cbok); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", cb0.center[j]); + fprintf(rf, ",%a,%d", cb0.radius, cb0ok); + fprintf(rf, ",%d,%d,%a,", kern_ok, kern_status, kern_sigma); + PRINT_COEFF(rf, kk.a); + fputc(',', rf); + PRINT_COEFF(rf, kk.b); + fputc(',', rf); + PRINT_COEFF(rf, kk.c); + fprintf(rf, ",%a,%a", route.log_alpha_p0_camera, route.log_alpha_p0); + fprintf(rf, ",%a", c->t0); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", c->obs[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", c->dir[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", c->c0[j]); + for (int j = 0; j < 3; ++j) + fprintf(rf, ",%a", c->v[j]); + fprintf(rf, ",%a,%a,%d,%s\n", c->R0, c->rr, c->model, + ray_reason_name(route.failure_reason)); + } + fclose(rf); +} + +/* ------------------------------------------------------------------ */ +/* microbench mode: coefficient assembly + entry_solve, serial. */ +static volatile double g_kernel_sink; + +static __attribute__((noinline)) double kernel_batch(const QuadKernelCase *cs, + int n, long reps) { + double acc = 0.0; + for (long r = 0; r < reps; ++r) { + for (int i = 0; i < n; ++i) { + /* Force the index through an opaque register so the compiler cannot + * hoist the pure coefficient solve out of the repetition loop or prove + * the loaded fixture invariant. Identical for every variant. */ + int idx = i; + __asm__ __volatile__("" : "+r"(idx) : : "memory"); + const QuadKernelCase *c = cs + idx; + EntryQuadratic k = entry_quadratic_coeffs(c->x, c->c, c->w, c->v, + c->R0, c->rr); + double s = 0.0; + acc += (double)entry_solve(&k, &s) + s * 1e-300; + } + } + return acc; +} + +static void mode_microbench(long target_calls, const char *out) { + const int n = quad_kernel_case_count; + long reps = target_calls / n; + if (reps < 1) + reps = 1; + /* Warm-up outside the timed region. */ + g_kernel_sink += kernel_batch(quad_kernel_cases, n, 1); + const double t0 = omp_get_wtime(); + const double acc = kernel_batch(quad_kernel_cases, n, reps); + const double t1 = omp_get_wtime(); + g_kernel_sink += acc; + const long calls = (long)n * reps; + const double seconds = t1 - t0; + FILE *f = fopen(out, "w"); + if (f == NULL) { + fprintf(stderr, "cannot open %s\n", out); + exit(2); + } + fprintf(f, + "{\"variant\":\"%s\",\"mode\":\"microbench\",\"cases\":%d," + "\"reps\":%ld,\"calls\":%ld,\"seconds\":%.9f,\"ns_per_call\":%.6f," + "\"sink\":%.17g}\n", + PROBE_VARIANT_NAME, n, reps, calls, seconds, + seconds * 1e9 / (double)calls, g_kernel_sink); + fclose(f); + fprintf(stdout, "microbench %s: %ld calls in %.6f s (%.2f ns/call)\n", + PROBE_VARIANT_NAME, calls, seconds, seconds * 1e9 / (double)calls); +} + +/* ------------------------------------------------------------------ */ +/* routebench mode: public asymptotic_route_camera, OpenMP static. */ +typedef struct { + FixtureContext *ctxs; + SpacetimeSource *srcs; + ObserverState *obss; + const QuadRouteCase *rcs; + int n; +} Preloaded; + +typedef struct { + long entry, escaped, inside, time_exhausted, invalid, unsupported, other; + long fallback_sum; +} RouteCounts; + +static volatile double g_route_sink; + +static __attribute__((noinline)) void route_batch(const Preloaded *p, long reps, + int threads, double *seconds, + RouteCounts *counts) { + const int n = p->n; + long entry = 0, escaped = 0, inside = 0, texh = 0, inv = 0, unsup = 0, + other = 0, fallback = 0; + const long total = (long)n * reps; + const double t0 = omp_get_wtime(); +#pragma omp parallel num_threads(threads) reduction(+ : entry, escaped, inside, texh, inv, unsup, other, fallback) + { +#pragma omp for schedule(static) + for (long k = 0; k < total; ++k) { + long kk = k; + __asm__ __volatile__("" : "+r"(kk) : : "memory"); + const int i = (int)(kk % n); + AsymptoticRoute route; + const AsymptoticStatus st = asymptotic_route_camera( + &p->srcs[i], &p->obss[i], p->rcs[i].dir, &route); + fallback += (long)route.entry_fallback_evaluations; + if (st == ASYMPTOTIC_OK) { + switch (route.kind) { + case ASYMPTOTIC_ROUTE_ENTRY: + ++entry; + break; + case ASYMPTOTIC_ROUTE_ESCAPED: + ++escaped; + break; + case ASYMPTOTIC_ROUTE_INSIDE: + ++inside; + break; + case ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED: + ++texh; + break; + default: + ++other; + break; + } + } else if (st == ASYMPTOTIC_UNSUPPORTED) { + ++unsup; + } else { + ++inv; + } + } + } + *seconds = omp_get_wtime() - t0; + counts->entry = entry; + counts->escaped = escaped; + counts->inside = inside; + counts->time_exhausted = texh; + counts->invalid = inv; + counts->unsupported = unsup; + counts->other = other; + counts->fallback_sum = fallback; +} + +static void mode_routebench(long target_calls, int threads, double max_seconds, + const char *out) { + const int n = quad_route_case_count; + Preloaded p; + p.n = n; + p.rcs = quad_route_cases; + p.ctxs = malloc(sizeof(FixtureContext) * (size_t)n); + p.srcs = malloc(sizeof(SpacetimeSource) * (size_t)n); + p.obss = malloc(sizeof(ObserverState) * (size_t)n); + if (!p.ctxs || !p.srcs || !p.obss) { + fprintf(stderr, "allocation failure\n"); + exit(2); + } + for (int i = 0; i < n; ++i) { + const QuadRouteCase *c = &quad_route_cases[i]; + p.ctxs[i] = (FixtureContext){.c0 = {c->c0[0], c->c0[1], c->c0[2]}, + .v = {c->v[0], c->v[1], c->v[2]}, + .R0 = c->R0, + .rr = c->rr, + .t0 = c->t0, + .valid_t_min = c->valid_t_min, + .model = c->model}; + p.srcs[i] = (SpacetimeSource){.ops = &fixture_ops, .context = &p.ctxs[i]}; + p.obss[i] = (ObserverState){0}; + p.obss[i].coordinate_time = c->t0; + for (int j = 0; j < 3; ++j) + p.obss[i].coordinate_position[j] = c->obs[j]; + p.obss[i].tetrad[0][0] = 1.0; + p.obss[i].tetrad[1][1] = 1.0; + p.obss[i].tetrad[2][2] = 1.0; + p.obss[i].tetrad[3][3] = 1.0; + } + + /* Calibrate with one pass, then size reps to the call/time budgets. */ + double cal_seconds = 0.0; + RouteCounts cal_counts; + route_batch(&p, 1, threads, &cal_seconds, &cal_counts); + const double per_call = cal_seconds / (double)n; + long reps = target_calls / n; + if (reps < 1) + reps = 1; + if (per_call > 0.0) { + const long by_time = (long)(max_seconds / (per_call * (double)n)); + if (by_time < 1) + reps = 1; + else if (reps > by_time) + reps = by_time; + } + + double seconds = 0.0; + RouteCounts counts; + route_batch(&p, reps, threads, &seconds, &counts); + g_route_sink += (double)counts.entry + (double)counts.escaped; + + FILE *f = fopen(out, "w"); + if (f == NULL) { + fprintf(stderr, "cannot open %s\n", out); + exit(2); + } + fprintf(f, + "{\"variant\":\"%s\",\"mode\":\"routebench\",\"cases\":%d," + "\"reps\":%ld,\"calls\":%ld,\"threads\":%d,\"cal_seconds\":%.9f," + "\"seconds\":%.9f,\"ns_per_call\":%.6f,\"entry\":%ld,\"escaped\":%ld," + "\"inside\":%ld,\"time_exhausted\":%ld,\"invalid\":%ld," + "\"unsupported\":%ld,\"other\":%ld,\"fallback_sum\":%ld}\n", + PROBE_VARIANT_NAME, n, reps, (long)n * reps, threads, cal_seconds, + seconds, seconds * 1e9 / (double)((long)n * reps), counts.entry, + counts.escaped, counts.inside, counts.time_exhausted, counts.invalid, + counts.unsupported, counts.other, counts.fallback_sum); + fclose(f); + fprintf(stdout, + "routebench %s: %ld calls in %.6f s on %d threads (%.2f ns/call)\n", + PROBE_VARIANT_NAME, (long)n * reps, seconds, threads, + seconds * 1e9 / (double)((long)n * reps)); + free(p.ctxs); + free(p.srcs); + free(p.obss); +} + +/* ------------------------------------------------------------------ */ +static void usage(const char *argv0) { + fprintf(stderr, + "usage:\n" + " %s accuracy --kernel-out K.csv --route-out R.csv\n" + " %s microbench --out F.json [--target-calls N]\n" + " %s routebench --out F.json [--target-calls N] [--threads T]" + " [--max-seconds S]\n", + argv0, argv0, argv0); +} + +static const char *arg_value(int argc, char **argv, const char *flag) { + for (int i = 1; i + 1 < argc; ++i) + if (strcmp(argv[i], flag) == 0) + return argv[i + 1]; + return NULL; +} + +int main(int argc, char **argv) { + print_environment(); + if (argc < 2) { + usage(argv[0]); + return 1; + } + if (strcmp(argv[1], "accuracy") == 0) { + const char *k = arg_value(argc, argv, "--kernel-out"); + const char *r = arg_value(argc, argv, "--route-out"); + if (!k || !r) { + usage(argv[0]); + return 1; + } + mode_accuracy(k, r); + return 0; + } + if (strcmp(argv[1], "microbench") == 0) { + const char *out = arg_value(argc, argv, "--out"); + const char *tc = arg_value(argc, argv, "--target-calls"); + if (!out) { + usage(argv[0]); + return 1; + } + mode_microbench(tc ? atol(tc) : 10000000L, out); + return 0; + } + if (strcmp(argv[1], "routebench") == 0) { + const char *out = arg_value(argc, argv, "--out"); + const char *tc = arg_value(argc, argv, "--target-calls"); + const char *th = arg_value(argc, argv, "--threads"); + const char *ms = arg_value(argc, argv, "--max-seconds"); + if (!out) { + usage(argv[0]); + return 1; + } + mode_routebench(tc ? atol(tc) : 1000000L, th ? atoi(th) : 4, + ms ? atof(ms) : 20.0, out); + return 0; + } + usage(argv[0]); + return 1; +} diff --git a/benchmarks/quadratic_precision/oracle.py b/benchmarks/quadratic_precision/oracle.py new file mode 100644 index 0000000..7c91124 --- /dev/null +++ b/benchmarks/quadratic_precision/oracle.py @@ -0,0 +1,571 @@ +#!/usr/bin/env python3 +"""High-precision reference for the quadratic entry benchmark. + +Inputs are exact IEEE-754 values, so they are represented exactly as +``fractions.Fraction``. The relative-distance quadratic coefficients, the +discriminant and the polynomial residuals are therefore *exact* rationals; only +the square root needs ``decimal`` (160 significant digits by default). This +avoids ``float.fromhex`` (which would silently drop a long-double coefficient to +53 bits) and keeps the sign of a near-zero discriminant exact. + +Reference classification mirrors the production contract: + * c < 0 -> camera already INSIDE (never an entry failure) + * c == 0 -> boundary slope b decides; b < 0 enters at once, b == 0 with a < 0 + enters later, otherwise no crossing + * c > 0 -> smallest positive root with inward slope; a missing/tangent root + is MISS. ``a == 0`` is the exact linear branch. +""" +from __future__ import annotations + +import csv +import json +import math +from decimal import Decimal, localcontext +from fractions import Fraction + +DEFAULT_PREC = 160 + + +def parse_hex(s: str) -> Fraction: + """Exact rational value of a C ``%a`` / ``%La`` hex float literal.""" + s = s.strip() + low = s.lower() + if "nan" in low: + raise ValueError("NaN input") + if "inf" in low: + return Fraction(0) # only used for flags; never present in fixtures + neg = False + if s and s[0] in "+-": + neg = s[0] == "-" + s = s[1:] + if s[:2].lower() == "0x": + s = s[2:] + mant, _, exp_s = s.partition("p") + if not exp_s: + mant, _, exp_s = s.partition("P") + exp = int(exp_s) if exp_s else 0 + ip, _, fp = mant.partition(".") + digits = (ip + fp) or "0" + val = Fraction(int(digits, 16), 1) + shift = exp - 4 * len(fp) + if shift >= 0: + val *= Fraction(2) ** shift + else: + val /= Fraction(2) ** (-shift) + return -val if neg else val + + +def frac_to_dec(fr: Fraction, prec: int = DEFAULT_PREC) -> Decimal: + with localcontext() as ctx: + ctx.prec = prec + return Decimal(fr.numerator) / Decimal(fr.denominator) + + +def hex_to_dec(s: str, prec: int = DEFAULT_PREC) -> Decimal: + return frac_to_dec(parse_hex(s), prec) + + +def classify(a: Fraction, b: Fraction, c: Fraction, prec: int, + R0: Fraction | None = None, rr: Fraction | None = None): + """Return (status, root_decimal_or_None, info). + + ``R0``/``rr`` (radius at the segment start and dR/dt) enforce the physical + positive-radius domain R(s) = R0 - rr*s > 0 along the past parameter. A + positive root outside that domain is not a physical entry: it is reported + as MISS with ``info['domain_clipped']`` set (shrinking worldtubes). + """ + info = {"a_zero": a == 0, "disc_sign": 0, "domain_clipped": False} + + def domain_ok(root: Fraction) -> bool: + if R0 is None: + return True + return (R0 - (rr if rr is not None else Fraction(0)) * root) > 0 + + if c < 0: + return "INSIDE", None, info + if c == 0: + if b < 0: + return "ENTER", Decimal(0), info + if b == 0: + if a < 0: + return "ENTER", Decimal(0), info + return "MISS", None, info + if a < 0: + root = -b / a + if root > 0 and domain_ok(root): + return "ENTER", frac_to_dec(root, prec), info + if root > 0: + info["domain_clipped"] = True + return "MISS", None, info + return "MISS", None, info + if a == 0: + if b < 0: + root = -c / b + if root > 0: + if domain_ok(root): + return "ENTER", frac_to_dec(root, prec), info + info["domain_clipped"] = True + return "MISS", None, info + disc = b * b - 4 * a * c + info["disc_sign"] = (disc > 0) - (disc < 0) + if disc < 0: + return "MISS", None, info + if disc == 0: + return "MISS", None, info # tangency is not a crossing + with localcontext() as ctx: + ctx.prec = prec + sd = frac_to_dec(disc, prec).sqrt() + ad = frac_to_dec(a, prec) + bd = frac_to_dec(b, prec) + r1 = (-bd - sd) / (2 * ad) + r2 = (-bd + sd) / (2 * ad) + pos = sorted(r for r in (r1, r2) if r > 0) + for r in pos: + if 2 * ad * r + bd < 0: # inward (outside -> inside) slope + rfr = _dec_to_frac_snapshot(r) + if domain_ok(rfr): + return "ENTER", r, info + info["domain_clipped"] = True + return "MISS", None, info + + +def _dec_to_frac_snapshot(d: Decimal) -> Fraction: + return Fraction(d) + + +def coeffs(x, c, w, v, R0, rr): + d = [x[i] - c[i] for i in range(3)] + q = [w[i] + v[i] for i in range(3)] + qq = sum(q[i] * q[i] for i in range(3)) + dq = sum(d[i] * q[i] for i in range(3)) + dd = sum(d[i] * d[i] for i in range(3)) + a = qq - rr * rr + b = 2 * (dq + R0 * rr) + cq = dd - R0 * R0 + return a, b, cq, qq + + +def load_cases(path): + data = json.loads(open(path).read()) + kernel = {} + for c in data["kernel"]: + kernel[c["id"]] = { + "category": c["category"], + "x": [parse_hex(z) for z in c["x"]], + "c": [parse_hex(z) for z in c["c"]], + "w": [parse_hex(z) for z in c["w"]], + "v": [parse_hex(z) for z in c["v"]], + "R0": parse_hex(c["R0"]), + "rr": parse_hex(c["rr"]), + } + route = {c["id"]: c for c in data["route"]} + return kernel, route + + +def kernel_reference(case, prec): + a, b, cq, qq = coeffs(case["x"], case["c"], case["w"], case["v"], + case["R0"], case["rr"]) + status, root, info = classify(a, b, cq, prec, case["R0"], case["rr"]) + dd = sum(xi * xi for xi in + (case["x"][i] - case["c"][i] for i in range(3))) + return { + "status": status, + "root": root, + "a": a, + "b": b, + "c": cq, + "qq": qq, + "dd": dd, + "R0": case["R0"], + "rr": case["rr"], + "a_zero": info["a_zero"], + "disc_sign": info["disc_sign"], + "domain_clipped": info["domain_clipped"], + } + + +def ulp_dec(cand: float) -> Decimal: + return frac_to_dec(Fraction(math.ulp(cand))) + + +def dec_of_float(f: float) -> Decimal: + return frac_to_dec(Fraction(f), 300) + + +def poly_resid(a: Fraction, b: Fraction, c: Fraction, s: Fraction) -> Fraction: + return a * s * s + b * s + c + + +def poly_scale(a: Fraction, b: Fraction, c: Fraction, s: Fraction) -> Fraction: + return abs(a) * s * s + abs(b) * s + abs(c) + + +def analyze_kernel_variant(rows, refs, prec): + """Return per-case records for one variant's kernel CSV.""" + out = [] + for r in rows: + cid = int(r["id"]) + ref = refs[cid] + status = int(r["status"]) + cand = float.fromhex(r["sigma"]) if r["sigma"] not in ("", "nan") else float("nan") + cand_fr = parse_hex(r["sigma"]) if r["sigma"] not in ("", "nan") else None + ak = parse_hex(r["a"]) + bk = parse_hex(r["b"]) + ck = parse_hex(r["c"]) + rec = { + "id": cid, + "category": r["category"], + "variant_status": status, + "ref_status": ref["status"], + "a_zero_ref": ref["a_zero"], + "a_zero_kernel": ak == 0, + "a_sign_ref": _sign(ref["a"]), + "a_sign_kernel": _sign(ak), + "disc_sign_ref": ref["disc_sign"], + "late": False, + "near_linear": False, + "near_boundary": False, + "ill_conditioned": False, + "domain_clipped": ref["domain_clipped"], + "root_err_ulps": None, + "root_rel_err": None, + "resid_kernel_rel": None, + "resid_ideal_rel": None, + "baseline_resid_rel": None, + } + # conditioning flags + if ref["status"] == "ENTER" and ref["root"] is not None: + s = ref["root"] + scale = max(abs(ref["R0"]), Fraction(1)) + rec["late"] = bool(s > frac_to_dec(scale, prec) * (10 ** 6)) + rec["near_linear"] = abs(ref["a"]) <= Fraction(1, 10**10) * max( + ref["qq"], ref["rr"] * ref["rr"], Fraction(1) + ) + r0sq = ref["R0"] * ref["R0"] + dd = ref["dd"] + rec["near_boundary"] = bool( + dd + r0sq != 0 + and abs(ref["c"]) <= Fraction(1, 10**8) * (dd + r0sq) + ) + rec["ill_conditioned"] = bool( + rec["near_linear"] or rec["late"] or rec["near_boundary"] + or r["category"].startswith("growing_out") + or r["category"] == "shrinking" + ) + # root-level comparison + if status == 1 and ref["status"] == "ENTER" and ref["root"] is not None and cand_fr is not None: + err = abs(frac_to_dec(cand_fr, prec) - ref["root"]) + u = ulp_dec(cand) + if u > 0: + rec["root_err_ulps"] = float(err / u) + if ref["root"] != 0: + rec["root_rel_err"] = float(err / abs(ref["root"])) + rec["resid_kernel_rel"] = _rel( + poly_resid(ak, bk, ck, cand_fr), ak, bk, ck, cand_fr + ) + rec["resid_ideal_rel"] = _rel( + poly_resid(ref["a"], ref["b"], ref["c"], cand_fr), + ref["a"], ref["b"], ref["c"], cand_fr, + ) + try: + nearest = float(ref["root"]) + nfr = Fraction(nearest) + rec["baseline_resid_rel"] = _rel( + poly_resid(ref["a"], ref["b"], ref["c"], nfr), + ref["a"], ref["b"], ref["c"], nfr, + ) + except (OverflowError, ValueError): + rec["baseline_resid_rel"] = None + out.append(rec) + return out + + +def _sign(fr: Fraction) -> int: + return (fr > 0) - (fr < 0) + + +def _rel(resid: Fraction, a, b, c, s) -> float: + scale = poly_scale(a, b, c, s) + if scale == 0: + return 0.0 + return float(abs(resid) / scale) + + +def route_reference(row, prec): + x = [parse_hex(row[f"canonx{i}"]) for i in range(3)] + w = [parse_hex(row[f"canonw{i}"]) for i in range(3)] + v = [parse_hex(row[f"v{i}"]) for i in range(3)] + center = [parse_hex(row[f"ct0{i}"]) for i in range(3)] + rr = parse_hex(row["rr"]) + R0 = parse_hex(row["rt0"]) if row["rt0"] not in ("", "nan") else parse_hex(row["R0"]) + a, b, cq, qq = coeffs(x, center, w, v, R0, rr) + status, root, info = classify(a, b, cq, prec, R0, rr) + return status, root, a, b, cq, qq, info["domain_clipped"] + + +STATUS_NAME = {-1: "INVALID", 0: "MISS", 1: "ENTRY", 2: "UNCERTAIN"} +KIND_NAME = { + 0: "INSIDE", + 1: "ENTRY", + 2: "ESCAPED", + 3: "TIME_RANGE_EXHAUSTED", + 4: "INVALID", +} + + +def analyze_route_variant(rows, prec, dbl_eps=2.220446049250313e-16): + out = [] + for r in rows: + (ref_status, ref_root, a, b, cq, qq, + ref_domain_clipped) = route_reference(r, prec) + status = int(r["status"]) + kind = int(r["kind"]) + F = _maybe_dec(r["F"]) + tol = _maybe_dec(r["tol"]) + kern_ok = r.get("kern_ok", "0") == "1" + kern_status = int(r["kern_status"]) if kern_ok else None + rec = { + "id": int(r["id"]), + "category": r["category"], + "ref_status": ref_status, + "ref_domain_clipped": ref_domain_clipped, + "status": status, + "kind": kind, + "kind_name": KIND_NAME.get(kind, "?"), + "failure_reason": int(r["failure_reason"]), + "fallback": int(r["fallback_evals"]), + "pi_match": int(r["pi_match"]), + "lcam_match": int(r["lcam_match"]), + "kern_status": kern_status if kern_ok else "", + "kern_sigma": r.get("kern_sigma", ""), + "F_over_tol": None, + "inside_mismatch": False, + "false_miss": False, + "false_candidate": False, + "unconfirmed": kind == 4 and r["reason_name"] == "ENTRY_UNCONFIRMED", + "kernel_false_miss": False, + "kernel_false_miss_escaped": False, + } + if ref_status == "INSIDE": + if kind != 0: + rec["inside_mismatch"] = True + elif ref_status == "ENTER": + if kind == 2: + rec["false_miss"] = True + elif ref_status == "MISS": + if kind == 1: + rec["false_candidate"] = True + if ref_status == "ENTER" and kern_ok and kern_status == 0: + rec["kernel_false_miss"] = True + if kind == 2: + rec["kernel_false_miss_escaped"] = True + if F is not None and tol is not None and tol > 0: + rec["F_over_tol"] = float(F / tol) + out.append(rec) + return out + + +def _maybe_dec(s): + if s in ("", "nan"): + return None + low = s.lower() + if "inf" in low: + return None + return hex_to_dec(s, 300) + + +def summarize_kernel(records_by_variant, refs): + """Aggregate counts, ULP distributions, and common-ENTRY intersections.""" + variants = list(records_by_variant.keys()) + by_id = {v: {} for v in variants} + for v in variants: + for rec in records_by_variant[v]: + by_id[v][rec["id"]] = rec + + ids = sorted(refs.keys()) + summary = {"variants": {}} + for v in variants: + recs = by_id[v] + counts = { + "enter": 0, "miss": 0, "uncertain": 0, "false_miss": 0, + "false_uncertain": 0, "false_candidate": 0, + "a_zero_ref": 0, "a_zero_kernel": 0, "a_sign_mismatch": 0, + "a_zero_collapse": 0, "a_spurious_nonzero": 0, + "domain_clipped_candidate": 0, + "ill_conditioned_enter": 0, + "near_boundary_enter": 0, + "domain_clipped_cases": 0, + "late_ref_enter": 0, + } + ulps = [] + rels = [] + for cid in ids: + rec = recs[cid] + st = rec["variant_status"] + if st == 1: + counts["enter"] += 1 + elif st == 0: + counts["miss"] += 1 + elif st == 2: + counts["uncertain"] += 1 + if rec["ref_status"] == "ENTER" and st == 0: + counts["false_miss"] += 1 + if rec["ref_status"] == "ENTER" and st == 2: + counts["false_uncertain"] += 1 + if rec["ref_status"] == "MISS" and st == 1: + if rec["domain_clipped"]: + counts["domain_clipped_candidate"] += 1 + else: + counts["false_candidate"] += 1 + if rec["a_zero_ref"]: + counts["a_zero_ref"] += 1 + if rec["a_zero_kernel"]: + counts["a_zero_kernel"] += 1 + if not rec["a_zero_ref"] and rec["a_zero_kernel"]: + counts["a_zero_collapse"] += 1 + if rec["a_zero_ref"] and not rec["a_zero_kernel"]: + counts["a_spurious_nonzero"] += 1 + if rec["ill_conditioned"] and rec["ref_status"] == "ENTER": + counts["ill_conditioned_enter"] += 1 + if rec["near_boundary"] and rec["ref_status"] == "ENTER": + counts["near_boundary_enter"] += 1 + if rec["domain_clipped"]: + counts["domain_clipped_cases"] += 1 + if (rec["a_sign_ref"] != 0 and rec["a_sign_kernel"] != 0 + and rec["a_sign_ref"] != rec["a_sign_kernel"]): + counts["a_sign_mismatch"] += 1 + if rec["late"]: + counts["late_ref_enter"] += 1 + if rec["root_err_ulps"] is not None: + ulps.append(rec["root_err_ulps"]) + if rec["root_rel_err"] is not None: + rels.append(rec["root_rel_err"]) + summary["variants"][v] = { + "counts": counts, + "root_ulp_all_ref_enter": _dist(ulps), + "root_rel_all_ref_enter": _dist(rels), + } + + # Common reference-ENTER subset that every variant solved as ENTRY. + common = [] + for cid in ids: + if refs[cid]["status"] != "ENTER": + continue + if all(by_id[v][cid]["variant_status"] == 1 for v in variants): + common.append(cid) + summary["common_enter_ids"] = len(common) + for v in variants: + vals = [] + vals_nonlate = [] + rels_nonlate = [] + for cid in common: + rec = by_id[v][cid] + e = rec["root_err_ulps"] + if e is not None: + vals.append(e) + if not rec["ill_conditioned"]: + vals_nonlate.append(e) + if rec["root_rel_err"] is not None: + rels_nonlate.append(rec["root_rel_err"]) + summary["variants"][v]["root_ulp_common_enter"] = _dist(vals) + summary["variants"][v]["root_ulp_common_enter_wellcond"] = _dist(vals_nonlate) + summary["variants"][v]["root_rel_common_enter_wellcond"] = _dist(rels_nonlate) + return summary + + +def _dist(vals): + if not vals: + return {"n": 0} + vs = sorted(vals) + n = len(vs) + def pct(p): + idx = min(n - 1, max(0, int(math.ceil(p * n)) - 1)) + return vs[idx] + return { + "n": n, + "min": vs[0], + "median": pct(0.5), + "p95": pct(0.95), + "max": vs[-1], + } + + +def summarize_route(records_by_variant): + summary = {"variants": {}} + for v, recs in records_by_variant.items(): + counts = { + "inside_ok": 0, "inside_mismatch": 0, "false_miss": 0, + "false_candidate": 0, "unconfirmed": 0, "fallback_used": 0, + "accepted_entry": 0, "F_over_tol_gt1": 0, "pi_mismatch": 0, + "lcam_mismatch": 0, "kernel_false_miss": 0, + "kernel_false_miss_escaped": 0, "ref_domain_clipped": 0, + } + fmax = 0.0 + by_cat = {} + for rec in recs: + cat = by_cat.setdefault(rec["category"], {"cases": 0, "fallback": 0, + "entry": 0, "escaped": 0}) + cat["cases"] += 1 + if rec["kind"] == 0 and rec["ref_status"] == "INSIDE": + counts["inside_ok"] += 1 + if rec["inside_mismatch"]: + counts["inside_mismatch"] += 1 + if rec["false_miss"]: + counts["false_miss"] += 1 + if rec["false_candidate"]: + counts["false_candidate"] += 1 + if rec["unconfirmed"]: + counts["unconfirmed"] += 1 + if rec["kernel_false_miss"]: + counts["kernel_false_miss"] += 1 + if rec["kernel_false_miss_escaped"]: + counts["kernel_false_miss_escaped"] += 1 + if rec["ref_domain_clipped"]: + counts["ref_domain_clipped"] += 1 + if rec["fallback"] > 0: + counts["fallback_used"] += 1 + cat["fallback"] += 1 + if rec["kind"] == 1: + counts["accepted_entry"] += 1 + cat["entry"] += 1 + if rec["F_over_tol"] is not None: + fmax = max(fmax, rec["F_over_tol"]) + if rec["F_over_tol"] > 1.0: + counts["F_over_tol_gt1"] += 1 + if rec["kind"] == 2: + cat["escaped"] += 1 + if not rec["pi_match"]: + counts["pi_mismatch"] += 1 + if not rec["lcam_match"]: + counts["lcam_mismatch"] += 1 + summary["variants"][v] = {"counts": counts, "max_F_over_tol": fmax, + "by_category": by_cat} + return summary + + +def precision_consistency(cases, ids, prec_a, prec_b): + """Compare reference classification and nearest-double root at two Decimal + precisions. Coefficients/discriminant are exact rationals, so only the + square-root precision can differ. Returns a list of mismatches.""" + mismatches = [] + for cid in sorted(ids): + ra = kernel_reference(cases[cid], prec_a) + rb = kernel_reference(cases[cid], prec_b) + reason = None + if ra["status"] != rb["status"]: + reason = "status" + elif ra["disc_sign"] != rb["disc_sign"]: + reason = "disc_sign" + elif ra["status"] == "ENTER": + try: + fa = float(ra["root"]) + fb = float(rb["root"]) + except (OverflowError, ValueError): + reason = "root_unrepresentable" + else: + if not (fa == fb or (math.isnan(fa) and math.isnan(fb))): + reason = "nearest_double_root" + if reason is not None: + mismatches.append({"id": cid, "reason": reason, + "status_a": ra["status"], "status_b": rb["status"]}) + return mismatches diff --git a/benchmarks/quadratic_precision/run.py b/benchmarks/quadratic_precision/run.py new file mode 100644 index 0000000..e3ecc8e --- /dev/null +++ b/benchmarks/quadratic_precision/run.py @@ -0,0 +1,458 @@ +#!/usr/bin/env python3 +"""End-to-end driver for the quadratic precision/performance benchmark. + +Self-contained: reads the current working-tree ``src/asymptotic.c``, generates +and builds four variants, runs the accuracy probe, evaluates the exact-decimal +reference, and times the kernel and the public pre-route. All raw artifacts +are written under ``--output-dir`` (default /tmp/opencode/quadratic-comparison). + +No production source, Makefile or git state is modified. +""" +from __future__ import annotations + +import argparse +import csv +import json +import os +import re +import subprocess +import sys + +# Do not leave __pycache__ inside the repository benchmark directory: the +# imported local modules (build/cases/oracle) must not be cached here. +sys.dont_write_bytecode = True + +from pathlib import Path + +HERE = Path(__file__).resolve().parent +sys.path.insert(0, str(HERE)) + +import build as buildmod # noqa: E402 +import cases as casesmod # noqa: E402 +import oracle # noqa: E402 + +VARIANTS = list(buildmod.VARIANTS) + + +def run_cmd(cmd, log_path: Path, env=None): + with open(log_path, "a") as fh: + fh.write("$ " + " ".join(str(c) for c in cmd) + "\n") + proc = subprocess.run(cmd, capture_output=True, text=True, check=False, env=env) + with open(log_path, "a") as fh: + if proc.stdout: + fh.write(proc.stdout) + if proc.stderr: + fh.write(proc.stderr) + fh.write(f"exit={proc.returncode}\n") + return proc + + +def read_csv(path: Path): + with open(path, newline="") as fh: + return list(csv.DictReader(fh)) + + +def run_timing(cmd, outp: Path, log_path: Path): + # Never ingest output left by a previous invocation, even after a failure. + outp.unlink(missing_ok=True) + proc = run_cmd(cmd, log_path) + if proc.returncode != 0 or not outp.is_file(): + raise RuntimeError(f"timing failed or produced no fresh output: {log_path}") + return json.loads(outp.read_text()) + + +def write_metrics_csv(path: Path, records, columns): + with open(path, "w", newline="") as fh: + w = csv.DictWriter(fh, fieldnames=columns, extrasaction="ignore") + w.writeheader() + for r in records: + w.writerow(r) + + +def write_curated(path: Path, rows, columns): + with open(path, "w", newline="") as fh: + w = csv.DictWriter(fh, fieldnames=columns, extrasaction="ignore") + w.writeheader() + for r in rows: + w.writerow(r) + + +def build_curated_kernel(kernel_records): + ref_ids = [r["id"] for r in kernel_records[VARIANTS[0]] if r["id"] < 1000] + by = {v: {r["id"]: r for r in kernel_records[v]} for v in VARIANTS} + rows = [] + for cid in sorted(ref_ids): + row = {"id": cid, "category": by[VARIANTS[0]][cid]["category"], + "ref": by[VARIANTS[0]][cid]["ref_status"]} + for v in VARIANTS: + row[f"{v}_status"] = by[v][cid]["variant_status"] + row[f"{v}_ulp"] = by[v][cid]["root_err_ulps"] + row[f"{v}_rel"] = by[v][cid]["root_rel_err"] + rows.append(row) + return rows + + +def build_curated_route(route_records): + ref_ids = [r["id"] for r in route_records[VARIANTS[0]] if r["id"] < 100] + by = {v: {r["id"]: r for r in route_records[v]} for v in VARIANTS} + rows = [] + for cid in sorted(ref_ids): + row = {"id": cid, "category": by[VARIANTS[0]][cid]["category"], + "ref": by[VARIANTS[0]][cid]["ref_status"]} + for v in VARIANTS: + row[f"{v}_kind"] = by[v][cid]["kind"] + row[f"{v}_fallback"] = by[v][cid]["fallback"] + row[f"{v}_F_over_tol"] = by[v][cid]["F_over_tol"] + rows.append(row) + return rows + + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--repo", default=str(HERE.parents[1])) + ap.add_argument("--output-dir", default="/tmp/opencode/quadratic-comparison") + ap.add_argument("--cc", default=os.environ.get("CC", "cc")) + ap.add_argument("--threads", type=int, default=4) + ap.add_argument("--precision", type=int, default=160) + ap.add_argument("--rounds", type=int, default=3) + ap.add_argument("--kernel-target-calls", type=int, default=10_000_000) + ap.add_argument("--route-target-calls", type=int, default=1_000_000) + ap.add_argument("--skip-build", action="store_true") + args = ap.parse_args() + + repo = Path(args.repo).resolve() + out = Path(args.output_dir).resolve() + (out / "raw").mkdir(parents=True, exist_ok=True) + (out / "cases").mkdir(parents=True, exist_ok=True) + (out / "logs").mkdir(parents=True, exist_ok=True) + import time as _time + t_start = _time.time() + phase_times = {} + + # ---------------------------------------------------------------- cases + t0 = _time.time() + kernel_cases, route_cases = casesmod.build() + casesmod.write_header(kernel_cases, route_cases, out / "cases" / "quadratic_cases.h") + casesmod.write_json(kernel_cases, route_cases, out / "cases" / "cases.json") + phase_times["case_gen"] = _time.time() - t0 + print(f"cases: kernel={len(kernel_cases)} route={len(route_cases)}") + + # ------------------------------------------------------------ environment + env_info = buildmod.environment(repo, args.cc, args.threads) + env_info["python_version"] = sys.version.split()[0] + env_info["invocation"] = " ".join([sys.executable, *sys.argv]) + (out / "environment.json").write_text(json.dumps(env_info, indent=2) + "\n") + (out / "raw" / "invocation.txt").write_text( + " ".join([sys.executable, *sys.argv]) + "\n" + ) + print(f"source sha256: {env_info['source_sha256'][:16]} cc: {env_info['cc_version']}") + + # ---------------------------------------------------------------- build + t0 = _time.time() + if not args.skip_build: + info = buildmod.build_all(repo, out, args.cc, buildmod.BASE_FLAGS, args.threads) + exes = {k: Path(v) for k, v in info["executables"].items()} + else: + manifest, reasons = buildmod.verify_manifest(repo, out, buildmod.BASE_FLAGS, args.cc) + if reasons: + print("refusing --skip-build: existing build does not match the " + f"current source/flags/fixtures ({', '.join(reasons)}); " + "rerun without --skip-build", file=sys.stderr) + sys.exit(1) + exes = {k: Path(v) for k, v in manifest["executables"].items()} + phase_times["build"] = _time.time() - t0 + for name, exe in exes.items(): + if not exe.exists(): + print(f"missing executable for {name}: {exe}", file=sys.stderr) + sys.exit(1) + + # -------------------------------------------------------------- accuracy + t0 = _time.time() + env_log = out / "logs" / "probe_environment.log" + env_log.write_text("") + for name in VARIANTS: + kout = out / "raw" / f"kernel_{name}.csv" + rout = out / "raw" / f"route_{name}.csv" + proc = run_cmd( + [str(exes[name]), "accuracy", "--kernel-out", str(kout), + "--route-out", str(rout)], + out / "logs" / f"accuracy_{name}.log", + ) + with open(env_log, "a") as fh: + fh.write(proc.stdout) + if proc.returncode != 0: + print(f"accuracy probe failed for {name}", file=sys.stderr) + sys.exit(1) + phase_times["accuracy"] = _time.time() - t0 + + # ---------------------------------------------------------------- oracle + t0 = _time.time() + kernel_ref_cases, route_cases_json = oracle.load_cases(out / "cases" / "cases.json") + refs = { + cid: oracle.kernel_reference(case, args.precision) + for cid, case in kernel_ref_cases.items() + } + kernel_records = {} + route_records = {} + for name in VARIANTS: + krows = read_csv(out / "raw" / f"kernel_{name}.csv") + kernel_records[name] = oracle.analyze_kernel_variant(krows, refs, args.precision) + rrows = read_csv(out / "raw" / f"route_{name}.csv") + route_records[name] = oracle.analyze_route_variant(rrows, args.precision) + + kcols = ["id", "category", "ref_status", "variant_status", "root_err_ulps", + "root_rel_err", + "resid_kernel_rel", "resid_ideal_rel", "baseline_resid_rel", + "a_zero_ref", "a_zero_kernel", "a_sign_ref", "a_sign_kernel", + "disc_sign_ref", "near_linear", "near_boundary", + "ill_conditioned", "domain_clipped", "late"] + for name in VARIANTS: + write_metrics_csv(out / "raw" / f"kernel_metrics_{name}.csv", + kernel_records[name], kcols) + rcols = ["id", "category", "ref_status", "ref_domain_clipped", "status", + "kind", "kind_name", "failure_reason", "fallback", "pi_match", + "lcam_match", "kern_status", "kern_sigma", "F_over_tol", + "inside_mismatch", "false_miss", "false_candidate", "unconfirmed", + "kernel_false_miss", "kernel_false_miss_escaped"] + for name in VARIANTS: + write_metrics_csv(out / "raw" / f"route_metrics_{name}.csv", + route_records[name], rcols) + + ksummary = oracle.summarize_kernel(kernel_records, refs) + rsummary = oracle.summarize_route(route_records) + consistency_ids = { + rec["id"] for rec in kernel_records[VARIANTS[0]] + if rec["id"] < 1000 or rec["near_linear"] or rec["late"] + or rec["near_boundary"] or rec["domain_clipped"] + } + precision_mismatches = oracle.precision_consistency( + kernel_ref_cases, consistency_ids, args.precision, 240 + ) + curated_kernel = build_curated_kernel(kernel_records) + curated_route = build_curated_route(route_records) + write_curated(out / "raw" / "curated_kernel.csv", curated_kernel, + ["id", "category", "ref"] + [f"{v}_{f}" for v in VARIANTS + for f in ("status", "ulp", "rel")]) + write_curated(out / "raw" / "curated_route.csv", curated_route, + ["id", "category", "ref"] + [f"{v}_{f}" for v in VARIANTS + for f in ("kind", "fallback", + "F_over_tol")]) + + # ------------------------------------------------------------- hardware FMA + phase_times["oracle"] = _time.time() - t0 + fma_info = verify_fma(out, args.threads) + + # ---------------------------------------------------------------- timing + t0 = _time.time() + timings = timed_runs(exes, out, args, kernel_target_calls=args.kernel_target_calls, + route_target_calls=args.route_target_calls) + lin_target = max(100_000, args.kernel_target_calls // 20) + timings["linearity"] = linearity_check(exes, out, lin_target) + phase_times["timing"] = _time.time() - t0 + + # ---------------------------------------------------------------- summary + summary = { + "source_sha256": env_info["source_sha256"], + "flags": buildmod.BASE_FLAGS, + "precision": args.precision, + "precision_consistency": { + "ids_checked": len(consistency_ids), + "mismatches": precision_mismatches, + }, + "threads_timing": args.threads, + "cases": {"kernel": len(kernel_cases), "route": len(route_cases)}, + "kernel_accuracy": ksummary, + "route_accuracy": rsummary, + "curated_kernel": curated_kernel, + "fma": fma_info, + "timing": timings, + "phase_seconds": phase_times, + "total_seconds": _time.time() - t_start, + } + (out / "summary.json").write_text(json.dumps(summary, indent=2) + "\n") + print_summary(summary) + print(f"\nartifacts under {out}") + + +def verify_fma(out: Path, threads: int): + info = {"objects": {}} + for name in VARIANTS: + obj = out / "build" / f"probe_{name}.o" + if not obj.exists(): + continue + proc = subprocess.run(["objdump", "-d", str(obj)], capture_output=True, + text=True, check=False) + n = sum(1 for line in proc.stdout.splitlines() + if "vfmadd" in line or "vfmsub" in line or "fmadd" in line) + nm = subprocess.run(["nm", "-u", str(obj)], capture_output=True, + text=True, check=False) + undef = sorted({tok for line in nm.stdout.splitlines() + for tok in line.split() + if tok.startswith("fma")}) + kd = subprocess.run(["objdump", "-dr", str(obj)], capture_output=True, + text=True, check=False) + kb_lines = [] + inside = False + for line in kd.stdout.splitlines(): + if re.match(r"^[0-9a-f]+ 0, + "undefined_fma_symbols": undef, + "entry_solve_x87": x87, + "entry_solve_vfmadd": vfma, + "entry_solve_fmal_calls": calls_fmal, + } + print(f"arith verify {name}: fused={n} undef={undef} " + f"entry_solve[x87={x87} vfmadd={vfma} fmal_calls={calls_fmal}]") + (out / "raw" / "fma_verify.json").write_text(json.dumps(info, indent=2) + "\n") + return info + + +def linearity_check(exes, out: Path, target_calls: int): + """Confirm the microbench time scales with the call count (no hoisting).""" + res = {} + for name in ("ld_fma", "double_fma"): + secs = [] + for mult in (1, 2): + outp = out / "raw" / f"linearity_{name}_{mult}.json" + d = run_timing([str(exes[name]), "microbench", "--out", str(outp), + "--target-calls", str(target_calls * mult)], + outp, out / "logs" / f"linearity_{name}_{mult}.log") + secs.append(d["seconds"]) + res[name] = {"t1": secs[0], "t2": secs[1], + "ratio": secs[1] / secs[0] if secs[0] > 0 else None} + return res + + +def timed_runs(exes, out: Path, args, kernel_target_calls, route_target_calls): + schedules = [] + fwd = list(VARIANTS) + schedules.append(fwd) + schedules.append(list(reversed(fwd))) + for i in range(2, args.rounds): + schedules.append(fwd if i % 2 == 0 else list(reversed(fwd))) + schedules = schedules[: args.rounds] + + results = {"microbench": [], "routebench": []} + + def order(seq): + return [exes[n] for n in seq] + + for rnd, seq in enumerate(schedules): + for name in seq: + outp = out / "raw" / f"microbench_{name}_r{rnd}.json" + d = run_timing([str(exes[name]), "microbench", "--out", str(outp), + "--target-calls", str(kernel_target_calls)], + outp, out / "logs" / f"microbench_{name}_r{rnd}.log") + d["round"] = rnd + results["microbench"].append(d) + for rnd, seq in enumerate(schedules): + for name in seq: + outp = out / "raw" / f"routebench_{name}_r{rnd}.json" + d = run_timing([str(exes[name]), "routebench", "--out", str(outp), + "--target-calls", str(route_target_calls), + "--threads", str(args.threads), "--max-seconds", "20"], + outp, out / "logs" / f"routebench_{name}_r{rnd}.log") + d["round"] = rnd + results["routebench"].append(d) + return results + + +def print_summary(summary): + print("\n=== kernel accuracy (per variant) ===") + print(f"{'variant':<18}{'ENTER':>7}{'MISS':>7}{'UNC':>6}{'falseMiss':>10}" + f"{'falseCand':>10}{'domClipCand':>12}{'a0ref':>7}{'a0k':>6}{'asign':>6}" + f"{'illcond':>8}{'nearBnd':>8}{'medULP':>10}{'p95ULP':>11}{'maxULP':>11}") + for v, s in summary["kernel_accuracy"]["variants"].items(): + c = s["counts"] + d = s["root_ulp_all_ref_enter"] + print(f"{v:<18}{c['enter']:>7}{c['miss']:>7}{c['uncertain']:>6}" + f"{c['false_miss']:>10}{c['false_candidate']:>10}" + f"{c['domain_clipped_candidate']:>12}" + f"{c['a_zero_ref']:>7}{c['a_zero_kernel']:>6}" + f"{c['a_sign_mismatch']:>6}{c['ill_conditioned_enter']:>8}" + f"{c['near_boundary_enter']:>8}" + f"{d.get('median', float('nan')):>10.4g}{d.get('p95', float('nan')):>11.4g}" + f"{d.get('max', float('nan')):>11.4g}") + print("well-conditioned common-ENTER ULP (excludes late / near-linear /" + " growing-out / shrinking):") + for v, s in summary["kernel_accuracy"]["variants"].items(): + d = s["root_ulp_common_enter_wellcond"] + r = s["root_rel_common_enter_wellcond"] + print(f" {v:<18}n={d.get('n',0):>5} median={d.get('median', float('nan')):>9.3g}" + f" p95={d.get('p95', float('nan')):>9.3g}" + f" max={d.get('max', float('nan')):>9.3g}" + f" |rel err median={r.get('median', float('nan')):>9.3g}" + f" p95={r.get('p95', float('nan')):>9.3g}") + print(f"common reference-ENTER ids solved ENTRY by all variants: " + f"{summary['kernel_accuracy']['common_enter_ids']}") + print("\n=== curated kernel cases (status/ULP) ===") + print(f"{'id':>4} {'category':<20} {'ref':<6} " + + " ".join(f"{v:>16}" for v in VARIANTS)) + for row in summary.get("curated_kernel", []): + cells = [] + for v in VARIANTS: + st = row.get(f"{v}_status") + u = row.get(f"{v}_ulp") + cells.append(f"{st}/{float(u):.3g}" if u not in (None, "") else f"{st}/-") + print(f"{row['id']:>4} {row['category']:<20} {row.get('ref',''):<6} " + + " ".join(f"{c:>16}" for c in cells)) + print("\n=== route accuracy (per variant) ===") + print(f"{'variant':<18}{'insideOK':>9}{'insideBad':>10}{'falseMiss':>10}" + f"{'falseCand':>10}{'unconf':>8}{'kernFM':>8}{'kernFMesc':>10}" + f"{'fallback':>9}{'F>tol':>7}{'maxF/tol':>10}") + for v, s in summary["route_accuracy"]["variants"].items(): + c = s["counts"] + print(f"{v:<18}{c['inside_ok']:>9}{c['inside_mismatch']:>10}" + f"{c['false_miss']:>10}{c['false_candidate']:>10}" + f"{c['unconfirmed']:>8}{c['kernel_false_miss']:>8}" + f"{c['kernel_false_miss_escaped']:>10}" + f"{c['fallback_used']:>9}" + f"{c['F_over_tol_gt1']:>7}{s['max_F_over_tol']:>10.3g}") + pc = summary.get("precision_consistency", {}) + print(f"reference precision {summary['precision']} vs 240: checked " + f"{pc.get('ids_checked')} curated/near-linear/late ids, " + f"mismatches={len(pc.get('mismatches', []))}") + print("\nroute fallback / entry / escaped by category " + "(ld_plain | ld_fma | double_fma):") + ra = summary["route_accuracy"]["variants"] + cats = sorted({c for s in ra.values() for c in s.get("by_category", {})}) + for cat in cats: + cells = [] + for v in ("ld_plain", "ld_fma", "double_fma"): + d = ra[v].get("by_category", {}).get(cat, {}) + cells.append(f"fb={d.get('fallback',0)} e={d.get('entry',0)} " + f"x={d.get('escaped',0)}") + print(f" {cat:<20} " + " | ".join(cells)) + print("\n=== timing (median ns/call over rounds) ===") + for mode in ("microbench", "routebench"): + per = {} + for rec in summary["timing"][mode]: + per.setdefault(rec["variant"], []).append(rec["ns_per_call"]) + for v, vals in per.items(): + vals.sort() + med = vals[len(vals) // 2] + print(f"{mode:<12}{v:<18}median={med:>10.2f} ns/call " + f"rounds={len(vals)}") + lin = summary["timing"].get("linearity", {}) + for v, d in lin.items(): + print(f"linearity {v:<18}t1={d['t1']:.4f}s t2={d['t2']:.4f}s " + f"ratio={d['ratio']:.3f} (expect ~2)") + + +if __name__ == "__main__": + main()