diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index cecbcb5..a4366e4 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -21,7 +21,7 @@ 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.794.0 AEB_REF=v0.325 sh
+            | AE_PIN=0.795.0 AEB_REF=v0.325 sh
           echo "$HOME/.local/bin" >> "$GITHUB_PATH"
 
       - name: Run the tests with ae
@@ -54,7 +54,7 @@ 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.794.0 AEB_REF=v0.325 sh
+            | AE_PIN=0.795.0 AEB_REF=v0.325 sh
           echo "$HOME/.local/bin" >> "$GITHUB_PATH"
 
       - name: Run the tests with ae
diff --git a/README.md b/README.md
index 7f4aa8f..ac3a83c 100644
--- a/README.md
+++ b/README.md
@@ -50,7 +50,7 @@ so a test written against the reference reads the same here.
 | `aephysics.joint_solver` | the seven joints solved: each kind's prepare (frames relative to the centres of mass, effective masses, spring softness), warm start and solve (rigid or soft constraints, speculative limits, springs, motors), the kinds' accessors (limits, springs, motors, current angles and translations, forces and torques), joint.c's dispatch with the constraint hertz clamp; beyond the reference, a spherical joint's cone may be centred off frame A's z and may be an ellipse of two half-angles (the circle stays the default, bit for bit) | `test_joint_solver.ae` (181 checks: test_joint.c's accessors on every kind, each kind solved by hand) | [ours alone](bench/RESULTS.md#joint_solver) |
 | `aephysics.solver` | the Soft Step over the workers, the reference's stages: the work in blocks (the bodies', each colour's joints, wide and scalar constraints) that workers claim by compare-and-swap, the stages published as sync bits by worker 0 -- the constraints prepared, per sub-step the velocities integrated (gravity, damping, the gyroscopic torque), warm start, solve, positions, relax, colour by colour with the overflow first on worker 0; restitution, the impulses stored; the bodies finalised in parallel (sleep velocities, move events, fast bodies swept for the time of impact, bounds and proxies), the joint and hit events, the trees refit, bullets, the sensors' hits, islands put to sleep | `test_solver.ae` (45 checks: free fall against the closed form, resting and sleeping, a stack, a bounce and its hit event, a joint event, a pendulum, a fast sphere stopped by a thin wall, a bullet, a sensor swept, two worlds bit for bit alike) | [1.2-1.3x the reference's step, stage by stage](bench/RESULTS.md#stage-by-stage) |
 | `aephysics.physics_world` | the world's face: the step (pairs, narrow phase, solve, sensors, events), the events read back, the settings and counters, queries over every shape (overlap of a box or a proxy, the mover's planes, ray, shape and mover casts, the closest ray hit) and against one body at a transform of the caller's, explosions | `test_physics_world.ae` (238 checks: test_world.c whole but the hull database, among it empty and recycled worlds, every setting, the worker count set between steps, contact, hit, move and sensor events, a compound child's hit events and materials, the overflow colour pile, the explosion near and far, the bullet and sleep flags, the compound count, an enlarged proxy destroyed; test_body_query.c whole: casts, overlaps of hull proxies and turned compounds, mover planes, and the mover's time of impact against 2,000 random sweeps; test_broad_phase.c's static rebuild and its stress test, a swarm through every proxy change with the pair set checked by brute force every step; test_id.c; the world queries; a wave pile stepping alike twice) | [the reference's own benchmark scenes with the same checksums: 1.2x its AVX2 build on the large pyramid, 1.15x on the many pyramids, 1.05x on the joint grid, ours at eight lanes in AVX2 with float records (#121)](bench/RESULTS.md#the-collide-pass-as-the-reference-keeps-it-2026-10-09-121); [the rest of its suite, nine scenes and the SAT runs, at 1.04-1.21x](bench/RESULTS.md#suite-the-references-benchmark-suite) |
-| `aephysics.human` | the ragdoll: twelve capsule bones on spherical joints with cone and twist limits and revolute joints with angle limits -- a person's ranges by default (shoulders that reach overhead on a cone centred off the rest, hips on an elliptical cone through 125 degrees ahead, 25 back, 45 out and 30 across, knees and elbows that never pass straight), the reference's tight ones on request -- soles that grip, a spring on every joint toward the reference pose, a motor whose torque limit is joint friction, filter joints for the limbs that clash; the align spring, kinematic anchors through motor or parallel joints and their targets driven over a step (the pose drive of an active ragdoll), velocity, kicks, bullets | `test_human.ae` (96 checks: the figure's shape, a fall to rest in one piece, the setters, a kick, standing under the align spring, the pose held on anchors, the anchors driven two metres and the figure walked behind them, the same drop twice bit for bit, a fitted shape, each arm reaching flexion 170, abduction 160 and extension 45 degrees standing and prone, the hips' cone through a person's four ranges and short of a split, the hinges stopping at straight and folding to their range, the soles' grip) | [the reference's rain benchmark, 300 ragdolls over 400 steps, at 0.9x](bench/RESULTS.md#human) |
+| `aephysics.human` | the ragdoll: twelve capsule bones on spherical joints with cone and twist limits and revolute joints with angle limits -- a person's ranges by default (shoulders that reach overhead on a cone centred off the rest, hips on an elliptical cone through 125 degrees ahead, 25 back, 45 out and 30 across, knees and elbows that never pass straight), the reference's tight ones on request -- soles that grip, and on request (`Human.ankles`, past the reference) the feet as bones of their own on ankles, an elliptical cone through 50 degrees down, 35 up and 15 aside with 15 of roll, so a crouched figure's soles lie flat, a spring on every joint toward the reference pose, a motor whose torque limit is joint friction, filter joints for the limbs that clash; the align spring, kinematic anchors through motor or parallel joints and their targets driven over a step (the pose drive of an active ragdoll), velocity, kicks, bullets | `test_human.ae` (115 checks: the figure's shape, a fall to rest in one piece, the setters, a kick, standing under the align spring, the pose held on anchors, the anchors driven two metres and the figure walked behind them, the same drop twice bit for bit, a fitted shape, each arm reaching flexion 170, abduction 160 and extension 45 degrees standing and prone, the hips' cone through a person's four ranges and short of a split, the hinges stopping at straight and folding to their range, the soles' grip, the ankles' four ranges measured through the joint and the soles flat under shins tipped 25 degrees over them) | [the reference's rain benchmark, 300 ragdolls over 400 steps, at 0.9x](bench/RESULTS.md#human) |
 | `aephysics` | the public API | |
 
 Deliberate choices:
@@ -77,11 +77,11 @@ Install the released toolchain with the shared, version-pinned installer
 
 ```sh
 curl -fsSL https://raw.githubusercontent.com/aether-lang-dev/aeb/main/get.sh \
-  | AE_PIN=0.794.0 AEB_REF=v0.325 sh
+  | AE_PIN=0.795.0 AEB_REF=v0.325 sh
 export PATH="$HOME/.local/bin:$PATH"
 ```
 
-This installs `ae` 0.794.0, the release CI pins, and `aeb` v0.325 under
+This installs `ae` 0.795.0, the release CI pins, and `aeb` v0.325 under
 `~/.local` (`PREFIX` overrides it), using release binaries where available.
 The code needs 0.788.0 or newer (the eight-wide lanes, `f32x8`).
 
@@ -166,7 +166,17 @@ 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
 cases, a hull manifold and a shape's fat bounds at 1e7 against the origin;
-56 checks) and `test_edge_cases.ae` (every
+56 checks), `test_restitution.ae` (the reference's test_restitution.c
+whole: the coefficient head on at three mass ratios and at sixteen impact
+phases, a ball on a column of dead balls with propagation off and on, flat
+and spinning landings, elastic drops from 40, 20 and 10 m, a cube dropped
+flat four times, the threshold, a resting stack that stays put, the
+published approach speed, the narrow-floor overshoot, the energy of
+thirteen elastic scenes, every step's impulse balanced on drops either side
+of continuous collision, and one, two and four workers hashing alike; its
+numbers are Box3D's own run's to every printed digit, the few that differ
+in the last digit differing between the reference's single and double
+builds too; 162 checks) and `test_edge_cases.ae` (every
 number a world takes refused when it is NaN, infinite or negative where
 it must not be, as the reference asserts it, the world left as it was;
 a world or a ragdoll of such numbers never made, a ragdoll made whole
diff --git a/aephysics/dynamics/module.ae b/aephysics/dynamics/module.ae
index 0ab436f..5b59207 100644
--- a/aephysics/dynamics/module.ae
+++ b/aephysics/dynamics/module.ae
@@ -68,7 +68,7 @@ exports (
     default_world_def, default_body_def, default_shape_def, null_body_id, null_shape_id,
     create_world, destroy_world, world_is_valid, get_world, world_body_count, world_shape_count, world_awake_body_count,
     world_set_gravity, world_get_gravity,
-    create_body, destroy_body, body_is_valid, get_body_full_id, get_body, get_body_sim, get_body_state, get_body_transform,
+    create_body, destroy_body, body_is_valid, get_body_full_id, get_body, get_body_sim, get_body_state, get_body_transform, hash_world_state,
     make_body_id, sync_body_flags, update_body_mass_data, is_body_awake, wake_body,
     body_get_position, body_get_rotation, body_get_transform, body_set_transform,
     body_get_local_point, body_get_world_point, body_get_local_vector, body_get_world_vector,
@@ -1750,6 +1750,43 @@ get_body_state(w: *World, body: *Body) -> *BodyState {
 }
 
 get_body_transform(w: *World, body_id: int) -> Transform { return get_body_sim(w, body_at(w, body_id)).transform }
+
+// The reference's b3HashWorldState: every live body's pose and, awake, its
+// velocities, folded as their bits by 64-bit FNV-1a in body order. Two runs
+// alike to the bit hash alike (worker parity, replay). The pose is held in
+// double here, so its bits are the double's, as the reference's large world
+// build mixes its positions. Aether's integers wrap (the language
+// reference), so the multiply is the FNV's own.
+hash_world_state(w: *World) -> long {
+    // 14695981039346656037 and 1099511628211, the FNV-1a basis and prime.
+    hash = 0 - 3750763034362895579
+    prime = 1099511628211 as long
+    i = 0
+    while i < w.bodies.count {
+        body = body_at(w, i)
+        if body.id == i {
+            xf = get_body_sim(w, body).transform
+            hash = (hash ^ native.f64_bits(xf.p.x)) * prime
+            hash = (hash ^ native.f64_bits(xf.p.y)) * prime
+            hash = (hash ^ native.f64_bits(xf.p.z)) * prime
+            hash = (hash ^ native.f64_bits(xf.q.x)) * prime
+            hash = (hash ^ native.f64_bits(xf.q.y)) * prime
+            hash = (hash ^ native.f64_bits(xf.q.z)) * prime
+            hash = (hash ^ native.f64_bits(xf.q.w)) * prime
+            state = get_body_state(w, body)
+            if state != null {
+                hash = (hash ^ native.f32_bits(state.v.x)) * prime
+                hash = (hash ^ native.f32_bits(state.v.y)) * prime
+                hash = (hash ^ native.f32_bits(state.v.z)) * prime
+                hash = (hash ^ native.f32_bits(state.w.x)) * prime
+                hash = (hash ^ native.f32_bits(state.w.y)) * prime
+                hash = (hash ^ native.f32_bits(state.w.z)) * prime
+            }
+        }
+        i = i + 1
+    }
+    return hash
+}
 make_body_id(w: *World, body_id: int) -> BodyId { return BodyId { index1: body_id + 1, world0: w.world_id, generation: body_at(w, body_id).generation } }
 
 // The body's flags onto its sim and state, never the transient ones,
diff --git a/aephysics/human/module.ae b/aephysics/human/module.ae
index b843fea..301c9e9 100644
--- a/aephysics/human/module.ae
+++ b/aephysics/human/module.ae
@@ -34,6 +34,7 @@ import aephysics.joint_solver
 exports (
     BONE_PELVIS, BONE_SPINE_01, BONE_SPINE_03, BONE_NECK, BONE_THIGH_L, BONE_CALF_L, BONE_THIGH_R, BONE_CALF_R,
     BONE_UPPER_ARM_L, BONE_LOWER_ARM_L, BONE_UPPER_ARM_R, BONE_LOWER_ARM_R, BONE_COUNT, FILTER_JOINT_COUNT,
+    BONE_FOOT_L, BONE_FOOT_R, MAX_BONE_COUNT, ANKLE_PLANTAR_DEGREES, ANKLE_DORSAL_DEGREES, ANKLE_SIDE_DEGREES, ANKLE_ROLL_DEGREES,
     HumanBone, Human, BoneSpec, bone_at, filter_joint_at, spec_at,
     human_reference_shape, human_anatomical_limits, human_tight_limits, human_chain_child, human_shape_head, human_shape_tail, human_shape_toe, human_fit_shape, create_human_shaped,
     create_human, destroy_human, human_set_velocity, human_apply_angular_impulse, human_set_joint_friction_torque,
@@ -42,6 +43,8 @@ exports (
 )
 
 extern calloc(count: int, size: int) -> ptr
+extern sin(x: float) -> float
+extern cos(x: float) -> float
 extern free(p: ptr)
 
 const BONE_PELVIS = 0
@@ -58,6 +61,20 @@ const BONE_UPPER_ARM_R = 10
 const BONE_LOWER_ARM_R = 11
 const BONE_COUNT = 12
 const FILTER_JOINT_COUNT = 3
+// With ankles (#152, past the reference, whose figure has no foot bodies)
+// the feet are bones of their own after the twelve, each on a joint below
+// its calf; without, the reference's twelve, the foot a shape on the calf.
+const BONE_FOOT_L = 12
+const BONE_FOOT_R = 13
+const MAX_BONE_COUNT = 14
+// The ankle's range, a loaded one's: the toes down 50 degrees and up 35 (an
+// elliptical cone about a centre tipped toward the toes), 15 to either
+// side, and the sole rolled 15 either way about the foot's line.
+const ANKLE_PLANTAR_DEGREES = 50.0
+const ANKLE_DORSAL_DEGREES = 35.0
+const ANKLE_SIDE_DEGREES = 15.0
+const ANKLE_ROLL_DEGREES = 15.0
+const ANKLE_FRICTION = 0.6
 
 // One bone: its body, the joint to its parent, the anchor a pose drive
 // hangs it from, the joint's frames on parent and bone, the bone's
@@ -85,13 +102,18 @@ struct HumanBone {
     cone_major_axis: Vec3          // in local_frame_a
 }
 
-// The ragdoll: a heap block of bones and filter joints (create allocates, destroy frees).
+// The ragdoll: a heap block of bones and filter joints (create allocates,
+// destroy frees). `ankles`, set before it is made, gives the feet bones of
+// their own (bone_count MAX_BONE_COUNT); left false it is the reference's
+// figure of BONE_COUNT bones.
 struct Human {
-    bones: ptr                     // HumanBone[BONE_COUNT]
+    bones: ptr                     // HumanBone[bone_count]
     filter_joints: ptr             // JointId[FILTER_JOINT_COUNT]
     filter_joint_count: int
     friction_torque: float
     is_spawned: bool
+    ankles: bool
+    bone_count: int
 }
 
 bone_at(h: *Human, index: int) -> *HumanBone { return (h.bones + index * sizeof(HumanBone)) as *HumanBone }
@@ -606,7 +628,7 @@ point_foot(spec: *BoneSpec, axis: Vec3, toe: Vec3, head: Vec3) {
 
 // A bone's body and shapes: the body at its frame placed by `placement`,
 // the capsule (and a calf's foot) with its material; the joint data on the bone.
-place_bone(h: *Human, world_id: WorldId, index: int, placement: Transform, spec: *BoneSpec, user_data: ptr, group_index: int) {
+place_bone(h: *Human, world_id: WorldId, index: int, placement: Transform, spec: *BoneSpec, user_data: ptr, group_index: int, foot_on_calf: bool) {
     bone = bone_at(h, index)
     bone.parent_index = spec.parent
     bone.reference_frame = spec.reference
@@ -623,7 +645,7 @@ place_bone(h: *Human, world_id: WorldId, index: int, placement: Transform, spec:
     sdef.base_material.rolling_resistance = 0.1
     bone.shape_id = dynamics.create_capsule_shape(bone.body_id, &sdef, spec.capsule)
     bone.foot_shape_id = dynamics.null_shape_id()
-    if spec.has_foot {
+    if spec.has_foot && foot_on_calf {
         fdef = dynamics.default_shape_def()
         fdef.base_material.friction = spec.foot_friction
         fdef.base_material.rolling_resistance = 0.1
@@ -640,6 +662,50 @@ place_bone(h: *Human, world_id: WorldId, index: int, placement: Transform, spec:
     bone.joint_friction = spec.joint_friction
 }
 
+// A calf's foot as a bone of its own: a body in the calf's frame carrying
+// the foot's capsule, on a spherical joint at the ankle (the far end of the
+// calf's capsule). The joint's frame looks down the foot (z from the ankle
+// to the toe) with x the way the toes go down (the shin's line across the
+// foot), so its elliptical cone reaches ANKLE_PLANTAR_DEGREES toward x and
+// ANKLE_DORSAL_DEGREES back (a centre tipped toward x by half their
+// difference), ANKLE_SIDE_DEGREES across, and its twist is the sole's roll.
+// The foot rests in the calf's frame, so both frames are one at rest.
+place_foot(h: *Human, world_id: WorldId, index: int, calf_index: int, placement: Transform, spec: *BoneSpec, user_data: ptr, group_index: int) {
+    bone = bone_at(h, index)
+    bone.parent_index = calf_index
+    bone.reference_frame = spec.reference
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.user_data = user_data
+    bdef.rotation = math.mul_quat(placement.q, spec.reference.q)
+    bdef.position = math.transform_point(placement, spec.reference.p)
+    bone.body_id = dynamics.create_body(world_id, &bdef)
+    fdef = dynamics.default_shape_def()
+    fdef.density = fdef.density * spec.density
+    fdef.filter.group_index = spec.group * group_index
+    fdef.base_material.friction = spec.foot_friction
+    fdef.base_material.rolling_resistance = 0.1
+    bone.shape_id = dynamics.create_capsule_shape(bone.body_id, &fdef, spec.foot)
+    bone.foot_shape_id = bone.shape_id
+    ankle = spec.capsule.center1
+    if math.length(spec.capsule.center2) > math.length(ankle) { ankle = spec.capsule.center2 }
+    z = math.normalize(math.sub(foot_toe(spec), ankle))
+    shin = math.normalize(ankle)
+    x = math.normalize(math.sub(shin, math.mul_sv(math.dot(shin, z), z)))
+    y = math.cross(z, x)
+    frame = Transform { p: ankle, q: math.make_quat_from_matrix(Matrix3 { cx: x, cy: y, cz: z }) }
+    bone.joint_kind = dynamics.JOINT_SPHERICAL
+    bone.local_frame_a = frame
+    bone.local_frame_b = frame
+    tilt = 0.5 * (ANKLE_PLANTAR_DEGREES - ANKLE_DORSAL_DEGREES) * math.DEG_TO_RAD
+    bone.cone_axis = math.vec3(sin(tilt), 0.0, cos(tilt))
+    bone.swing_limit = 0.5 * (ANKLE_PLANTAR_DEGREES + ANKLE_DORSAL_DEGREES) * math.DEG_TO_RAD
+    bone.swing_minor_limit = ANKLE_SIDE_DEGREES * math.DEG_TO_RAD
+    bone.cone_major_axis = math.vec3_zero()
+    bone.twist_limit = Vec2 { x: 0.0 - ANKLE_ROLL_DEGREES * math.DEG_TO_RAD, y: ANKLE_ROLL_DEGREES * math.DEG_TO_RAD }
+    bone.joint_friction = ANKLE_FRICTION
+}
+
 // The reference's figure standing at `position`.
 create_human(h: *Human, world_id: WorldId, position: Vec3, friction_torque: float, hertz: float, damping_ratio: float, group_index: int, user_data: ptr) {
     shape = human_reference_shape()
@@ -655,11 +721,15 @@ create_human_shaped(h: *Human, world_id: WorldId, placement: Transform, shape: p
     if h.is_spawned { return }
     if dynamics.world_is_valid(world_id) == false || shape == null || math.is_valid_world_transform(placement) == false { return }
     if valid_nonnegative(friction_torque) == false || valid_nonnegative(hertz) == false || valid_nonnegative(damping_ratio) == false { return }
-    h.bones = calloc(BONE_COUNT, sizeof(HumanBone))
+    // Ankles need a foot on each calf to stand on.
+    if h.ankles && (spec_at(shape, BONE_CALF_L).has_foot == false || spec_at(shape, BONE_CALF_R).has_foot == false) { return }
+    h.bone_count = BONE_COUNT
+    if h.ankles { h.bone_count = MAX_BONE_COUNT }
+    h.bones = calloc(h.bone_count, sizeof(HumanBone))
     h.filter_joints = calloc(FILTER_JOINT_COUNT, sizeof(JointId))
     h.filter_joint_count = 0
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         bone.body_id = dynamics.null_body_id()
         bone.anchor_id = dynamics.null_body_id()
@@ -683,14 +753,18 @@ create_human_shaped(h: *Human, world_id: WorldId, placement: Transform, shape: p
     h.friction_torque = friction_torque
     i = 0
     while i < BONE_COUNT {
-        place_bone(h, world_id, i, placement, spec_at(shape, i), user_data, group_index)
+        place_bone(h, world_id, i, placement, spec_at(shape, i), user_data, group_index, h.ankles == false)
         i = i + 1
     }
+    if h.ankles {
+        place_foot(h, world_id, BONE_FOOT_L, BONE_CALF_L, placement, spec_at(shape, BONE_CALF_L), user_data, group_index)
+        place_foot(h, world_id, BONE_FOOT_R, BONE_CALF_R, placement, spec_at(shape, BONE_CALF_R), user_data, group_index)
+    }
     // A bone the world refused (a full world, a shape spec out of range)
     // would leave a figure missing a limb that every later call reaches
     // for: what was made is taken down and the figure stays unspawned.
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         if dynamics.body_is_valid(bone_at(h, i).body_id) == false {
             take_down(h)
             return
@@ -701,7 +775,7 @@ create_human_shaped(h: *Human, world_id: WorldId, placement: Transform, shape: p
     // The joints: a spring toward the reference pose when asked, the
     // motor's torque limit as friction, the limits as the rig allows.
     i = 1
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         parent = bone_at(h, bone.parent_index)
         bone.local_frame_a.q = math.normalize_quat(bone.local_frame_a.q)
@@ -777,14 +851,14 @@ take_down(h: *Human) {
         i = i + 1
     }
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if dynamics.joint_is_valid(bone.joint_id) { dynamics.destroy_joint(bone.joint_id, false) }
         if dynamics.joint_is_valid(bone.anchor_joint_id) { dynamics.destroy_joint(bone.anchor_joint_id, false) }
         i = i + 1
     }
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if dynamics.body_is_valid(bone.body_id) { dynamics.destroy_body(bone.body_id) }
         if dynamics.body_is_valid(bone.anchor_id) { dynamics.destroy_body(bone.anchor_id) }
@@ -801,7 +875,7 @@ take_down(h: *Human) {
 human_set_velocity(h: *Human, velocity: Vec3) {
     if h.is_spawned == false { return }
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if dynamics.body_is_valid(bone.body_id) { dynamics.body_set_linear_velocity(bone.body_id, velocity) }
         i = i + 1
@@ -818,7 +892,7 @@ human_set_joint_friction_torque(h: *Human, torque: float) {
     if h.is_spawned == false || valid_nonnegative(torque) == false { return }
     h.friction_torque = torque
     i = 1
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if bone.joint_kind == dynamics.JOINT_REVOLUTE {
             joint_solver.revolute_joint_set_max_motor_torque(bone.joint_id, bone.joint_friction * torque)
@@ -832,7 +906,7 @@ human_set_joint_friction_torque(h: *Human, torque: float) {
 human_set_joint_spring_hertz(h: *Human, hertz: float) {
     if h.is_spawned == false { return }
     i = 1
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if bone.joint_kind == dynamics.JOINT_REVOLUTE {
             joint_solver.revolute_joint_set_spring_hertz(bone.joint_id, hertz)
@@ -846,7 +920,7 @@ human_set_joint_spring_hertz(h: *Human, hertz: float) {
 human_set_joint_damping_ratio(h: *Human, damping_ratio: float) {
     if h.is_spawned == false { return }
     i = 1
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if bone.joint_kind == dynamics.JOINT_REVOLUTE {
             joint_solver.revolute_joint_set_spring_damping_ratio(bone.joint_id, damping_ratio)
@@ -903,7 +977,7 @@ human_create_motor_anchors(h: *Human, world_id: WorldId) {
     mdef.max_spring_force = math.MAX_FLOAT
     mdef.max_spring_torque = math.MAX_FLOAT
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         xf = dynamics.body_get_transform(bone.body_id)
         adef.position = xf.p
@@ -929,7 +1003,7 @@ human_create_parallel_anchors(h: *Human, world_id: WorldId) {
     pdef.damping_ratio = 1.0
     pdef.max_torque = 800.0
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         xf = dynamics.body_get_transform(bone.body_id)
         adef.position = xf.p
@@ -958,7 +1032,7 @@ human_has_anchors(h: *Human) -> bool {
 // as the animation does, colliding as the ragdoll it is. Nothing happens
 // for a bone without an anchor.
 human_set_anchor_target(h: *Human, index: int, target: Transform, time_step: float) {
-    if h.is_spawned == false || index < 0 || index >= BONE_COUNT { return }
+    if h.is_spawned == false || index < 0 || index >= h.bone_count { return }
     bone = bone_at(h, index)
     if dynamics.body_is_valid(bone.anchor_id) == false { return }
     dynamics.body_set_target_transform(bone.anchor_id, target, time_step, true)
@@ -969,7 +1043,7 @@ human_set_anchor_target(h: *Human, index: int, target: Transform, time_step: flo
 human_destroy_anchors(h: *Human) {
     if h.is_spawned == false { return }
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         bone = bone_at(h, i)
         if dynamics.joint_is_valid(bone.anchor_joint_id) { dynamics.destroy_joint(bone.anchor_joint_id, true) }
         if dynamics.body_is_valid(bone.anchor_id) { dynamics.destroy_body(bone.anchor_id) }
@@ -982,7 +1056,7 @@ human_destroy_anchors(h: *Human) {
 human_set_bullet(h: *Human, flag: bool) {
     if h.is_spawned == false { return }
     i = 0
-    while i < BONE_COUNT {
+    while i < h.bone_count {
         dynamics.body_set_bullet(bone_at(h, i).body_id, flag)
         i = i + 1
     }
diff --git a/aephysics/native/aephysics_inline.h b/aephysics/native/aephysics_inline.h
index 50a1afe..2091b85 100644
--- a/aephysics/native/aephysics_inline.h
+++ b/aephysics/native/aephysics_inline.h
@@ -70,4 +70,20 @@ static inline int aephysics_load_u16(const void *p, long long index)
     return value;
 }
 
+// A single's and a double's bits, unsigned in the single's case (the
+// world state hash, which mixes every value as its bits).
+static inline long long aephysics_f32_bits(float value)
+{
+    unsigned int bits;
+    __builtin_memcpy(&bits, &value, 4);
+    return bits;
+}
+
+static inline long long aephysics_f64_bits(double value)
+{
+    long long bits;
+    __builtin_memcpy(&bits, &value, 8);
+    return bits;
+}
+
 #endif
diff --git a/aephysics/native/module.ae b/aephysics/native/module.ae
index d3b42f6..09b3154 100644
--- a/aephysics/native/module.ae
+++ b/aephysics/native/module.ae
@@ -11,7 +11,7 @@ exports (
     MAX_WORKERS, worker_index, set_worker_index, claim_slots, release_slots, enter_traversal, leave_traversal, yield_thread, pause, processor_count, scratch,
     atomic_load, atomic_store, atomic_add, atomic_or, atomic_cas,
     semaphore_new, semaphore_free, semaphore_wait, semaphore_signal,
-    name_copy, name_text, name_free, prefetch, load_word, load_u8, load_u16
+    name_copy, name_text, name_free, prefetch, load_word, load_u8, load_u16, f32_bits, f64_bits
 )
 
 // The most workers a world can have: every module's scratch is a block
@@ -40,6 +40,8 @@ extern aephysics_prefetch(p: ptr) @c_import
 extern aephysics_load_word(p: ptr, index: long) -> long @c_import
 extern aephysics_load_u8(p: ptr, index: long) -> int @c_import
 extern aephysics_load_u16(p: ptr, index: long) -> int @c_import
+extern aephysics_f32_bits(value: f32) -> long @c_import
+extern aephysics_f64_bits(value: float) -> long @c_import
 extern aephysics_name_copy(s: string) -> ptr
 extern aephysics_name_text(name: ptr) -> string
 extern free(p: ptr)
@@ -58,6 +60,10 @@ load_u8(p: ptr, index: long) -> int { return aephysics_load_u8(p, index) }
 // every platform aephysics runs on).
 load_u16(p: ptr, index: long) -> int { return aephysics_load_u16(p, index) }
 
+// A value's bits: a single's as an unsigned 32-bit value, a double's whole.
+f32_bits(value: f32) -> long { return aephysics_f32_bits(value) }
+f64_bits(value: float) -> long { return aephysics_f64_bits(value) }
+
 // A name the world keeps: a plain copy of `s` (null for none), its text
 // (empty for none), and its release.
 name_copy(s: string) -> ptr { return aephysics_name_copy(s) }
diff --git a/tests/test_human.ae b/tests/test_human.ae
index ff372c0..bb1438d 100644
--- a/tests/test_human.ae
+++ b/tests/test_human.ae
@@ -20,6 +20,7 @@ extern calloc(count: int, size: int) -> ptr
 extern exit(code: int)
 extern free(p: ptr)
 extern atan2(y: float, x: float) -> float
+extern asin(x: float) -> float
 extern cos(x: float) -> float
 extern sin(x: float) -> float
 
@@ -813,6 +814,218 @@ test_grip() {
     dynamics.destroy_world(world_id)
 }
 
+// #152: the feet as bones of their own. With ankles the figure has 14
+// bones, each foot a body on a spherical joint below its calf carrying the
+// foot's capsule (the calf keeps only its shin); every joint starts closed
+// and inside its limits; the ankle's ranges, measured through the joint
+// with the calf held and a steady torque on the foot, are a person's (the
+// toes down 50 degrees and up 35, 15 to the side, the sole rolled 15); the
+// figure falls and comes to rest whole; and a shape whose calf has no foot
+// is refused. The default figure stays the reference's twelve (every
+// other test here).
+ankle_frame(h: *Human, foot_index: int) -> Quat {
+    bone = human.bone_at(h, foot_index)
+    calf = human.bone_at(h, bone.parent_index)
+    return math.mul_quat(dynamics.body_get_rotation(calf.body_id), bone.local_frame_a.q)
+}
+
+// The foot's line against the ankle's frame on the calf: its swing toward
+// the frame's x (toes down) and toward its y (the side), in degrees.
+ankle_swing(h: *Human, foot_index: int) -> Vec2 {
+    bone = human.bone_at(h, foot_index)
+    frame = ankle_frame(h, foot_index)
+    foot_z = math.rotate_vector(math.mul_quat(dynamics.body_get_rotation(bone.body_id), bone.local_frame_b.q), math.vec3_axis_z())
+    along = math.dot(foot_z, math.rotate_vector(frame, math.vec3_axis_z()))
+    toward_x = math.dot(foot_z, math.rotate_vector(frame, math.vec3_axis_x()))
+    toward_y = math.dot(foot_z, math.rotate_vector(frame, math.vec3_axis_y()))
+    return Vec2 { x: atan2(toward_x, along) * math.RAD_TO_DEG, y: atan2(toward_y, along) * math.RAD_TO_DEG }
+}
+
+// The left foot pushed by a steady torque (in the ankle's frame) against
+// its limit, from the rest, the calf held still; what it settles at.
+pushed_foot(torque: Vec3) -> Vec3 {
+    wdef = dynamics.default_world_def()
+    wdef.gravity = math.vec3_zero()
+    world_id = dynamics.create_world(&wdef)
+    h = Human { bones: null, filter_joints: null, filter_joint_count: 0, friction_torque: 0.0, is_spawned: false }
+    h.ankles = true
+    human.create_human(&h, world_id, math.vec3(0.0, 1.0, 0.0), 0.0, 0.0, 0.0, 1, null)
+    calf = human.bone_at(&h, human.BONE_CALF_L)
+    dynamics.body_set_type(calf.body_id, dynamics.BODY_STATIC)
+    foot = human.bone_at(&h, human.BONE_FOOT_L)
+    frame = ankle_frame(&h, human.BONE_FOOT_L)
+    world_torque = math.rotate_vector(frame, torque)
+    i = 0
+    while i < 240 {
+        dynamics.body_apply_torque(foot.body_id, world_torque, true)
+        physics_world.step(world_id, 1.0 / 60.0, 4)
+        i = i + 1
+    }
+    swing = ankle_swing(&h, human.BONE_FOOT_L)
+    roll = joint_solver.spherical_joint_get_twist_angle(foot.joint_id) * math.RAD_TO_DEG
+    human.destroy_human(&h)
+    dynamics.destroy_world(world_id)
+    return math.vec3(swing.x, swing.y, roll)
+}
+
+test_ankles() {
+    world_id = world_with_ground()
+    h = Human { bones: null, filter_joints: null, filter_joint_count: 0, friction_torque: 0.0, is_spawned: false }
+    h.ankles = true
+    human.create_human(&h, world_id, math.vec3(0.0, 1.0, 0.0), 5.0, 1.0, 0.7, 1, null)
+    ensure("ankles: spawned with 14 bones", h.is_spawned && h.bone_count == human.MAX_BONE_COUNT)
+    counters = physics_world.get_counters(world_id)
+    ensure("ankles: 15 bodies, 16 joints (${counters.body_count}, ${counters.joint_count})", counters.body_count == 15 && counters.joint_count == 16)
+    shapes_block = calloc(4, sizeof(ShapeId))
+    shapes = shapes_block as ShapeId[]
+    side = 0
+    while side < 2 {
+        calf_index = human.BONE_CALF_L
+        foot_index = human.BONE_FOOT_L
+        if side == 1 {
+            calf_index = human.BONE_CALF_R
+            foot_index = human.BONE_FOOT_R
+        }
+        calf = human.bone_at(&h, calf_index)
+        foot = human.bone_at(&h, foot_index)
+        ensure("ankles (side ${side}): the foot hangs from its calf on a spherical joint",
+               foot.parent_index == calf_index && dynamics.joint_is_valid(foot.joint_id) && dynamics.joint_get_type(foot.joint_id) == dynamics.JOINT_SPHERICAL)
+        ensure("ankles (side ${side}): the calf keeps only its shin", dynamics.body_get_shapes(calf.body_id, shapes, 4) == 1 && calf.foot_shape_id.index1 == 0)
+        ensure("ankles (side ${side}): the foot body carries the foot's capsule", dynamics.body_get_shapes(foot.body_id, shapes, 4) == 1 &&
+               dynamics.shape_get_type(shapes[0]) == shape.SHAPE_CAPSULE && foot.foot_shape_id.index1 == foot.shape_id.index1 &&
+               math.abs_float(dynamics.shape_get_capsule(shapes[0]).radius - 0.07) < 0.000001)
+        side = side + 1
+    }
+    free(shapes_block)
+    worst_gap = 0.0
+    worst_angle = 0.0
+    i = 1
+    while i < h.bone_count {
+        joint = human.bone_at(&h, i).joint_id
+        worst_gap = math.max_float(worst_gap, dynamics.joint_get_linear_separation(joint))
+        worst_angle = math.max_float(worst_angle, dynamics.joint_get_angular_separation(joint))
+        i = i + 1
+    }
+    ensure("ankles: every joint closed and inside its limits at the rest (${worst_gap} m, ${worst_angle} rad)", worst_gap < 0.001 && worst_angle < 0.001)
+    // It falls and comes to rest whole.
+    s = 0
+    while s < 300 {
+        physics_world.step(world_id, 1.0 / 60.0, 4)
+        s = s + 1
+    }
+    whole = true
+    i = 0
+    while i < h.bone_count {
+        p = dynamics.body_get_position(human.bone_at(&h, i).body_id)
+        if math.is_valid_vec3(p) == false || p.y < 0.0 - 0.2 { whole = false }
+        i = i + 1
+    }
+    ensure("ankles: fallen and at rest, every bone finite and above the ground", whole)
+    worst_gap = 0.0
+    i = 1
+    while i < h.bone_count {
+        worst_gap = math.max_float(worst_gap, dynamics.joint_get_linear_separation(human.bone_at(&h, i).joint_id))
+        i = i + 1
+    }
+    ensure("ankles: the joints held through the fall (${worst_gap} m)", worst_gap < 0.01)
+    human.destroy_human(&h)
+    ensure("ankles: destroyed", h.is_spawned == false && physics_world.get_counters(world_id).body_count == 1)
+    dynamics.destroy_world(world_id)
+    // The ranges, each pushed to its limit through the joint.
+    toes_down = pushed_foot(math.vec3(0.0, 0.5, 0.0))
+    toes_up = pushed_foot(math.vec3(0.0, 0.0 - 0.5, 0.0))
+    aside = pushed_foot(math.vec3(0.0 - 0.5, 0.0, 0.0))
+    rolled = pushed_foot(math.vec3(0.0, 0.0, 0.5))
+    println("human: ankle ranges: toes down ${toes_down.x}, up ${toes_up.x}, aside ${aside.y}, rolled ${rolled.z} degrees")
+    ensure("ankles: the toes down to ${human.ANKLE_PLANTAR_DEGREES} (${toes_down.x})", math.abs_float(toes_down.x - human.ANKLE_PLANTAR_DEGREES) < 2.0)
+    ensure("ankles: the toes up to ${human.ANKLE_DORSAL_DEGREES} (${toes_up.x})", math.abs_float(toes_up.x + human.ANKLE_DORSAL_DEGREES) < 2.0)
+    ensure("ankles: about ${human.ANKLE_SIDE_DEGREES} to the side (${aside.y})", math.abs_float(aside.y - human.ANKLE_SIDE_DEGREES) < 2.0)
+    ensure("ankles: the sole rolled ${human.ANKLE_ROLL_DEGREES} (${rolled.z})", math.abs_float(math.abs_float(rolled.z) - human.ANKLE_ROLL_DEGREES) < 2.0)
+    // Ankles need a foot on each calf: refused without.
+    footless = human.human_reference_shape()
+    human.spec_at(footless, human.BONE_CALF_L).has_foot = false
+    world_id = world_with_ground()
+    h = Human { bones: null, filter_joints: null, filter_joint_count: 0, friction_torque: 0.0, is_spawned: false }
+    h.ankles = true
+    human.create_human_shaped(&h, world_id, math.transform_identity(), footless, 5.0, 1.0, 0.7, 1, null)
+    ensure("ankles: a calf without a foot is refused", h.is_spawned == false && physics_world.get_counters(world_id).body_count == 1)
+    dynamics.destroy_world(world_id)
+    free(footless)
+}
+
+// The issue's own measure (#152): a crouch tips the shin forward over the
+// foot, 25 degrees at a 50 degree knee bend. With the figure held still
+// and each calf turned forward about its ankle by `degrees`, the free feet
+// settle on the ground: within the ankle's range the soles lie as they did
+// standing (the reference's foot, fixed to the calf, would tilt by the
+// whole angle); past it the sole lifts by what the range cannot take.
+foot_tilt(h: *Human, foot_index: int) -> float {
+    foot = human.bone_at(h, foot_index)
+    c = dynamics.shape_get_capsule(foot.shape_id)
+    xf = Transform { p: dynamics.body_get_position(foot.body_id), q: dynamics.body_get_rotation(foot.body_id) }
+    axis = math.normalize(math.sub(math.transform_point(xf, c.center2), math.transform_point(xf, c.center1)))
+    return asin(math.abs_float(axis.y)) * math.RAD_TO_DEG
+}
+
+crouched_soles(degrees: float) -> Vec2 {
+    world_id = world_with_ground()
+    h = Human { bones: null, filter_joints: null, filter_joint_count: 0, friction_torque: 0.0, is_spawned: false }
+    h.ankles = true
+    human.create_human(&h, world_id, math.vec3(0.0, 0.01, 0.0), 0.0, 0.0, 0.0, 1, null)
+    rest_l = foot_tilt(&h, human.BONE_FOOT_L)
+    rest_r = foot_tilt(&h, human.BONE_FOOT_R)
+    i = 0
+    while i < human.BONE_COUNT {
+        dynamics.body_set_type(human.bone_at(&h, i).body_id, dynamics.BODY_STATIC)
+        i = i + 1
+    }
+    side = 0
+    while side < 2 {
+        calf_index = human.BONE_CALF_L
+        foot_index = human.BONE_FOOT_L
+        if side == 1 {
+            calf_index = human.BONE_CALF_R
+            foot_index = human.BONE_FOOT_R
+        }
+        calf = human.bone_at(&h, calf_index)
+        foot = human.bone_at(&h, foot_index)
+        frame = ankle_frame(&h, foot_index)
+        q = dynamics.body_get_rotation(calf.body_id)
+        p = dynamics.body_get_position(calf.body_id)
+        ankle = math.add(p, math.rotate_vector(q, foot.local_frame_a.p))
+        lateral = math.rotate_vector(frame, math.vec3_axis_y())
+        toe_way = math.rotate_vector(frame, math.vec3_axis_z())
+        forward = math.normalize(math.vec3(toe_way.x, 0.0, toe_way.z))
+        // The turn that takes the knee forward, over the toes.
+        turn = math.make_quat_from_axis_angle(lateral, degrees * math.DEG_TO_RAD)
+        if math.dot(math.sub(math.rotate_vector(turn, math.sub(p, ankle)), math.sub(p, ankle)), forward) < 0.0 {
+            turn = math.make_quat_from_axis_angle(lateral, 0.0 - degrees * math.DEG_TO_RAD)
+        }
+        dynamics.body_set_transform(calf.body_id, math.add(ankle, math.rotate_vector(turn, math.sub(p, ankle))), math.mul_quat(turn, q))
+        side = side + 1
+    }
+    s = 0
+    while s < 180 {
+        physics_world.step(world_id, 1.0 / 60.0, 4)
+        s = s + 1
+    }
+    tilt = Vec2 { x: foot_tilt(&h, human.BONE_FOOT_L) - rest_l, y: foot_tilt(&h, human.BONE_FOOT_R) - rest_r }
+    human.destroy_human(&h)
+    dynamics.destroy_world(world_id)
+    return tilt
+}
+
+test_ankle_crouch() {
+    shallow = crouched_soles(25.0)
+    ensure("crouch: shins 25 degrees over the feet, the soles lie as standing (${shallow.x}, ${shallow.y} degrees)",
+           math.abs_float(shallow.x) < 3.0 && math.abs_float(shallow.y) < 3.0)
+    deep = crouched_soles(45.0)
+    over = 45.0 - human.ANKLE_DORSAL_DEGREES
+    println("human: crouched soles against standing: shins at 25 degrees ${shallow.x} and ${shallow.y}, at 45 degrees ${deep.x} and ${deep.y}")
+    ensure("crouch: shins 45 degrees over, the soles lift by what the 35 cannot take (${deep.x}, ${deep.y} against ${over})",
+           math.abs_float(deep.x - over) < 4.0 && math.abs_float(deep.y - over) < 4.0)
+}
+
 main() {
     g_slab = hull.make_box_hull(20.0, 1.0, 20.0) as ptr
     test_build_and_fall()
@@ -826,6 +1039,8 @@ main() {
     test_hinges()
     test_hips()
     test_grip()
+    test_ankles()
+    test_ankle_crouch()
     hull.destroy_hull(g_slab as *HullData)
     println("human: ${checks} checks")
     if failures == 0 {
diff --git a/tests/test_restitution.ae b/tests/test_restitution.ae
new file mode 100644
index 0000000..b0a3cbc
--- /dev/null
+++ b/tests/test_restitution.ae
@@ -0,0 +1,913 @@
+// The reference's test_restitution.c (e77352c, #144): restitution measured
+// across scenes, each subtest printing its measurements before its gates,
+// every subtest run so one bad scene does not hide the rest. Head-on balls
+// at three density ratios, a ground bounce at sixteen impact phases, an
+// impactor on a column of dead balls with propagation off and on, a flat
+// landing and a spinning one, elastic drops from 40, 20 and 10 m, a cube
+// dropped flat four times, the threshold, a resting stack that must not
+// wake, the published normal velocity, the narrow-floor overshoot, the
+// impulse balance of a drop at every step, the energy of thirteen elastic
+// scenes, and the same scene at one, two and four workers.
+
+import std.string
+import aephysics.math
+import aephysics.basics
+import aephysics.hull
+import aephysics.manifold
+import aephysics.dynamics
+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
+// A flat landing converges its coupling in four stage iterations and the
+// default is two, so the tests that read a flat landing's coefficient or
+// spin pin the count.
+const RESTITUTION_ITERATIONS = 4
+// Well above the restitution threshold, so the bounce is always armed.
+const IMPACT_SPEED = 5.0
+
+var failures = 0
+var checks = 0
+
+ensure(name: string, ok: bool) {
+    checks = checks + 1
+    if !ok {
+        println("restitution: FAIL ${name}")
+        failures = failures + 1
+    }
+}
+
+// The hulls a scene's shapes hold, freed after its world.
+var g_hulls: ptr = null
+var g_hull_count = 0
+
+kept_hull(h: *HullData) -> *HullData {
+    hulls = g_hulls as ptr[]
+    hulls[g_hull_count] = h as ptr
+    g_hull_count = g_hull_count + 1
+    return h
+}
+
+end_world(world_id: WorldId) {
+    dynamics.destroy_world(world_id)
+    hulls = g_hulls as ptr[]
+    i = 0
+    while i < g_hull_count {
+        hull.destroy_hull(hulls[i] as *HullData)
+        i = i + 1
+    }
+    g_hull_count = 0
+}
+
+make_world_iterations(gravity_y: float, iterations: int) -> WorldId {
+    wdef = dynamics.default_world_def()
+    wdef.gravity = math.vec3(0.0, gravity_y, 0.0)
+    wdef.enable_sleep = false
+    wdef.restitution_iterations = iterations
+    return dynamics.create_world(&wdef)
+}
+
+make_world(gravity_y: float) -> WorldId {
+    wdef = dynamics.default_world_def()
+    wdef.gravity = math.vec3(0.0, gravity_y, 0.0)
+    wdef.enable_sleep = false
+    return dynamics.create_world(&wdef)
+}
+
+// Ground with its top at y = 0.
+make_ground(world_id: WorldId, restitution: float) {
+    bdef = dynamics.default_body_def()
+    bdef.position = math.vec3(0.0, 0.0 - 1.0, 0.0)
+    ground = dynamics.create_body(world_id, &bdef)
+    sdef = dynamics.default_shape_def()
+    sdef.base_material.friction = 0.0
+    sdef.base_material.restitution = restitution
+    dynamics.create_hull_shape(ground, &sdef, kept_hull(hull.make_box_hull(40.0, 1.0, 40.0)))
+}
+
+make_ball(world_id: WorldId, x: float, y: float, velocity_y: float, restitution: float) -> BodyId {
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(x, y, 0.0)
+    bdef.linear_velocity = math.vec3(0.0, velocity_y, 0.0)
+    body = dynamics.create_body(world_id, &bdef)
+    sdef = dynamics.default_shape_def()
+    sdef.base_material.friction = 0.0
+    sdef.base_material.restitution = restitution
+    dynamics.create_sphere_shape(body, &sdef, manifold.sphere(math.vec3_zero(), 0.5))
+    return body
+}
+
+steps(world_id: WorldId, count: int) {
+    i = 0
+    while i < count {
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        i = i + 1
+    }
+}
+
+// --- the coefficient --------------------------------------------------------------------------------------------
+
+var g_momentum_error = 0.0
+
+// Two free spheres closing head on with no gravity: the coefficient is the
+// ratio of relative normal speeds, the mass-free definition.
+measure_head_on(restitution: float, density_b: float) -> float {
+    world_id = make_world(0.0)
+    sdef = dynamics.default_shape_def()
+    sdef.density = 1.0
+    sdef.base_material.friction = 0.0
+    sdef.base_material.restitution = restitution
+    ball = manifold.sphere(math.vec3_zero(), 0.5)
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(0.0 - 1.0, 0.0, 0.0)
+    bdef.linear_velocity = math.vec3(0.5 * IMPACT_SPEED, 0.0, 0.0)
+    a = dynamics.create_body(world_id, &bdef)
+    dynamics.create_sphere_shape(a, &sdef, ball)
+    bdef.position = math.vec3(1.0, 0.0, 0.0)
+    bdef.linear_velocity = math.vec3(0.0 - 0.5 * IMPACT_SPEED, 0.0, 0.0)
+    b = dynamics.create_body(world_id, &bdef)
+    sdef.density = density_b
+    dynamics.create_sphere_shape(b, &sdef, ball)
+    mass_a = dynamics.body_get_mass(a)
+    mass_b = dynamics.body_get_mass(b)
+    momentum0 = mass_a * dynamics.body_get_linear_velocity(a).x + mass_b * dynamics.body_get_linear_velocity(b).x
+    steps(world_id, 120)
+    va = dynamics.body_get_linear_velocity(a).x
+    vb = dynamics.body_get_linear_velocity(b).x
+    momentum1 = mass_a * va + mass_b * vb
+    end_world(world_id)
+    g_momentum_error = (momentum1 - momentum0) / (mass_a + mass_b)
+    return (vb - va) / IMPACT_SPEED
+}
+
+test_head_on() {
+    restitutions = [ 0.0, 0.25, 0.5, 0.75, 1.0 ]
+    densities = [ 1.0, 10.0, 100.0 ]
+    i = 0
+    while i < 3 {
+        worst = 0.0
+        worst_momentum = 0.0
+        j = 0
+        while j < 5 {
+            measured = measure_head_on(restitutions[j], densities[i])
+            error = math.abs_float(measured - restitutions[j])
+            worst = math.max_float(worst, error)
+            if math.abs_float(g_momentum_error) > math.abs_float(worst_momentum) { worst_momentum = g_momentum_error }
+            ensure("head on: density ${densities[i]} e ${restitutions[j]} -> ${measured}", error <= 0.02)
+            j = j + 1
+        }
+        println("restitution: head on density ${densities[i]} worst error ${worst} momentum drift ${worst_momentum}")
+        ensure("head on: density ${densities[i]} momentum kept (${worst_momentum})", math.abs_float(worst_momentum) <= 0.0001)
+        i = i + 1
+    }
+}
+
+// A ball driven into the ground with no gravity; the gap moves where in
+// the step the impact lands.
+measure_ground_bounce(restitution: float, gap: float) -> float {
+    world_id = make_world(0.0)
+    make_ground(world_id, 0.0)
+    ball = make_ball(world_id, 0.0, 0.5 + gap, 0.0 - IMPACT_SPEED, restitution)
+    steps(world_id, 60)
+    speed = dynamics.body_get_linear_velocity(ball).y
+    end_world(world_id)
+    return speed / IMPACT_SPEED
+}
+
+test_phase() {
+    restitutions = [ 0.25, 0.5, 0.9 ]
+    j = 0
+    while j < 3 {
+        low = math.MAX_FLOAT
+        high = 0.0 - math.MAX_FLOAT
+        sum = 0.0
+        i = 0
+        while i < 16 {
+            gap = ((i as float) / 16.0) * IMPACT_SPEED * TIME_STEP
+            measured = measure_ground_bounce(restitutions[j], gap)
+            low = math.min_float(low, measured)
+            high = math.max_float(high, measured)
+            sum = sum + measured
+            i = i + 1
+        }
+        mean = sum / 16.0
+        println("restitution: phase e ${restitutions[j]} -> mean ${mean} spread ${high - low} [${low}, ${high}]")
+        ensure("phase: e ${restitutions[j]} spread ${high - low} and mean ${mean} within 0.05", high - low <= 0.05 && math.abs_float(mean - restitutions[j]) <= 0.05)
+        j = j + 1
+    }
+}
+
+// An impactor driven onto a column of dead balls on the ground: a
+// supported target has infinite effective mass along the normal, so the
+// coefficient must not depend on the column. Without propagation the
+// stage sweeps only the contacts with restitution and the column gives way
+// a step late; with it the ground answers within the stage. Never above
+// the free coefficient, never collapsed.
+measure_supported(restitution: float, supports: int, propagate: bool) -> float {
+    wdef = dynamics.default_world_def()
+    wdef.gravity = math.vec3_zero()
+    wdef.enable_sleep = false
+    wdef.enable_restitution_propagation = propagate
+    wdef.restitution_iterations = RESTITUTION_ITERATIONS
+    world_id = dynamics.create_world(&wdef)
+    make_ground(world_id, 0.0)
+    i = 0
+    while i < supports {
+        make_ball(world_id, 0.0, 0.5 + 1.0 * (i as float), 0.0, 0.0)
+        i = i + 1
+    }
+    y = 0.5 + 1.0 * (supports as float) + 0.5 * IMPACT_SPEED * TIME_STEP
+    ball = make_ball(world_id, 0.0, y, 0.0 - IMPACT_SPEED, restitution)
+    steps(world_id, 60)
+    speed = dynamics.body_get_linear_velocity(ball).y
+    end_world(world_id)
+    return speed / IMPACT_SPEED
+}
+
+test_supported() {
+    restitutions = [ 0.5, 0.9 ]
+    propagate = 0
+    while propagate <= 1 {
+        floor_ratio = 0.5
+        if propagate == 1 { floor_ratio = 0.75 }
+        j = 0
+        while j < 2 {
+            n = 0
+            while n <= 3 {
+                measured = measure_supported(restitutions[j], n, propagate == 1)
+                println("restitution: supported e ${restitutions[j]} supports ${n} propagation ${propagate} -> ${measured}")
+                ensure("supported: e ${restitutions[j]} on ${n} propagation ${propagate}, ${measured} not above e nor below ${floor_ratio} e",
+                       measured <= restitutions[j] + 0.01 && measured >= floor_ratio * restitutions[j])
+                n = n + 1
+            }
+            j = j + 1
+        }
+        propagate = propagate + 1
+    }
+}
+
+// A flat box landing on four corners at once with no gravity.
+measure_flat(restitution: float, sub_steps: int, spin_z: float) -> Vec3 {
+    world_id = make_world_iterations(0.0, RESTITUTION_ITERATIONS)
+    make_ground(world_id, 0.0)
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(0.0, 0.25 + 0.5 * IMPACT_SPEED * TIME_STEP, 0.0)
+    bdef.linear_velocity = math.vec3(0.0, 0.0 - IMPACT_SPEED, 0.0)
+    bdef.angular_velocity = math.vec3(0.0, 0.0, spin_z)
+    box = dynamics.create_body(world_id, &bdef)
+    sdef = dynamics.default_shape_def()
+    sdef.density = 1.0
+    sdef.base_material.friction = 0.0
+    sdef.base_material.restitution = restitution
+    dynamics.create_hull_shape(box, &sdef, kept_hull(hull.make_box_hull(1.0, 0.25, 1.0)))
+    i = 0
+    while i < 60 {
+        physics_world.step(world_id, TIME_STEP, sub_steps)
+        i = i + 1
+    }
+    v = dynamics.body_get_linear_velocity(box)
+    w = dynamics.body_get_angular_velocity(box)
+    end_world(world_id)
+    // The speed ratio, the spin's size, and the spin about z.
+    return math.vec3(v.y / IMPACT_SPEED, math.length(w), w.z)
+}
+
+// One sweep over four coplanar points leaves the box spinning; four stage
+// iterations converge the coupling. Only the shipping sub-step count gates.
+test_flat_landing() {
+    restitutions = [ 0.5, 0.9 ]
+    counts = [ 4, 8 ]
+    j = 0
+    while j < 2 {
+        k = 0
+        while k < 2 {
+            m = measure_flat(restitutions[j], counts[k], 0.0)
+            println("restitution: flat e ${restitutions[j]} substeps ${counts[k]} -> ${m.x} spin ${m.y}")
+            if counts[k] == SUB_STEP_COUNT {
+                ensure("flat landing: e ${restitutions[j]} -> ${m.x} within 0.01, spin ${m.y} under 0.03", math.abs_float(m.x - restitutions[j]) <= 0.01 && m.y <= 0.03)
+            }
+            k = k + 1
+        }
+        j = j + 1
+    }
+}
+
+// A flat box carrying spin, perfectly elastic: the linear and angular
+// velocity come back reversed.
+test_spin() {
+    spins = [ 0.0, 0.5, 1.0, 2.0 ]
+    j = 0
+    while j < 4 {
+        m = measure_flat(1.0, SUB_STEP_COUNT, spins[j])
+        speed = m.x * IMPACT_SPEED
+        println("restitution: spin in ${spins[j]} -> vy ${speed} (want ${IMPACT_SPEED}) wz ${m.z} (want ${0.0 - spins[j]})")
+        ensure("spin: ${spins[j]} reversed (vy ${speed}, wz ${m.z})", math.abs_float(speed - IMPACT_SPEED) <= 0.05 && math.abs_float(m.z + spins[j]) <= 0.05)
+        j = j + 1
+    }
+}
+
+// --- drops --------------------------------------------------------------------------------------------------------
+
+// A perfectly elastic ball dropped under gravity; its apexes. Only heights
+// where continuous collision engages: it lands the ball on the surface, so
+// the apex holds.
+test_drop() {
+    heights = [ 40.0, 20.0, 10.0 ]
+    apexes_block = calloc(6, 8)
+    apexes = apexes_block as float[]
+    j = 0
+    while j < 3 {
+        world_id = make_world(0.0 - 10.0)
+        make_ground(world_id, 0.0)
+        ball = make_ball(world_id, 0.0, heights[j], 0.0, 1.0)
+        count = 0
+        previous = 0.0
+        i = 0
+        while i < 4000 && count < 6 {
+            physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+            speed = dynamics.body_get_linear_velocity(ball).y
+            if previous > 0.0 && speed <= 0.0 {
+                apexes[count] = dynamics.body_get_position(ball).y
+                count = count + 1
+            }
+            previous = speed
+            i = i + 1
+        }
+        end_world(world_id)
+        line = "restitution: drop ${heights[j]} ->"
+        highest = 0.0
+        k = 0
+        while k < count {
+            line = string.concat(line, " ${apexes[k]}")
+            highest = math.max_float(highest, apexes[k])
+            k = k + 1
+        }
+        println(line)
+        ensure("drop ${heights[j]}: six apexes (${count})", count == 6)
+        if count == 6 {
+            ensure("drop ${heights[j]}: no apex above the first, the first not above the drop, the last above 0.8 of the first",
+                   highest <= 1.001 * apexes[0] && apexes[0] <= 1.001 * heights[j] && apexes[count - 1] >= 0.8 * apexes[0])
+        }
+        j = j + 1
+    }
+    free(apexes_block)
+}
+
+// A cube dropped flat from 10 m with a tight continuous safety factor,
+// elastic and frictionless: it comes back to the drop with no spin, four
+// times.
+test_single_box() {
+    world_id = make_world_iterations(0.0 - 10.0, RESTITUTION_ITERATIONS)
+    make_ground(world_id, 0.0)
+    sdef = dynamics.default_shape_def()
+    sdef.density = 1.0
+    sdef.base_material.friction = 0.0
+    sdef.base_material.restitution = 1.0
+    drop = 10.0
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(0.0, drop, 0.0)
+    bdef.safety_factor = 0.01
+    box = dynamics.create_body(world_id, &bdef)
+    dynamics.create_hull_shape(box, &sdef, kept_hull(hull.make_cube_hull(0.5)))
+    first_spin = 0.0
+    apexes_block = calloc(8, 8)
+    apexes = apexes_block as float[]
+    count = 0
+    previous = 0.0
+    i = 0
+    while i < 1200 && count < 4 {
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        speed = dynamics.body_get_linear_velocity(box).y
+        if count == 0 && previous <= 0.0 && speed > 0.0 { first_spin = math.length(dynamics.body_get_angular_velocity(box)) }
+        if previous > 0.0 && speed <= 0.0 {
+            apexes[count] = dynamics.body_get_position(box).y
+            apexes[4 + count] = math.length(dynamics.body_get_angular_velocity(box))
+            count = count + 1
+        }
+        previous = speed
+        i = i + 1
+    }
+    end_world(world_id)
+    line = "restitution: single box first bounce spin ${first_spin} apexes"
+    k = 0
+    while k < count {
+        line = string.concat(line, " ${apexes[k]} (spin ${apexes[4 + k]})")
+        k = k + 1
+    }
+    println(line)
+    ensure("single box: four apexes (${count})", count == 4)
+    ensure("single box: no spin off the first bounce (${first_spin})", first_spin <= 0.01)
+    k = 0
+    while k < count {
+        ensure("single box: apex ${k} at the drop (${apexes[k]})", apexes[k] >= 0.95 * drop && apexes[k] <= 1.001 * drop)
+        k = k + 1
+    }
+    free(apexes_block)
+}
+
+// Below the threshold no bounce; above it, nearly all of it.
+test_threshold() {
+    wdef = dynamics.default_world_def()
+    threshold = wdef.restitution_threshold
+    scales = [ 0.5, 0.9, 1.5, 4.0 ]
+    j = 0
+    while j < 4 {
+        speed = scales[j] * threshold
+        world_id = make_world(0.0)
+        make_ground(world_id, 0.0)
+        ball = make_ball(world_id, 0.0, 0.5 + 0.5 * speed * TIME_STEP, 0.0 - speed, 1.0)
+        steps(world_id, 60)
+        ratio = dynamics.body_get_linear_velocity(ball).y / speed
+        end_world(world_id)
+        println("restitution: threshold ${scales[j]}x -> ${ratio}")
+        if scales[j] < 1.0 { ensure("threshold: ${scales[j]}x does not bounce (${ratio})", ratio <= 0.05) }
+        else { ensure("threshold: ${scales[j]}x bounces (${ratio})", ratio >= 0.9) }
+        j = j + 1
+    }
+}
+
+// High restitution does not wake a settled stack back up.
+test_resting() {
+    world_id = make_world(0.0 - 10.0)
+    make_ground(world_id, 0.9)
+    ids_block = calloc(10, sizeof(BodyId))
+    ids = ids_block as BodyId[]
+    sdef = dynamics.default_shape_def()
+    sdef.density = 1.0
+    sdef.base_material.friction = 0.6
+    sdef.base_material.restitution = 0.9
+    cube = kept_hull(hull.make_cube_hull(0.5))
+    i = 0
+    while i < 10 {
+        bdef = dynamics.default_body_def()
+        bdef.kind = dynamics.BODY_DYNAMIC
+        bdef.position = math.vec3(0.0, 0.5 + 1.0 * (i as float), 0.0)
+        ids[i] = dynamics.create_body(world_id, &bdef)
+        dynamics.create_hull_shape(ids[i], &sdef, cube)
+        i = i + 1
+    }
+    steps(world_id, 240)
+    settled = dynamics.body_get_position(ids[9]).y
+    drift = 0.0
+    peak = 0.0
+    i = 0
+    while i < 300 {
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        drift = math.max_float(drift, math.abs_float(dynamics.body_get_position(ids[9]).y - settled))
+        j = 0
+        while j < 10 {
+            peak = math.max_float(peak, math.length(dynamics.body_get_linear_velocity(ids[j])))
+            j = j + 1
+        }
+        i = i + 1
+    }
+    end_world(world_id)
+    free(ids_block)
+    println("restitution: resting drift ${drift} peak speed ${peak}")
+    ensure("resting: the stack does not drift (${drift})", drift < 0.01)
+    ensure("resting: nor wake (${peak})", peak < 0.05)
+}
+
+// --- what a manifold publishes --------------------------------------------------------------------------------
+
+// The approach speed is published only for contacts with hit events: zero
+// otherwise, so a reader can tell "not measured" from "measured zero".
+first_touch_normal_velocity(hit_events: bool) -> float {
+    world_id = make_world(0.0)
+    make_ground(world_id, 0.0)
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(0.0, 0.6, 0.0)
+    ball = dynamics.create_body(world_id, &bdef)
+    sdef = dynamics.default_shape_def()
+    sdef.base_material.friction = 0.0
+    sdef.base_material.restitution = 0.0
+    sdef.enable_hit_events = hit_events
+    dynamics.create_sphere_shape(ball, &sdef, manifold.sphere(math.vec3_zero(), 0.5))
+    data_block = calloc(1, sizeof(ContactData))
+    data = data_block as ContactData[]
+    value = math.MAX_FLOAT
+    i = 0
+    while i < 20 {
+        dynamics.body_set_linear_velocity(ball, math.vec3(0.0, 0.0 - IMPACT_SPEED, 0.0))
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        if dynamics.body_get_contact_data(ball, data, 1) == 1 && data[0].manifold_count > 0 {
+            m = data[0].manifolds as *Manifold
+            if m.point_count > 0 {
+                value = m.point0.normal_velocity as float
+                break
+            }
+        }
+        i = i + 1
+    }
+    free(data_block)
+    end_world(world_id)
+    return value
+}
+
+test_normal_velocity() {
+    quiet = first_touch_normal_velocity(false)
+    published = first_touch_normal_velocity(true)
+    println("restitution: normal velocity hit events off ${quiet} on ${published} (closing at ${0.0 - IMPACT_SPEED})")
+    ensure("normal velocity: zero without hit events (${quiet})", quiet == 0.0)
+    ensure("normal velocity: the approach with them (${published})", published <= 0.0 - 0.8 * IMPACT_SPEED)
+}
+
+// --- energy ---------------------------------------------------------------------------------------------------
+
+// Mechanical energy as symplectic Euler keeps it: kinetic plus potential
+// plus half a sub-step of gravity's power, so free flight is flat to
+// round-off and a gain or a loss at a bounce is real.
+measure_energy(body: BodyId, world_id: WorldId) -> float {
+    md = dynamics.body_get_mass_data(body)
+    v = dynamics.body_get_linear_velocity(body)
+    gravity = math.mul_sv(dynamics.body_get_gravity_scale(body), physics_world.get_gravity(world_id))
+    h = TIME_STEP / (SUB_STEP_COUNT as float)
+    w_local = math.inv_rotate_vector(dynamics.body_get_rotation(body), dynamics.body_get_angular_velocity(body))
+    kinetic = 0.5 * md.mass * math.dot(v, v) + 0.5 * math.dot(w_local, math.mul_mv(md.inertia, w_local))
+    potential = 0.0 - md.mass * math.dot(gravity, dynamics.body_get_world_center(body))
+    phase = 0.5 * h * md.mass * math.dot(v, gravity)
+    return kinetic + potential + phase
+}
+
+struct Overshoot {
+    first_apex: float
+    first_energy: float
+    peak_energy: float
+    speed_ratio: float
+    bounces: int
+}
+
+// The reference's Restitution Overshoot sample: a unit cube dropped flat
+// on a floor narrower than itself, perfectly elastic, the contact points
+// inboard of the corners. It can come back no higher than it fell and
+// gain no energy.
+measure_overshoot(continuous: bool) -> Overshoot {
+    world_id = make_world_iterations(0.0 - 10.0, RESTITUTION_ITERATIONS)
+    fdef = dynamics.default_body_def()
+    fdef.position = math.vec3(0.0, 0.0 - 0.25, 0.0)
+    floor_body = dynamics.create_body(world_id, &fdef)
+    fsdef = dynamics.default_shape_def()
+    dynamics.create_hull_shape(floor_body, &fsdef, kept_hull(hull.make_box_hull(0.375, 0.25, 0.375)))
+    drop = 10.0
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(0.0, drop, 0.0)
+    if continuous { bdef.safety_factor = 0.1 }
+    box = dynamics.create_body(world_id, &bdef)
+    sdef = dynamics.default_shape_def()
+    sdef.base_material.restitution = 1.0
+    dynamics.create_hull_shape(box, &sdef, kept_hull(hull.make_cube_hull(0.5)))
+    r = Overshoot { first_apex: 0.0, first_energy: 0.0, peak_energy: 0.0, speed_ratio: 0.0, bounces: 0 }
+    start = measure_energy(box, world_id)
+    impact = 0.0
+    rebound = 0.0
+    bounce_height = 0.0
+    previous = 0.0
+    i = 0
+    while i < 600 && r.bounces < 2 {
+        before = dynamics.body_get_linear_velocity(box).y
+        height_before = dynamics.body_get_position(box).y
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        speed = dynamics.body_get_linear_velocity(box).y
+        if rebound == 0.0 && before < 0.0 && speed > 0.0 {
+            impact = 0.0 - before
+            rebound = speed
+            bounce_height = height_before
+        }
+        if rebound != 0.0 {
+            ratio = measure_energy(box, world_id) / start
+            r.peak_energy = math.max_float(r.peak_energy, ratio)
+            if previous > 0.0 && speed <= 0.0 {
+                if r.bounces == 0 {
+                    r.first_apex = dynamics.body_get_position(box).y
+                    r.first_energy = ratio
+                }
+                r.bounces = r.bounces + 1
+            }
+            previous = speed
+        }
+        i = i + 1
+    }
+    end_world(world_id)
+    if impact > 0.0 { r.speed_ratio = rebound / impact }
+    println("restitution: overshoot continuous ${continuous} apex ${r.first_apex} of ${drop} at y ${bounce_height}, speed ratio ${r.speed_ratio}, energy ${r.first_energy} then ${r.peak_energy}")
+    return r
+}
+
+// Never more speed back than received, in both runs. With continuous
+// collision the cube lands on the surface and comes in under its start;
+// without, it lands 0.1 m deep and climbs back through that depth for
+// free, so only its first bounce is gated.
+test_overshoot() {
+    discrete = measure_overshoot(false)
+    continuous = measure_overshoot(true)
+    ensure("overshoot: no more speed back than received (${discrete.speed_ratio}, ${continuous.speed_ratio})",
+           discrete.speed_ratio <= 1.00001 && continuous.speed_ratio <= 1.00001)
+    ensure("overshoot: continuous, no higher than the drop and no energy gained (${continuous.first_apex}, ${continuous.peak_energy})",
+           continuous.first_apex <= 10.0 && continuous.peak_energy <= 1.001)
+    ensure("overshoot: discrete, its first bounce within its bounds (${discrete.first_apex}, ${discrete.first_energy})",
+           discrete.first_apex <= 10.15 && discrete.first_energy <= 1.02)
+    ensure("overshoot: the bounces counted (${discrete.bounces}, ${continuous.bounces})", discrete.bounces >= 2 && continuous.bounces >= 1)
+}
+
+struct ImpulseRun {
+    worst_error: float
+    approach_speed: float
+    first_separation: float
+    bounce_impulse: float
+    contact_steps: int
+    toi_steps: int
+    toi_impulse_steps: int
+}
+
+// A ball dropped under gravity with the contact's impulse read every step:
+// gravity taken out, the change of momentum over a step is what the
+// contact applied, and the total normal impulse must say so. A time of
+// impact step is left out: no contact was solved on it.
+measure_drop_impulse(restitution: float, drop: float) -> ImpulseRun {
+    world_id = make_world(0.0 - 10.0)
+    make_ground(world_id, 0.0)
+    ball = make_ball(world_id, 0.0, 0.5 + drop, 0.0, restitution)
+    w = dynamics.get_world(world_id.index1 - 1)
+    mass = dynamics.body_get_mass(ball)
+    gravity_y = physics_world.get_gravity(world_id).y
+    r = ImpulseRun { worst_error: 0.0, approach_speed: 0.0, first_separation: 0.0, bounce_impulse: 0.0, contact_steps: 0, toi_steps: 0, toi_impulse_steps: 0 }
+    touched = false
+    bouncing = false
+    data_block = calloc(4, sizeof(ContactData))
+    data = data_block as ContactData[]
+    i = 0
+    while i < 240 {
+        before = dynamics.body_get_linear_velocity(ball).y
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        after = dynamics.body_get_linear_velocity(ball).y
+        measured = 0.0
+        separation = 0.0
+        contacts = dynamics.body_get_contact_data(ball, data, 4)
+        c = 0
+        while c < contacts {
+            m = 0
+            while m < data[c].manifold_count {
+                mf = (data[c].manifolds + m * sizeof(Manifold)) as *Manifold
+                p = 0
+                while p < mf.point_count {
+                    mp = dynamics.manifold_point_at(mf, p)
+                    measured = measured + (mp.total_normal_impulse as float)
+                    separation = mp.separation as float
+                    p = p + 1
+                }
+                m = m + 1
+            }
+            c = c + 1
+        }
+        flags = dynamics.get_body_sim(w, dynamics.get_body_full_id(w, ball)).flags
+        if (flags & dynamics.HAD_TIME_OF_IMPACT) != 0 {
+            r.toi_steps = r.toi_steps + 1
+            if measured != 0.0 { r.toi_impulse_steps = r.toi_impulse_steps + 1 }
+        } else {
+            expected = mass * (after - before) - mass * gravity_y * TIME_STEP
+            r.worst_error = math.max_float(r.worst_error, math.abs_float(measured - expected))
+            if contacts > 0 {
+                if touched == false {
+                    touched = true
+                    bouncing = true
+                    r.approach_speed = 0.0 - before
+                    r.first_separation = separation
+                }
+                if bouncing {
+                    r.bounce_impulse = r.bounce_impulse + measured
+                    r.contact_steps = r.contact_steps + 1
+                }
+            } else {
+                bouncing = false
+            }
+        }
+        i = i + 1
+    }
+    free(data_block)
+    end_world(world_id)
+    return r
+}
+
+// Heights either side of the continuous threshold. The step balance is the
+// gate; the bounce's total is bracketed by the reversal at the coefficient
+// and the weight carried while the contact lasted.
+test_impulse() {
+    restitutions = [ 0.0, 0.25, 0.5, 0.75, 1.0 ]
+    heights = [ 1.0, 5.0, 20.0, 45.0 ]
+    world_id = make_world(0.0 - 10.0)
+    probe = make_ball(world_id, 0.0, 0.5, 0.0, 0.0)
+    mass = dynamics.body_get_mass(probe)
+    fast_speed = dynamics.body_get_safety_factor(probe) * dynamics.body_get_min_extent(probe) / TIME_STEP
+    end_world(world_id)
+    weight_impulse = mass * 10.0 * TIME_STEP
+    j = 0
+    while j < 4 {
+        impact_speed = sqrt(20.0 * heights[j])
+        expect_toi = impact_speed > fast_speed
+        tolerance = 0.00001 * mass * (10.0 + impact_speed)
+        k = 0
+        while k < 5 {
+            r = measure_drop_impulse(restitutions[k], heights[j])
+            reversal = (1.0 + restitutions[k]) * mass * r.approach_speed
+            lower = reversal - weight_impulse
+            upper = reversal + weight_impulse * (r.contact_steps as float)
+            slack = (r.contact_steps as float) * tolerance
+            println("restitution: impulse drop ${heights[j]} e ${restitutions[k]} toi ${r.toi_steps} at ${r.first_separation} -> worst step error ${r.worst_error}, bounce ${r.bounce_impulse} in [${lower}, ${upper}] over ${r.contact_steps} steps")
+            ensure("impulse: drop ${heights[j]} e ${restitutions[k]}, every step's impulse balanced (${r.worst_error} against ${tolerance})", r.worst_error <= tolerance)
+            ensure("impulse: drop ${heights[j]} e ${restitutions[k]}, contact steps and none on a time of impact", r.contact_steps > 0 && r.toi_impulse_steps == 0)
+            ensure("impulse: drop ${heights[j]} e ${restitutions[k]}, the bounce bracketed", r.bounce_impulse >= lower - slack && r.bounce_impulse <= upper + slack)
+            ensure("impulse: drop ${heights[j]} e ${restitutions[k]}, continuous collision as the speed says (${r.toi_steps > 0} at ${impact_speed} against ${fast_speed})", (r.toi_steps > 0) == expect_toi)
+            k = k + 1
+        }
+        j = j + 1
+    }
+}
+
+struct EnergyScene {
+    safety_factor: float
+    friction: float
+    angle: float
+    spin: float
+    height: float
+    velocity_y: float
+    half_width: float
+    half_height: float
+    bullet: bool
+    stack: bool
+    step_count: int
+}
+
+var g_peak_step = 0
+
+// The peak energy over the run against the start, every dynamic body: an
+// elastic bounce may return at most what the contact took.
+peak_energy(scene: EnergyScene) -> float {
+    world_id = make_world(0.0 - 10.0)
+    gdef = dynamics.default_body_def()
+    gdef.position = math.vec3(0.0, 0.0 - 1.0, 0.0)
+    ground = dynamics.create_body(world_id, &gdef)
+    gsdef = dynamics.default_shape_def()
+    gsdef.base_material.friction = scene.friction
+    dynamics.create_hull_shape(ground, &gsdef, kept_hull(hull.make_box_hull(40.0, 1.0, 40.0)))
+    sdef = dynamics.default_shape_def()
+    sdef.density = 1.0
+    sdef.base_material.friction = scene.friction
+    sdef.base_material.restitution = 1.0
+    box = kept_hull(hull.make_box_hull(scene.half_width, scene.half_height, scene.half_width))
+    ids_block = calloc(2, sizeof(BodyId))
+    ids = ids_block as BodyId[]
+    count = 0
+    if scene.stack {
+        bdef = dynamics.default_body_def()
+        bdef.kind = dynamics.BODY_DYNAMIC
+        bdef.position = math.vec3(0.0, scene.half_height, 0.0)
+        bdef.safety_factor = scene.safety_factor
+        ids[count] = dynamics.create_body(world_id, &bdef)
+        dynamics.create_hull_shape(ids[count], &sdef, box)
+        count = count + 1
+    }
+    bdef = dynamics.default_body_def()
+    bdef.kind = dynamics.BODY_DYNAMIC
+    bdef.position = math.vec3(0.0, scene.height, 0.0)
+    bdef.rotation = math.make_quat_from_axis_angle(math.vec3_axis_z(), scene.angle)
+    bdef.angular_velocity = math.vec3(0.0, 0.0, scene.spin)
+    bdef.linear_velocity = math.vec3(0.0, scene.velocity_y, 0.0)
+    bdef.safety_factor = scene.safety_factor
+    bdef.is_bullet = scene.bullet
+    ids[count] = dynamics.create_body(world_id, &bdef)
+    dynamics.create_hull_shape(ids[count], &sdef, box)
+    count = count + 1
+    start = 0.0
+    k = 0
+    while k < count {
+        start = start + measure_energy(ids[k], world_id)
+        k = k + 1
+    }
+    peak = start
+    g_peak_step = 0 - 1
+    i = 0
+    while i < scene.step_count {
+        physics_world.step(world_id, TIME_STEP, SUB_STEP_COUNT)
+        energy = 0.0
+        k = 0
+        while k < count {
+            energy = energy + measure_energy(ids[k], world_id)
+            k = k + 1
+        }
+        if energy > peak {
+            peak = energy
+            g_peak_step = i
+        }
+        i = i + 1
+    }
+    end_world(world_id)
+    free(ids_block)
+    return peak / start
+}
+
+energy_scene(safety_factor: float, friction: float, angle: float, spin: float, height: float, velocity_y: float, half_width: float, half_height: float, bullet: bool, stack: bool, step_count: int) -> EnergyScene {
+    return EnergyScene { safety_factor: safety_factor, friction: friction, angle: angle, spin: spin, height: height, velocity_y: velocity_y,
+                         half_width: half_width, half_height: half_height, bullet: bullet, stack: stack, step_count: step_count }
+}
+
+// Perfectly elastic, with friction, dropped every way that once gained
+// energy: the energy never rises above the start (0.5% with continuous
+// collision, 2% discrete).
+test_energy() {
+    names = [ "spinning edge", "spinning edge fast", "tumble", "tumble discrete", "flat with spin", "flat frictionless", "edge",
+              "two box", "two box discrete", "tall flat frictionless", "tall flat friction", "narrow flat frictionless", "narrow flat friction" ]
+    scenes_block = calloc(13, sizeof(EnergyScene))
+    scenes = scenes_block as EnergyScene[]
+    scenes[0] = energy_scene(0.01, 0.6, 0.0 - 1.43, 7.9, 0.7, 0.0 - 1.9, 0.5, 0.5, false, false, 300)
+    scenes[1] = energy_scene(0.01, 0.6, 0.0 - 0.8, 12.0, 3.0, 0.0 - 8.0, 0.5, 0.5, false, false, 300)
+    scenes[2] = energy_scene(0.01, 0.6, 0.0 - 0.3, 3.0, 8.0, 0.0, 0.5, 0.5, false, false, 600)
+    scenes[3] = energy_scene(0.5, 0.6, 0.0 - 0.3, 3.0, 8.0, 0.0, 0.5, 0.5, false, false, 600)
+    scenes[4] = energy_scene(0.01, 0.6, 0.0, 2.0, 2.0, 0.0 - 5.0, 0.5, 0.5, false, false, 300)
+    scenes[5] = energy_scene(0.01, 0.0, 0.0, 0.0, 10.0, 0.0, 0.5, 0.5, false, false, 1200)
+    scenes[6] = energy_scene(0.01, 0.6, 0.0 - 0.2, 0.0, 5.0, 0.0, 0.5, 0.5, false, false, 600)
+    scenes[7] = energy_scene(0.02, 0.6, 0.0, 0.0, 4.0, 0.0, 0.5, 0.5, true, true, 1500)
+    scenes[8] = energy_scene(0.5, 0.6, 0.0, 0.0, 4.0, 0.0, 0.5, 0.5, false, true, 1500)
+    scenes[9] = energy_scene(0.01, 0.0, 0.0, 0.0, 10.0, 0.0, 0.5, 1.0, false, false, 1200)
+    scenes[10] = energy_scene(0.01, 0.6, 0.0, 0.0, 10.0, 0.0, 0.5, 1.0, false, false, 1200)
+    scenes[11] = energy_scene(0.01, 0.0, 0.0, 0.0, 10.0, 0.0, 0.25, 1.0, false, false, 1200)
+    scenes[12] = energy_scene(0.01, 0.6, 0.0, 0.0, 10.0, 0.0, 0.25, 1.0, false, false, 1200)
+    i = 0
+    while i < 13 {
+        ratio = peak_energy(scenes[i])
+        limit = 0.02
+        if scenes[i].safety_factor < 0.1 { limit = 0.005 }
+        println("restitution: energy ${names[i]} peak ${ratio} at step ${g_peak_step}")
+        ensure("energy: ${names[i]} never above the start (${ratio}, within ${limit})", ratio <= 1.0 + limit)
+        i = i + 1
+    }
+    free(scenes_block)
+}
+
+// --- workers ----------------------------------------------------------------------------------------------------
+
+// The per-point state the restitution stage keeps survives the split
+// across workers: one, two and four, hashed alike.
+worker_scene(workers: int) -> long {
+    wdef = dynamics.default_world_def()
+    wdef.gravity = math.vec3(0.0, 0.0 - 10.0, 0.0)
+    wdef.enable_sleep = false
+    wdef.worker_count = workers
+    world_id = dynamics.create_world(&wdef)
+    make_ground(world_id, 0.6)
+    i = 0
+    while i < 20 {
+        make_ball(world_id, 0.0 - 10.0 + 1.05 * (i as float), 3.0 + 0.13 * (i as float), 0.0, 0.6)
+        i = i + 1
+    }
+    steps(world_id, 200)
+    hash = dynamics.hash_world_state(dynamics.get_world(world_id.index1 - 1))
+    end_world(world_id)
+    return hash
+}
+
+test_worker_parity() {
+    one = worker_scene(1)
+    two = worker_scene(2)
+    four = worker_scene(4)
+    println("restitution: worker hashes ${one} ${two} ${four}")
+    ensure("workers: two alike to one", one == two)
+    ensure("workers: four alike to one", one == four)
+}
+
+main() {
+    g_hulls = calloc(16, 8)
+    test_head_on()
+    test_phase()
+    test_supported()
+    test_flat_landing()
+    test_spin()
+    test_drop()
+    test_single_box()
+    test_threshold()
+    test_resting()
+    test_normal_velocity()
+    test_overshoot()
+    test_energy()
+    test_impulse()
+    test_worker_parity()
+    free(g_hulls)
+
+    println("restitution: ${checks} checks")
+    if failures == 0 {
+        println("restitution: all checks passed")
+    } else {
+        println("restitution: ${failures} failure(s)")
+        exit(1)
+    }
+}