Windows and -O0 in CI; std.math's sqrt; SIMD widths compared where the AVX2 copy is plain

.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

.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,69 @@ 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): 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
+    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. 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
+
+      # 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 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
+        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

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
@@ -239,12 +239,14 @@ 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
-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
+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). 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.

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 }

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)

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

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))

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,15 +1797,15 @@ 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
     delta = math.compute_cos_sin(theta)
     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

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

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)