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 diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 768f770..51d50f6 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,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 diff --git a/README.md b/README.md index bb235a1..c4bdf11 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 @@ -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. 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 99a28e6..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) @@ -651,7 +650,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 +661,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) 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)