From 699c5bf4ca5fb1fa6dc9c9c66a722b6d7b61c3bd Mon Sep 17 00:00:00 2001 From: Nicolas Maman Date: Tue, 20 Oct 2026 13:58:45 -0300 Subject: [PATCH 1/6] Compare the SIMD widths where the AVX2 copy is plain code (#184) test_simd_widths checked contact_solver_wide8.lanes_native() before comparing the widths. Aether does not promise a native copy: GCC on Windows at -O0 (aether#3023) and the AETHER_GCC_COMPAT=0 path compile it as plain code, where native8() is false and the copy still has to give the same bits. The check failed those builds although they are correct. The scenes are now compared at both widths in every build, and the test says when the copy is plain. Measured with Aether at bd667cb1 on Windows: 29 checks pass at -O3; at -O0 the widths agree bit for bit and the note prints. --- README.md | 2 +- tests/test_determinism.ae | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index bb235a1..57f4819 100644 --- a/README.md +++ b/README.md @@ -170,7 +170,7 @@ compared bit for bit; the query spawn sleeps a step after the reference's with its 59 query hits, the mesh drop two; and its rolling mix of spheres, capsules and boxes with rolling resistance, whose trace folds the manifolds' impulses in, the scene of its SIMD width test, which here steps every scene at width -4 and at the AVX2 copy's 8 in one process, bit for bit; 30 checks) and `test_large_world.ae` (a stack, a bullet and the +4 and at the AVX2 copy's 8 in one process, bit for bit; 29 checks) and `test_large_world.ae` (a stack, a bullet and the origin-relative queries at x = 0 and at x = 1e7 agree, the whole engine being in doubles; 43 checks), `test_collision.ae` (the reference's test_collision.c: boxes valid and not, a ray against a box in its twelve diff --git a/tests/test_determinism.ae b/tests/test_determinism.ae index 99a28e6..a8f1aa0 100644 --- a/tests/test_determinism.ae +++ b/tests/test_determinism.ae @@ -651,7 +651,10 @@ rolling_mix(workers: int) -> Outcome { // Box3D's SIMDWidthTest: every scene stepped by the wide stages' copy for // AVX2 and by the build's own copy, in one process, gives the same sleep // step, checksum and query results. Where the CPU has no AVX2 there is -// one width and nothing to compare. +// one width and nothing to compare. Some builds compile the copy as plain +// code (GCC on Windows at -O0, aether#3023, and the AETHER_GCC_COMPAT=0 +// path): it still has to give the same bits, so the scenes are compared +// there too (#184). test_simd_widths() { basics.set_simd_width(0) native_width = basics.get_simd_width() @@ -659,7 +662,9 @@ test_simd_widths() { println("determinism: simd widths skipped, this CPU runs one width") return } - ensure("simd widths: the AVX2 copy holds eight lanes a register", contact_solver_wide8.lanes_native()) + if contact_solver_wide8.lanes_native() == false { + println("determinism: simd widths: this build compiles the AVX2 copy as plain code") + } basics.set_simd_width(4) a4 = falling_ragdolls(1) b4 = wave_pile(1) From 902e781e422bc22c4699b1b53f6e826c49fb39ef Mon Sep 17 00:00:00 2001 From: Nicolas Maman Date: Tue, 20 Oct 2026 14:00:36 -0300 Subject: [PATCH 2/6] A Windows leg, with the determinism scenes at -O0 (#184) ae3d and Siltreach run on Windows, and ae3d builds its suites at -O0, but CI here built only Linux and macOS at -O3. Two bugs reached ae3d first: GCC on Win64 crashing on the AVX2 copies at -O0 (aether#3023), and four workers parting at -O0, where every sqrt is MinGW's library call and that rounds twice on the main thread's x87 precision. The new job runs the suite on windows-latest under MSYS2 UCRT64, with the release fetched as its .zip (get.sh looks for a .tar.gz), then test_determinism at -O0 against the golden trace. The pins move to one workflow-level env, read by get.sh, so a bump is one line. The README names the leg, and says five scenes, not four: the rolling mix joined the trace and the count stayed. --- .github/workflows/test.yaml | 76 ++++++++++++++++++++++++++++++++----- README.md | 8 ++-- 2 files changed, 72 insertions(+), 12 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 768f770..264c9b1 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -7,6 +7,15 @@ on: - main pull_request: +# Pin: 0.835.0, the latest release: its lane moves (aether#2550), bit +# casts and reductions (#2603) and `@target("avx2") import` (#2607) are +# what the wide gather, the wide separating axis test and the run-time +# width are written with. A reproducible bootstrap version: the floor +# this code needs, verified by every job below. +env: + AE_PIN: 0.835.0 + AEB_REF: v0.328 + jobs: test: name: aephysics tests @@ -15,15 +24,9 @@ jobs: - name: Checkout aephysics uses: actions/checkout@v4 - # Pin: 0.835.0, the latest release: its lane moves (aether#2550), bit - # casts and reductions (#2603) and `@target("avx2") import` (#2607) are - # what the wide gather, the wide separating axis test and the run-time - # width are written with. A reproducible bootstrap - # version: the floor this code needs, verified on Linux by this job. - name: Install the released Aether toolchain and aeb run: | - curl -fsSL https://raw.githubusercontent.com/aether-lang-dev/aeb/main/get.sh \ - | AE_PIN=0.835.0 AEB_REF=v0.328 sh + curl -fsSL https://raw.githubusercontent.com/aether-lang-dev/aeb/main/get.sh | sh echo "$HOME/.local/bin" >> "$GITHUB_PATH" - name: Run the tests with ae @@ -55,9 +58,64 @@ jobs: - name: Install the released Aether toolchain and aeb run: | - curl -fsSL https://raw.githubusercontent.com/aether-lang-dev/aeb/main/get.sh \ - | AE_PIN=0.835.0 AEB_REF=v0.328 sh + curl -fsSL https://raw.githubusercontent.com/aether-lang-dev/aeb/main/get.sh | sh echo "$HOME/.local/bin" >> "$GITHUB_PATH" - name: Run the tests with ae run: scripts/test.sh + + # Windows, where ae3d and Siltreach run, and -O0, which ae3d builds its + # suites at. Neither was built here, and two bugs reached ae3d first + # (#184): GCC on Win64 crashing on the AVX2 copies at -O0 + # (aether#3023), and four workers parting at -O0 because MinGW's + # library sqrt rounds twice on the main thread's x87 setting. + windows: + name: aephysics tests (Windows x86-64) + runs-on: windows-latest + timeout-minutes: 45 + defaults: + run: + shell: msys2 {0} + steps: + - name: Checkout aephysics + uses: actions/checkout@v4 + + # UCRT64: the release's libaether.a links the Universal CRT, and does + # not link against an msvcrt gcc. + - name: Set up MSYS2 + uses: msys2/setup-msys2@v2 + with: + msystem: UCRT64 + install: >- + unzip + mingw-w64-ucrt-x86_64-gcc + + # The Windows release is a .zip, and get.sh looks for a .tar.gz (it + # would build the compiler from source instead), so it is fetched + # here. + - name: Install the released Aether toolchain + run: | + mkdir -p target + curl -fsSL -o target/aether.zip \ + "https://github.com/aether-lang-dev/aether/releases/download/v$AE_PIN/aether-$AE_PIN-windows-x86_64.zip" + mkdir -p "$HOME/aether" && unzip -q target/aether.zip -d "$HOME/aether" + "$HOME/aether/bin/ae" version + + # PATH is set in each step: setup-msys2's shell does not read what + # GITHUB_PATH adds. + - name: Run the tests with ae + run: | + export PATH="$HOME/aether/bin:$PATH" + scripts/test.sh + + # Unoptimized, every sqrt is a library call and every function holds + # its lanes on the stack: the scenes keep their bits over four + # workers, and the trace matches the golden one. + - name: The determinism scenes at -O0 + env: + AETHER_AEPHYSICS_CFLAGS: -O0 + run: | + export PATH="$HOME/aether/bin:$PATH" + ae build tests/test_determinism.ae -o target/test_determinism + target/test_determinism + scripts/golden_trace.sh check diff --git a/README.md b/README.md index 57f4819..2790c87 100644 --- a/README.md +++ b/README.md @@ -239,9 +239,11 @@ The step is deterministic across platforms and worker counts, not only across runs: `AEPHYSICS_TRACE=1 target/test_determinism` prints every body's checksum after every step with its bits, and `scripts/golden_trace.sh check` (run by `scripts/test.sh`) diffs that -against `tests/golden/determinism_trace.txt` on every step of the four -scenes. CI checks it on Linux, on Linux with `-mavx2 -mfma` and on macOS -arm64; Windows (MinGW) gives the same bits. The engine's own maths (the +against `tests/golden/determinism_trace.txt` on every step of the five +scenes. CI checks it on Linux, on Linux with `-mavx2 -mfma`, on macOS +arm64 and on Windows x86-64 (MinGW), there also unoptimized (#184): at +`-O0` every `sqrt` is the C library's, and the library has to round it +exactly, as the instruction does. The engine's own maths (the reference's cos, sin and atan2 approximations, `sqrt`, `floor` and `remainder` as the only libm calls) and `-ffp-contract=off` are what make that hold: without the flag GCC and clang fuse multiply-adds wherever the From 41eaddd93be9ab0b4966530c7bdd56f976ec27bb Mon Sep 17 00:00:00 2001 From: Nicolas Maman Date: Tue, 20 Oct 2026 14:26:42 -0300 Subject: [PATCH 3/6] The engine's double square root is std.math's, exact at every -O (#184) At -O0 every sqrt is a library call, and MinGW's library sqrt is the x87 one at the 64-bit precision a MinGW main thread runs with: it rounds twice, one ulp off for 253 of 1,000,000 inputs. From -O1 up GCC inlines sqrtsd and the library is never reached, so only unoptimized builds on Windows parted from the golden trace (wave pile from step 10). -O3 with -fno-builtin-sqrt gives exactly the -O0 bits, which is the proof. math.sqrt_float now calls std.math's sqrt, which runs in Aether's optimized runtime as the instruction, and every module, test and bench calls math.sqrt_float: the per-module `extern sqrt` lines are gone, so a bare sqrt no longer compiles. libm's sqrtf stays: its 64-bit x87 result rounds to a float once, and no float of 304 million differs from sqrtss. The call costs nothing measurable (bench/RESULTS.md: three scenes, alternated runs, within their spread) and the bits are the same: the suite and the golden trace pass at -O3 with Aether 0.835.0. bench/RESULTS.md also drops the slice indexing from what is left in the SAT: the Aether session measured every bounds test folded and the samples on aether_slice.h as skid from the loads (aether#3020). --- README.md | 12 +++++----- aephysics/capsule/module.ae | 7 +++--- aephysics/contact_solver/module.ae | 5 ++--- aephysics/distance/module.ae | 6 ++--- aephysics/dynamics/module.ae | 17 +++++++-------- aephysics/hull/module.ae | 5 ++--- aephysics/human/module.ae | 4 ++-- aephysics/joint_solver/module.ae | 5 ++--- aephysics/math/module.ae | 35 ++++++++++++++++++++---------- aephysics/sphere/module.ae | 5 ++--- bench/RESULTS.md | 29 +++++++++++++++++++++++-- bench/hull.ae | 5 ++--- bench/suite.ae | 4 ++-- tests/test_determinism.ae | 9 ++++---- tests/test_distance.ae | 3 +-- tests/test_height_field.ae | 1 - tests/test_hull.ae | 9 ++++---- tests/test_manifold.ae | 1 - tests/test_mesh.ae | 1 - tests/test_physics_world.ae | 5 ++--- tests/test_restitution.ae | 3 +-- tests/test_shape.ae | 7 +++--- tests/test_soak.ae | 3 +-- tests/test_triangle_manifold.ae | 3 +-- 24 files changed, 101 insertions(+), 83 deletions(-) diff --git a/README.md b/README.md index 2790c87..c4bdf11 100644 --- a/README.md +++ b/README.md @@ -241,12 +241,12 @@ every body's checksum after every step with its bits, and `scripts/golden_trace.sh check` (run by `scripts/test.sh`) diffs that against `tests/golden/determinism_trace.txt` on every step of the five scenes. CI checks it on Linux, on Linux with `-mavx2 -mfma`, on macOS -arm64 and on Windows x86-64 (MinGW), there also unoptimized (#184): at -`-O0` every `sqrt` is the C library's, and the library has to round it -exactly, as the instruction does. The engine's own maths (the -reference's cos, sin and atan2 approximations, `sqrt`, `floor` and -`remainder` as the only libm calls) and `-ffp-contract=off` are what make -that hold: without the flag GCC and clang fuse multiply-adds wherever the +arm64 and on Windows x86-64 (MinGW), there also unoptimized (#184). The +engine's own maths (the reference's cos, sin and atan2 approximations; +`sqrtf`, `floor` and `remainder` as the only libm calls; std.math's +`sqrt`, exact at every `-O` level, where MinGW's library `sqrt` rounds +twice at `-O0`) and `-ffp-contract=off` are what make that hold: without +the flag GCC and clang fuse multiply-adds wherever the target has FMA, and arm64 parted from Linux on the first step (#125). A program that builds aephysics into itself must pass `-ffp-contract=off` in its own `[build] cflags`, as `ae` reads only the root manifest's. diff --git a/aephysics/capsule/module.ae b/aephysics/capsule/module.ae index 4bc80bd..91796fa 100644 --- a/aephysics/capsule/module.ae +++ b/aephysics/capsule/module.ae @@ -21,7 +21,6 @@ exports ( ray_cast_capsule, shape_cast_capsule, collide_mover_and_capsule ) -extern sqrt(x: float) -> float // The capsule's mass: a cylinder and a sphere (the two caps) about the // midpoint, the caps carried out to the ends by Steiner, then turned @@ -93,7 +92,7 @@ ray_cast_capsule(c: Capsule, origin: Vec3, translation: Vec3, max_fraction: floa } // The ray's origin from the first centre, projected onto the axis. s = math.sub(origin, c1) - length = sqrt(length_squared) + length = math.sqrt_float(length_squared) axis = math.mul_sv(1.0 / length, d) u = math.dot(s, axis) c_axis = math.mul_sv(u, axis) @@ -137,7 +136,7 @@ ray_cast_capsule(c: Capsule, origin: Vec3, translation: Vec3, max_fraction: floa // Casting away from the axis, or never closing to the radius. if beta >= 0.0 || disc < 0.0 { return output } // The near root in the form that avoids the cancellation as the ray nears parallel. - tr = gamma / (0.0 - beta + sqrt(disc)) + tr = gamma / (0.0 - beta + math.sqrt_float(disc)) } else { // The closest points of the infinite ray and the infinite axis. inv_det = 1.0 / det @@ -151,7 +150,7 @@ ray_cast_capsule(c: Capsule, origin: Vec3, translation: Vec3, max_fraction: floa g2 = math.length_squared(g) if g2 > r * r { return output } // The ray into the cylinder, relative to the closest point (like the sphere). - h = sqrt((r * r - g2) * inv_det) + h = math.sqrt_float((r * r - g2) * inv_det) tr = t2 - h } if tr < 0.0 || max_fraction * ray_length < tr { return output } diff --git a/aephysics/contact_solver/module.ae b/aephysics/contact_solver/module.ae index b781495..e1d890f 100644 --- a/aephysics/contact_solver/module.ae +++ b/aephysics/contact_solver/module.ae @@ -32,7 +32,6 @@ exports ( color_constraint_counts, count_color_constraints ) -extern sqrt(x: float) -> float const MIN_FRICTION_WEIGHT = 0.0000000001 const NULL_INDEX = 0 - 1 @@ -512,7 +511,7 @@ solve_contacts(ctx: *StepContext, constraints: ptr, count: int, use_bias: bool) max_rolling = rolling_resistance * total_normal_impulse mag_sqr = math.dot(constraint.rolling_impulse, constraint.rolling_impulse) if mag_sqr > max_rolling * max_rolling + math.EPSILON { - constraint.rolling_impulse = math.mul_sv(max_rolling / sqrt(mag_sqr), constraint.rolling_impulse) + constraint.rolling_impulse = math.mul_sv(max_rolling / math.sqrt_float(mag_sqr), constraint.rolling_impulse) } delta_rolling = math.sub(constraint.rolling_impulse, old_rolling) w_a = math.sub(w_a, math.mul_mv(i_a, delta_rolling)) @@ -532,7 +531,7 @@ solve_contacts(ctx: *StepContext, constraints: ptr, count: int, use_bias: bool) max_friction = friction * total_normal_impulse length_squared = math.dot2(new_friction, new_friction) if length_squared > max_friction * max_friction { - scale = max_friction / sqrt(length_squared) + scale = max_friction / math.sqrt_float(length_squared) new_friction = math.mul_sv2(scale, new_friction) } delta_friction = math.sub2(new_friction, constraint.friction_impulse) diff --git a/aephysics/distance/module.ae b/aephysics/distance/module.ae index 4b86355..ead8c06 100644 --- a/aephysics/distance/module.ae +++ b/aephysics/distance/module.ae @@ -719,7 +719,7 @@ segment_proxy_distance(input: *DistanceInput, cache: *SimplexCache) -> DistanceO gap_sq = math.dot(gap, gap) // Touching cores: no normal, no radii, as GJK's overlap. if gap_sq < 1000.0 * FLOAT_MIN { return output } - length = sqrt(gap_sq) + length = math.sqrt_float(gap_sq) normal = math.mul_sv(1.0 / length, gap) output.distance = length output.normal = normal @@ -916,7 +916,7 @@ triangle_point_distance(input: *DistanceInput, cache: *SimplexCache, out: *Dista } cache.metric = 0.0 if count == 2 { cache.metric = math.distance(feature_point(near.f0, t0, t1, t2), feature_point(near.f1, t0, t1, t2)) } - if count == 3 { cache.metric = sqrt(area_sq) / 2.0 } + if count == 3 { cache.metric = math.sqrt_float(area_sq) / 2.0 } if triangle_is_a { write_output(out, input, near.point, p, true) } else { write_output(out, input, p, near.point, true) } return true } @@ -971,7 +971,7 @@ write_output(out: *DistanceOutput, input: *DistanceInput, point_a: Vec3, point_b gap = math.sub(point_b, point_a) gap_sq = math.dot(gap, gap) if separate == false || gap_sq < 1000.0 * FLOAT_MIN { return } - length = sqrt(gap_sq) + length = math.sqrt_float(gap_sq) normal = math.mul_sv(1.0 / length, gap) out.distance = length out.normal = normal diff --git a/aephysics/dynamics/module.ae b/aephysics/dynamics/module.ae index f8e8621..8a3ced6 100644 --- a/aephysics/dynamics/module.ae +++ b/aephysics/dynamics/module.ae @@ -166,7 +166,6 @@ exports ( extern memset(block: ptr, value: int, size: int) -> ptr extern memcpy(dst: ptr, src: ptr, size: int) -> ptr extern memcmp(a: ptr, b: ptr, size: int) -> int -extern sqrt(x: float) -> float const NULL_INDEX = 0 - 1 const MAX_WORLDS = 128 @@ -390,7 +389,7 @@ struct ShapeId { generation: int } -default_friction(friction_a: float, user_a: long, friction_b: float, user_b: long) -> float { return sqrt(friction_a * friction_b) } +default_friction(friction_a: float, user_a: long, friction_b: float, user_b: long) -> float { return math.sqrt_float(friction_a * friction_b) } default_restitution(restitution_a: float, user_a: long, restitution_b: float, user_b: long) -> float { return math.max_float(restitution_a, restitution_b) } default_world_def() -> WorldDef { @@ -2201,7 +2200,7 @@ compute_shape_margin(g: *Shape) -> float { max_extent_sqr = math.max_float(max_extent_sqr, math.distance_squared(points[i], g.hull.center)) i = i + 1 } - margin = sqrt(max_extent_sqr) + margin = math.sqrt_float(max_extent_sqr) } else { // Static-only shapes: the broad phase uses the speculative distance for static proxies. return MAX_AABB_MARGIN @@ -6007,7 +6006,7 @@ try_recycle_contact(w: *World, c: *Contact, sim_a: *BodySim, sim_b: *BodySim, st max_extent = math.max_vec3(max_extent_a, max_extent_b) dist_squared = math.distance_squared(xf.p, xfc.p) if angular_distance <= CONTACT_RECYCLE_ANGULAR_DISTANCE || dist_squared >= tolerance * tolerance { return false } - slack = tolerance - sqrt(dist_squared) + slack = tolerance - math.sqrt_float(dist_squared) qr = math.inv_mul_quat(xfc.q, xf.q) a = math.abs_vec3(math.qv(qr)) arc = math.vec3(a.y * max_extent.z + a.z * max_extent.y, a.z * max_extent.x + a.x * max_extent.z, a.x * max_extent.y + a.y * max_extent.x) @@ -7079,7 +7078,7 @@ cone_reach(s: *SphericalJoint, b: Vec3) -> ConeReach { // The ellipse's half-angle in a direction: the polar radius of an ellipse // of half-angles a (along the major axis) and b (across it). cone_radius(a: float, b: float, along_major: float, along_minor: float) -> float { - return a * b / sqrt((b * along_major) * (b * along_major) + (a * along_minor) * (a * along_minor)) + return a * b / math.sqrt_float((b * along_major) * (b * along_major) + (a * along_minor) * (a * along_minor)) } // How far frame B's z is inside the cone, radians: the limit's separation @@ -7394,13 +7393,13 @@ get_joint_reaction(w: *World, sim: *JointSim, inv_time_step: float) -> JointReac } else if kind == JOINT_WHEEL { wheel = wheel_joint_of(sim) axial = wheel.suspension_spring_impulse + wheel.lower_suspension_impulse - wheel.upper_suspension_impulse - linear_impulse = sqrt(wheel.linear_impulse.x * wheel.linear_impulse.x + wheel.linear_impulse.y * wheel.linear_impulse.y + axial * axial) + linear_impulse = math.sqrt_float(wheel.linear_impulse.x * wheel.linear_impulse.x + wheel.linear_impulse.y * wheel.linear_impulse.y + axial * axial) // The spin with the lock (unsteered) or the steering, on their // (orthogonal at rest) axes; the spin alone missed a wheel wrenched sideways. lock_x = wheel.angular_impulse.x lock_y = wheel.angular_impulse.y if wheel.enable_steering { lock_y = wheel.steering_spring_impulse + wheel.lower_steering_impulse - wheel.upper_steering_impulse } - angular_impulse = sqrt(lock_x * lock_x + lock_y * lock_y + wheel.spin_impulse * wheel.spin_impulse) + angular_impulse = math.sqrt_float(lock_x * lock_x + lock_y * lock_y + wheel.spin_impulse * wheel.spin_impulse) } return JointReaction { force: linear_impulse * inv_time_step, torque: angular_impulse * inv_time_step } } @@ -7462,7 +7461,7 @@ line_separation(dp: Vec3, axis: Vec3, enable_limit: bool, lower: float, upper: f if translation < lower { limit = lower - translation } if upper < translation { limit = translation - upper } } - return sqrt(across + limit * limit) + return math.sqrt_float(across + limit * limit) } // The same as an angle (the wheel joint's is a todo in the reference: zero @@ -7531,7 +7530,7 @@ wheel_angular_separation(wheel: *WheelJoint, q_a: Quat, q_b: Quat, rel_q: Quat) matrix_a = math.make_matrix_from_quat(q_a) matrix_b = math.make_matrix_from_quat(q_b) toward = math.clamp_float(math.abs_float(math.dot(matrix_b.cz, matrix_a.cx)), 0.0, 1.0) - lean = math.atan2_approx(toward, sqrt(1.0 - toward * toward)) + lean = math.atan2_approx(toward, math.sqrt_float(1.0 - toward * toward)) excess = 0.0 if wheel.enable_steering_limit { angle = math.atan2_approx(0.0 - math.dot(matrix_b.cz, matrix_a.cy), math.dot(matrix_b.cz, matrix_a.cz)) diff --git a/aephysics/hull/module.ae b/aephysics/hull/module.ae index 34810c1..18fe30f 100644 --- a/aephysics/hull/module.ae +++ b/aephysics/hull/module.ae @@ -45,7 +45,6 @@ extern memcpy(dst: ptr, src: ptr, size: int) -> ptr extern memmove(dst: ptr, src: ptr, size: int) -> ptr extern memset(block: ptr, value: int, size: int) -> ptr extern memcmp(a: ptr, b: ptr, size: int) -> int -extern sqrt(x: float) -> float const NULL_INDEX = 0 - 1 const MAX_HULL_VERTICES = 128 @@ -1798,7 +1797,7 @@ create_rock(radius: float) -> *HullData { count = 10 block = basics.alloc(count * sizeof(Vec3)) points = block as Vec3[] - phi = (1.0 + sqrt(5.0)) / 2.0 + phi = (1.0 + math.sqrt_float(5.0)) / 2.0 theta = 2.0 * math.PI / phi cosine = 1.0 sine = 0.0 @@ -1806,7 +1805,7 @@ create_rock(radius: float) -> *HullData { i = 0 while i < count { z = 1.0 - (2.0 * (i as float) + 1.0) / (count as float) - radius_xy = sqrt(1.0 - z * z) + radius_xy = math.sqrt_float(1.0 - z * z) points[i] = math.vec3(radius * radius_xy * cosine, radius * radius_xy * sine, radius * z) c0 = cosine s0 = sine diff --git a/aephysics/human/module.ae b/aephysics/human/module.ae index 03fbed4..dabd495 100644 --- a/aephysics/human/module.ae +++ b/aephysics/human/module.ae @@ -580,8 +580,8 @@ aim_ellipse(shape: ptr, index: int, up: Vec3, forward: Vec3, ahead: float, back: minor_degrees = hs iteration = 0 while iteration < 32 { - minor_degrees = hs / sqrt(math.max_float(1.0 - (cf * cf) / (major_degrees * major_degrees), 0.01)) - major_degrees = hf / sqrt(math.max_float(1.0 - (cs * cs) / (minor_degrees * minor_degrees), 0.01)) + minor_degrees = hs / math.sqrt_float(math.max_float(1.0 - (cf * cf) / (major_degrees * major_degrees), 0.01)) + major_degrees = hf / math.sqrt_float(math.max_float(1.0 - (cs * cs) / (minor_degrees * minor_degrees), 0.01)) iteration = iteration + 1 } spec.swing_degrees = major_degrees diff --git a/aephysics/joint_solver/module.ae b/aephysics/joint_solver/module.ae index 6903328..1552ae3 100644 --- a/aephysics/joint_solver/module.ae +++ b/aephysics/joint_solver/module.ae @@ -87,7 +87,6 @@ exports ( draw_joint ) -extern sqrt(x: float) -> float const NULL_INDEX = 0 - 1 @@ -875,7 +874,7 @@ solve_parallel_joint(base: *JointSim, ctx: *StepContext) { y: 0.0 - j.softness.mass_scale * sol.y - j.softness.impulse_scale * old_impulse.y } j.perp_impulse = math.add2(old_impulse, delta) if math.dot2(j.perp_impulse, j.perp_impulse) > max_impulse * max_impulse { - j.perp_impulse = math.mul_sv2(max_impulse / sqrt(math.dot2(j.perp_impulse, j.perp_impulse)), j.perp_impulse) + j.perp_impulse = math.mul_sv2(max_impulse / math.sqrt_float(math.dot2(j.perp_impulse, j.perp_impulse)), j.perp_impulse) } delta = math.sub2(j.perp_impulse, old_impulse) angular = math.blend2(delta.x, perp_x, delta.y, perp_y) @@ -1660,7 +1659,7 @@ prepare_spherical_joint(base: *JointSim, ctx: *StepContext) { } if j.enable_twist_limit { rel_q = math.inv_mul_quat(j.frame_a.q, j.frame_b.q) - tan_half = sqrt((rel_q.x * rel_q.x + rel_q.y * rel_q.y) / (rel_q.z * rel_q.z + rel_q.w * rel_q.w)) + tan_half = math.sqrt_float((rel_q.x * rel_q.x + rel_q.y * rel_q.y) / (rel_q.z * rel_q.z + rel_q.w * rel_q.w)) swing_axis = math.normalize(math.cross(cone_axis, twist_axis)) perp = math.cross(swing_axis, cone_axis) jacobian = math.mul_add(cone_axis, tan_half, perp) diff --git a/aephysics/math/module.ae b/aephysics/math/module.ae index 47ae306..90db442 100644 --- a/aephysics/math/module.ae +++ b/aephysics/math/module.ae @@ -21,6 +21,7 @@ // is_normalized, ...) a caller checks where it matters, and the functions // return what the original returns on the degenerate case. import std.string +import std.math as stdmath import aephysics.native exports ( @@ -70,9 +71,11 @@ exports ( vec3_text, quat_text ) -extern sqrt(x: float) -> float // libm's single-precision square root: the one the f32 family uses, so a -// length is computed at the width it is stored at. +// length is computed at the width it is stored at. Exact at every -O +// level: the x87 sqrt MinGW's library has at -O0 keeps 64 bits, enough +// for a float to round once (no float of 304 million differs from the +// SSE instruction's). extern sqrtf(x: f32) -> f32 extern remainder(x: float, y: float) -> float @@ -230,7 +233,15 @@ clamp_float(a: float, lower: float, upper: float) -> float { return a } lerp_float(a: float, b: float, alpha: float) -> float { return (1.0 - alpha) * a + alpha * b } -sqrt_float(a: float) -> float { return sqrt(a) } +// The square root of every double in the engine, rounded exactly, as +// IEEE 754 has it and the golden trace needs at every -O level (#184). +// It is std.math's, which runs in Aether's optimized runtime as the +// instruction. libm's own, bound here before, is the instruction too +// once GCC inlines it, but at -O0 every call reaches the library, and +// MinGW's is the x87 sqrt at 64-bit precision, which rounds twice: one +// ulp off for 1 input in 4,000. The call costs nothing measurable +// (bench/RESULTS.md). +sqrt_float(a: float) -> float { return stdmath.sqrt(a) } // An approximate arctangent in [-pi, pi], hand coded for cross-platform // determinism (the standard library's is not). Accurate to about 0.0023 @@ -292,7 +303,7 @@ compute_cos_sin(radians: float) -> CosSin { s = 16.0 * x * (PI - x) / (5.0 * pi2 - 4.0 * x * (PI - x)) } - mag = sqrt(s * s + c * c) + mag = sqrt_float(s * s + c * c) inv_mag = 0.0 if mag > 0.0 { inv_mag = 1.0 / mag } return CosSin { cosine: c * inv_mag, sine: s * inv_mag } @@ -312,7 +323,7 @@ sub(a: Vec3, b: Vec3) -> Vec3 { return Vec3 { x: a.x - b.x, y: a.y - b.y, z: a.z mul(a: Vec3, b: Vec3) -> Vec3 { return Vec3 { x: a.x * b.x, y: a.y * b.y, z: a.z * b.z } } neg(a: Vec3) -> Vec3 { return Vec3 { x: 0.0 - a.x, y: 0.0 - a.y, z: 0.0 - a.z } } dot(a: Vec3, b: Vec3) -> float { return a.x * b.x + a.y * b.y + a.z * b.z } -length(v: Vec3) -> float { return sqrt(dot(v, v)) } +length(v: Vec3) -> float { return sqrt_float(dot(v, v)) } length_squared(a: Vec3) -> float { return a.x * a.x + a.y * a.y + a.z * a.z } distance(a: Vec3, b: Vec3) -> float { return length(sub(b, a)) } distance_squared(a: Vec3, b: Vec3) -> float { return length_squared(sub(b, a)) } @@ -321,7 +332,7 @@ distance_squared(a: Vec3, b: Vec3) -> float { return length_squared(sub(b, a)) } normalize(a: Vec3) -> Vec3 { ls = a.x * a.x + a.y * a.y + a.z * a.z if ls > 1000.0 * TINY { - s = 1.0 / sqrt(ls) + s = 1.0 / sqrt_float(ls) return Vec3 { x: s * a.x, y: s * a.y, z: s * a.z } } return vec3_zero() @@ -332,7 +343,7 @@ normalize(a: Vec3) -> Vec3 { get_length_and_normalize(length_out: float[], a: Vec3) -> Vec3 { ls = a.x * a.x + a.y * a.y + a.z * a.z if ls > 1000.0 * TINY { - l = sqrt(ls) + l = sqrt_float(ls) length_out[0] = l s = 1.0 / l return Vec3 { x: s * a.x, y: s * a.y, z: s * a.z } @@ -505,7 +516,7 @@ negate_quat(q: Quat) -> Quat { return quat_vs(neg(qv(q)), 0.0 - q.w) } normalize_quat(q: Quat) -> Quat { ls = dot_quat(q, q) if ls > 1000.0 * TINY { - s = 1.0 / sqrt(ls) + s = 1.0 / sqrt_float(ls) return quat_vs(Vec3 { x: s * q.x, y: s * q.y, z: s * q.z }, s * q.w) } return quat_identity() @@ -520,7 +531,7 @@ make_quat_from_axis_angle(axis: Vec3, radians: float) -> Quat { // The axis of a unit quaternion, and in `radians_out[0]` its angle; a zero // axis for no rotation. get_axis_angle(radians_out: float[], q: Quat) -> Vec3 { - l = sqrt(q.x * q.x + q.y * q.y + q.z * q.z) + l = sqrt_float(q.x * q.x + q.y * q.y + q.z * q.z) radians_out[0] = 2.0 * atan2_approx(l, q.w) if l > 0.0 { inv = 1.0 / l @@ -530,7 +541,7 @@ get_axis_angle(radians_out: float[], q: Quat) -> Vec3 { } get_quat_angle(q: Quat) -> float { - l = sqrt(q.x * q.x + q.y * q.y + q.z * q.z) + l = sqrt_float(q.x * q.x + q.y * q.y + q.z * q.z) return 2.0 * atan2_approx(l, q.w) } @@ -581,8 +592,8 @@ get_twist_angle(q: Quat) -> float { // Swing angle for a cone limit, in [0, pi]; polarity cannot matter since // every term is squared. get_swing_angle(q: Quat) -> float { - x = sqrt(q.z * q.z + q.w * q.w) - y = sqrt(q.x * q.x + q.y * q.y) + x = sqrt_float(q.z * q.z + q.w * q.w) + y = sqrt_float(q.x * q.x + q.y * q.y) return 2.0 * atan2_approx(y, x) } diff --git a/aephysics/sphere/module.ae b/aephysics/sphere/module.ae index 767cbeb..9e51f1c 100644 --- a/aephysics/sphere/module.ae +++ b/aephysics/sphere/module.ae @@ -25,7 +25,6 @@ exports ( sphere_proxy, capsule_proxy, segment_proxy, length_scratch ) -extern sqrt(x: float) -> float // Closer than this is an overlap (a tenth of the linear slop). const OVERLAP_SLOP = 0.0005 @@ -117,7 +116,7 @@ ray_cast_sphere(s: Sphere, origin: Vec3, translation: Vec3, max_fraction: float) c = math.mul_add(shifted, t, d) cc = math.dot(c, c) if cc > rr { return output } - h = sqrt(rr - cc) + h = math.sqrt_float(rr - cc) fraction = t - h if fraction < 0.0 || max_fraction * length < fraction { // The intersection is off the segment; an origin inside still counts. @@ -149,7 +148,7 @@ ray_cast_hollow_sphere(s: Sphere, origin: Vec3, translation: Vec3, max_fraction: r = s.radius rr = r * r if cc > rr { return output } - h = sqrt(rr - cc) + h = math.sqrt_float(rr - cc) fraction = t - h if fraction < 0.0 { fraction = t + h } if fraction < 0.0 { return output } diff --git a/bench/RESULTS.md b/bench/RESULTS.md index f3fd790..eed74b0 100644 --- a/bench/RESULTS.md +++ b/bench/RESULTS.md @@ -1336,8 +1336,11 @@ The convex pile (5,120 hulls), ms per step, the best of three alternated pairs: collide 7.65 to 7.07, step 11.32 to 10.74. Box3D in the same period: collide 5.81-5.88, step 8.94-9.07, so the collide pass went from 1.31x to 1.20x. What is left in the SAT is the scalar work around the -lanes (the candidate edges' inscribed-sphere tests, a face's bias) and the -slice indexing the arrays go through. +lanes (the candidate edges' inscribed-sphere tests, a face's bias). The +slice indexing the arrays go through costs nothing: GCC folds every bounds +test once the views inline, and the profile samples on aether_slice.h's +line are skid from the element loads, which a raw pointer pays too +(aether#3020, measured by the Aether session). The wide contact stages gather and scatter their bodies through 8x8 register transposes (Box3D's b3GatherBodies, b3ScatterBodies). The second @@ -1352,6 +1355,28 @@ a default build (no -mavx2), large pyramid, best of four: width 4 runs 1.23), the speed of the -mavx2 build. test_determinism steps every scene at both widths in one process, bit for bit. +### std.math's square root (2026-10-20, #184) + +At -O0 every `sqrt` is a library call, and MinGW's library `sqrt` is the +x87 one at the 64-bit precision a MinGW main thread runs with. It rounds +twice: 253 of 1,000,000 random inputs come out one ulp off the +instruction, so the -O0 trace left the golden one at step 10 of the wave +pile. The engine's double square root is now std.math's, which runs in +Aether's optimized runtime as the instruction. It is a call that GCC +cannot inline into the solvers, where libm's became one `sqrtsd`. The +default build, ms a step, three alternated runs each: + +| scene | libm `sqrt`, inlined | std.math `sqrt`, called | +|---|---|---| +| large pyramid | 6.00 / 6.03 / 6.34 | 6.16 / 6.00 / 6.50 | +| many pyramids | 13.55 / 13.16 / 14.93 | 14.35 / 12.67 / 14.44 | +| joint grid | 11.38 / 12.26 / 12.05 | 11.33 / 11.86 / 11.95 | + +Within the runs' own spread (the machine was in use), and the best of +each is the same or better with the call. The bits are the same. libm's +`sqrtf` stays: its 64-bit x87 result rounds to a float once, and no +float of 304 million differs from `sqrtss`. + ## Build flags (inlining) The generated C of every Aether function is a plain `static` function, diff --git a/bench/hull.ae b/bench/hull.ae index cfc5231..f9009de 100644 --- a/bench/hull.ae +++ b/bench/hull.ae @@ -11,7 +11,6 @@ extern calloc(count: int, size: int) -> ptr extern free(p: ptr) extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float clock() -> long { return os.now_monotonic_ns() } ms(ns: long) -> float { return (ns as float) / 1000000.0 } @@ -31,8 +30,8 @@ fill_sphere(points: Vec3[], count: int) { u1 = unit_random() u2 = 2.0 * math.PI * unit_random() u3 = 2.0 * math.PI * unit_random() - s1 = sqrt(1.0 - u1) - s2 = sqrt(u1) + s1 = math.sqrt_float(1.0 - u1) + s2 = math.sqrt_float(u1) points[i] = math.vec3(s1 * sin(u2), s1 * cos(u2), s2 * sin(u3)) i = i + 1 } diff --git a/bench/suite.ae b/bench/suite.ae index 058aabe..dc8da64 100644 --- a/bench/suite.ae +++ b/bench/suite.ae @@ -465,7 +465,7 @@ spinner() { // The grid starts low, where the cylinder is narrow: the chord there bounds it. spacing = 0.55 y_start = 2.0 - chord = sqrt(radius * radius - (y_center - y_start) * (y_center - y_start)) + chord = math.sqrt_float(radius * radius - (y_center - y_start) * (y_center - y_start)) x_limit = chord - 2.0 * spacing z_limit = half_depth - 2.0 * sphere_radius x = 0.0 - x_limit @@ -677,7 +677,7 @@ random_quat() -> Quat { u3 = random_float_range(0.0, 2.0 * math.PI) cs2 = math.compute_cos_sin(u2) cs3 = math.compute_cos_sin(u3) - return math.quat_vs(math.vec3(sqrt(1.0 - u1) * cs2.sine, sqrt(1.0 - u1) * cs2.cosine, sqrt(u1) * cs3.sine), sqrt(u1) * cs3.cosine) + return math.quat_vs(math.vec3(math.sqrt_float(1.0 - u1) * cs2.sine, math.sqrt_float(1.0 - u1) * cs2.cosine, math.sqrt_float(u1) * cs3.sine), math.sqrt_float(u1) * cs3.cosine) } sat_hull(kind: int) -> *HullData { diff --git a/tests/test_determinism.ae b/tests/test_determinism.ae index a8f1aa0..fd56f1d 100644 --- a/tests/test_determinism.ae +++ b/tests/test_determinism.ae @@ -33,7 +33,6 @@ import aephysics.parallel extern calloc(count: int, size: int) -> ptr extern exit(code: int) extern free(p: ptr) -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -118,8 +117,8 @@ random_unit_vector() -> Vec3 { u1 = random_float_range(0.0, 1.0) u2 = random_float_range(0.0, 2.0 * math.PI) u3 = random_float_range(0.0, 2.0 * math.PI) - sqrt1_minus_u1 = sqrt(1.0 - u1) - sqrt_u1 = sqrt(u1) + sqrt1_minus_u1 = math.sqrt_float(1.0 - u1) + sqrt_u1 = math.sqrt_float(u1) cs2 = math.compute_cos_sin(u2) cs3 = math.compute_cos_sin(u3) return math.vec3(sqrt1_minus_u1 * cs2.sine, sqrt1_minus_u1 * cs2.cosine, sqrt_u1 * cs3.sine) @@ -129,8 +128,8 @@ random_quat() -> Quat { u1 = random_float_range(0.0, 1.0) u2 = random_float_range(0.0, 2.0 * math.PI) u3 = random_float_range(0.0, 2.0 * math.PI) - sqrt1_minus_u1 = sqrt(1.0 - u1) - sqrt_u1 = sqrt(u1) + sqrt1_minus_u1 = math.sqrt_float(1.0 - u1) + sqrt_u1 = math.sqrt_float(u1) cs2 = math.compute_cos_sin(u2) cs3 = math.compute_cos_sin(u3) return math.quat_vs(math.vec3(sqrt1_minus_u1 * cs2.sine, sqrt1_minus_u1 * cs2.cosine, sqrt_u1 * cs3.sine), sqrt_u1 * cs3.cosine) diff --git a/tests/test_distance.ae b/tests/test_distance.ae index a7741fe..2ed039e 100644 --- a/tests/test_distance.ae +++ b/tests/test_distance.ae @@ -17,7 +17,6 @@ extern exit(code: int) extern free(p: ptr) extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -136,7 +135,7 @@ test_rounded() { input.transform.p = math.vec3(3.0, 3.0, 3.0) cache = distance.empty_cache() o = distance.shape_distance(&input, &cache, null, 0) - small("sphere-box: corner distance", o.distance - (sqrt(12.0) - 0.5), 0.00001) + small("sphere-box: corner distance", o.distance - (math.sqrt_float(12.0) - 0.5), 0.00001) small("sphere-box: corner witness", o.point_a.x + o.point_a.y + o.point_a.z - 3.0, 0.00001) // Overlapping: distance zero, no normal. diff --git a/tests/test_height_field.ae b/tests/test_height_field.ae index a98db08..6b6e034 100644 --- a/tests/test_height_field.ae +++ b/tests/test_height_field.ae @@ -23,7 +23,6 @@ extern calloc(count: int, size: int) -> ptr extern exit(code: int) extern free(p: ptr) extern memcmp(a: ptr, b: ptr, size: int) -> int -extern sqrt(x: float) -> float var failures = 0 var checks = 0 diff --git a/tests/test_hull.ae b/tests/test_hull.ae index 8b9ba10..055510e 100644 --- a/tests/test_hull.ae +++ b/tests/test_hull.ae @@ -18,7 +18,6 @@ extern free(p: ptr) extern memcmp(a: ptr, b: ptr, size: int) -> int extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -119,8 +118,8 @@ test_tetrahedron() { ensure("tetrahedron: 4 faces", h.face_count == 4) euler("tetrahedron", h) expected_volume = 1.0 / 6.0 - expected_area = 1.5 + 0.5 * sqrt(3.0) - expected_inner = 0.25 / sqrt(3.0) + expected_area = 1.5 + 0.5 * math.sqrt_float(3.0) + expected_inner = 0.25 / math.sqrt_float(3.0) small("tetrahedron: volume", h.volume - expected_volume, 0.00001) small("tetrahedron: surface area", h.surface_area - expected_area, 0.00001) small("tetrahedron: inner radius", h.inner_radius - expected_inner, 0.00001) @@ -335,8 +334,8 @@ fill_sphere_sample(points: Vec3[], count: int, seed: int) { u1 = unit_random() u2 = 2.0 * math.PI * unit_random() u3 = 2.0 * math.PI * unit_random() - s1 = sqrt(1.0 - u1) - s2 = sqrt(u1) + s1 = math.sqrt_float(1.0 - u1) + s2 = math.sqrt_float(u1) points[i] = math.vec3(s1 * sin(u2), s1 * cos(u2), s2 * sin(u3)) i = i + 1 } diff --git a/tests/test_manifold.ae b/tests/test_manifold.ae index c2bf58a..6c8d638 100644 --- a/tests/test_manifold.ae +++ b/tests/test_manifold.ae @@ -21,7 +21,6 @@ extern exit(code: int) extern free(p: ptr) extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 diff --git a/tests/test_mesh.ae b/tests/test_mesh.ae index 4b18bc1..2dae29c 100644 --- a/tests/test_mesh.ae +++ b/tests/test_mesh.ae @@ -22,7 +22,6 @@ extern exit(code: int) extern free(p: ptr) extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 diff --git a/tests/test_physics_world.ae b/tests/test_physics_world.ae index 65ec352..eb8f1ad 100644 --- a/tests/test_physics_world.ae +++ b/tests/test_physics_world.ae @@ -32,7 +32,6 @@ extern calloc(count: int, size: int) -> ptr extern free(p: ptr) extern cos(x: float) -> float extern sin(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -828,9 +827,9 @@ capture_friction(friction_a: float, id_a: long, friction_b: float, id_b: long) - if id_a == (101 as long) || id_b == (101 as long) { g_saw_child0 = true } if id_a == (202 as long) || id_b == (202 as long) { g_saw_child1 = true - g_mixed_friction = sqrt(friction_a * friction_b) + g_mixed_friction = math.sqrt_float(friction_a * friction_b) } - return sqrt(friction_a * friction_b) + return math.sqrt_float(friction_a * friction_b) } // TestCompoundContactMaterials: a unit-friction sphere struck down onto diff --git a/tests/test_restitution.ae b/tests/test_restitution.ae index b0a3cbc..d568bd3 100644 --- a/tests/test_restitution.ae +++ b/tests/test_restitution.ae @@ -20,7 +20,6 @@ import aephysics.physics_world extern calloc(count: int, size: int) -> ptr extern exit(code: int) extern free(p: ptr) -extern sqrt(x: float) -> float const TIME_STEP = 0.016666666666666666 const SUB_STEP_COUNT = 4 @@ -713,7 +712,7 @@ test_impulse() { weight_impulse = mass * 10.0 * TIME_STEP j = 0 while j < 4 { - impact_speed = sqrt(20.0 * heights[j]) + impact_speed = math.sqrt_float(20.0 * heights[j]) expect_toi = impact_speed > fast_speed tolerance = 0.00001 * mass * (10.0 + impact_speed) k = 0 diff --git a/tests/test_shape.ae b/tests/test_shape.ae index 38e580a..a34e304 100644 --- a/tests/test_shape.ae +++ b/tests/test_shape.ae @@ -29,7 +29,6 @@ extern exit(code: int) extern free(p: ptr) extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -476,7 +475,7 @@ sphere_hit_error(s: Sphere, origin: Vec3, translation: Vec3, point: Vec3) -> flo d = math.normalize(translation) b = math.dot(sv, d) cc = math.dot(sv, sv) - r * r - t = 0.0 - b - sqrt(b * b - cc) + t = 0.0 - b - math.sqrt_float(b * b - cc) return math.distance(point, math.mul_add(origin, t, d)) } @@ -489,7 +488,7 @@ capsule_hit_error(c: Capsule, origin: Vec3, translation: Vec3, point: Vec3) -> f qa = math.dot(dp, dp) qb = 2.0 * math.dot(sp, dp) qc = math.dot(sp, sp) - c.radius * c.radius - tau = (0.0 - qb - sqrt(qb * qb - 4.0 * qa * qc)) / (2.0 * qa) + tau = (0.0 - qb - math.sqrt_float(qb * qb - 4.0 * qa * qc)) / (2.0 * qa) return math.distance(point, math.mul_add(origin, tau, d)) } @@ -784,7 +783,7 @@ test_dispatch() { small("mesh centroid", math.length(shape.get_shape_centroid(&sh)), 0.000001) extent = shape.compute_shape_extent(&sh, math.vec3_zero()) small("mesh max extent", extent.max_extent.x - 2.0, 0.000001) - small("mesh min extent", extent.min_extent - sqrt(8.0), 0.000001) + small("mesh min extent", extent.min_extent - math.sqrt_float(8.0), 0.000001) extent = shape.compute_shape_extent(&sh, math.vec3(1.0, 0.0, 0.0)) small("mesh max extent off centre", extent.max_extent.x - 3.0, 0.000001) ensure("mesh has no mass", shape.compute_shape_mass(&sh).mass == 0.0) diff --git a/tests/test_soak.ae b/tests/test_soak.ae index fb6a1cd..5d29b5e 100644 --- a/tests/test_soak.ae +++ b/tests/test_soak.ae @@ -26,7 +26,6 @@ import aephysics.dynamics import aephysics.physics_world extern exit(code: int) -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -255,7 +254,7 @@ soak(seed: int) { start = energy_of(bodies) // The fastest any body may go: thrown at up to 2 m/s per axis from the // top layer, it falls the whole height; spin adds its rim's speed. - speed_cap = 1.25 * sqrt(12.0 + 2.0 * GRAVITY * (top + 1.0)) + 6.0 + speed_cap = 1.25 * math.sqrt_float(12.0 + 2.0 * GRAVITY * (top + 1.0)) + 6.0 peak = start fastest = 0.0 lowest = 0.0 diff --git a/tests/test_triangle_manifold.ae b/tests/test_triangle_manifold.ae index 1d9fead..a3b9ea6 100644 --- a/tests/test_triangle_manifold.ae +++ b/tests/test_triangle_manifold.ae @@ -21,7 +21,6 @@ extern exit(code: int) extern free(p: ptr) extern sin(x: float) -> float extern cos(x: float) -> float -extern sqrt(x: float) -> float var failures = 0 var checks = 0 @@ -369,7 +368,7 @@ test_sphere_and_closest() { triangle_manifold.collide_triangle_and_sphere(&m, 8, v1, v2, v3, s) ensure("sphere edge: one point", m.point_count == 1) ensure("sphere edge: feature", m.feature == triangle_manifold.FEATURE_EDGE3) - if m.point_count == 1 { small("sphere edge: separation", manifold.point(&m, 0).separation - (sqrt(0.045) - 0.25), 0.00001) } + if m.point_count == 1 { small("sphere edge: separation", manifold.point(&m, 0).separation - (math.sqrt_float(0.045) - 0.25), 0.00001) } // Below the plane: culled. Far above: nothing. s = manifold.sphere(math.vec3(0.0, 0.0 - 0.1, 0.0), 0.25) m = manifold.local_manifold(points) From 8c34ca68d766c5eacf5b13c7a1d8b031807dafaa Mon Sep 17 00:00:00 2001 From: Nicolas Maman Date: Tue, 20 Oct 2026 14:26:42 -0300 Subject: [PATCH 4/6] The Windows leg installs diffutils for the golden trace's diff (#184) MSYS2's base install has no diff, and the first run of the leg failed there after every test had passed. --- .github/workflows/test.yaml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 264c9b1..cad67d7 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -81,12 +81,14 @@ jobs: uses: actions/checkout@v4 # UCRT64: the release's libaether.a links the Universal CRT, and does - # not link against an msvcrt gcc. + # not link against an msvcrt gcc. diffutils for the golden trace's + # diff, which MSYS2's base install lacks. - name: Set up MSYS2 uses: msys2/setup-msys2@v2 with: msystem: UCRT64 install: >- + diffutils unzip mingw-w64-ucrt-x86_64-gcc From 91d9d853057ee4c0a05073780bcdfed854a65bb3 Mon Sep 17 00:00:00 2001 From: Nicolas Maman Date: Tue, 20 Oct 2026 14:57:12 -0300 Subject: [PATCH 5/6] LF checkouts for the golden trace and the scripts (#184) The Windows runner's git checks text out as CRLF, so the golden trace differed from the build's on every line (1,1250c1,1250) although the bits agree: a UCRT64 build here matches it bit for bit. .gitattributes keeps tests/golden/*.txt and the shell scripts LF on every checkout. --- .gitattributes | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 .gitattributes diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..f8deca0 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,6 @@ +# LF on every checkout, Windows' included, where git turns text into CRLF +# by default. scripts/golden_trace.sh compares the golden trace line by +# line with what the tests print, and every line differed on the Windows +# runner (#184); the scripts run under bash. +tests/golden/*.txt text eol=lf +*.sh text eol=lf From 8f44cc5baaa5734bc17a2e59769445f2a93da8b2 Mon Sep 17 00:00:00 2001 From: Nicolas Maman Date: Tue, 20 Oct 2026 15:16:10 -0300 Subject: [PATCH 6/6] The Windows leg's comment says what it can and cannot catch (#184) ae build passes -fstack-protector-strong, and with it GCC on Win64 does not crash on the AVX2 copies at -O0 (measured on the release's C for test_contact_solver_wide: plain -O0 crashes, adding that flag alone compiles). So the leg runs -O0 on Aether 0.835.0 and cannot see aether#3023, which shows in builds that call gcc themselves without it, as ae3d's do. With std.math's sqrt nothing in a step depends on the x87 precision, so the -O0 scenes keep their bits without aether#3027. --- .github/workflows/test.yaml | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index cad67d7..51d50f6 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -66,9 +66,12 @@ jobs: # Windows, where ae3d and Siltreach run, and -O0, which ae3d builds its # suites at. Neither was built here, and two bugs reached ae3d first - # (#184): GCC on Win64 crashing on the AVX2 copies at -O0 - # (aether#3023), and four workers parting at -O0 because MinGW's - # library sqrt rounds twice on the main thread's x87 setting. + # (#184): four workers parting at -O0, because MinGW's library sqrt + # rounds twice at the main thread's x87 precision (the engine now takes + # std.math's), and GCC on Win64 crashing on the AVX2 copies at -O0 + # (aether#3023). This leg cannot see the crash: `ae build` passes + # -fstack-protector-strong, which keeps GCC clear of it, and a build + # that calls gcc itself without it (ae3d's) is where it shows. windows: name: aephysics tests (Windows x86-64) runs-on: windows-latest @@ -110,9 +113,9 @@ jobs: export PATH="$HOME/aether/bin:$PATH" scripts/test.sh - # Unoptimized, every sqrt is a library call and every function holds - # its lanes on the stack: the scenes keep their bits over four - # workers, and the trace matches the golden one. + # Unoptimized, every maths function is a library call and every + # function holds its lanes on the stack: the scenes keep their bits + # over four workers, and the trace matches the golden one. - name: The determinism scenes at -O0 env: AETHER_AEPHYSICS_CFLAGS: -O0