The reference's benchmark suite and missing tests: nine scenes, SAT runs, test_collision, test_hash, the hash reads every byte (#141, #145)

aephysics/basics/module.ae

@@ -370,10 +370,14 @@ map_bytes(m: *LongMap) -> int { return m.capacity * sizeof(MapItem) }
 
 // --- a hash over a block --------------------------------------------------------
 
-// This engine's hash of a block of 8-byte words (the count rounded down);
-// never zero. The reference uses rapidhash for the same purpose: identity
-// of hull and mesh data.
+// This engine's hash of a block of bytes; never zero. The reference uses
+// rapidhash for the same purpose: identity of hull and mesh data. Every
+// byte and the length reach it (test_hash.c's properties, #145): the whole
+// words through four lanes, the tail of fewer than eight bytes as one more
+// word, the byte count folded in at the end.
 hash_bytes(block: ptr, byte_count: int) -> long {
+    // An empty block has nothing to mix: the reserved value, the reference's.
+    if byte_count <= 0 { return 1 as long }
     // Four lanes take every fourth word, so their mixes (each a chain of
     // multiplies) overlap instead of waiting on one another: a word at a
     // time through one chain was most of building a box hull.
@@ -400,8 +404,19 @@ hash_bytes(block: ptr, byte_count: int) -> long {
         a = mix32(mix32(a ^ (w & m32)) ^ lsr(w, 32))
         i = i + 1
     }
-    // The lanes and the length folded into two halves of the result.
-    low = mix32(mix32(a ^ (count as long)) ^ b)
+    // The bytes after the last whole word, little-endian into one word.
+    tail_count = byte_count - 8 * count
+    if tail_count > 0 {
+        w = 0 as long
+        k = tail_count - 1
+        while k >= 0 {
+            w = (w << 8) | (native.load_u8(block, (8 * count + k) as long) as long)
+            k = k - 1
+        }
+        b = mix32(mix32(b ^ (w & m32)) ^ lsr(w, 32))
+    }
+    // The lanes and the length in bytes folded into two halves of the result.
+    low = mix32(mix32(a ^ (byte_count as long)) ^ b)
     high = mix32(mix32(c ^ low) ^ d)
     // high kept to 31 bits: a shift into the sign bit is undefined in C, and
     // GCC compiled it differently where the hash was inlined.

aephysics/hull/module.ae

@@ -30,7 +30,7 @@ exports (
     hull_vertices, hull_points, hull_edges, hull_planes, hull_faces,
     create_hull, destroy_hull, clone_hull, clone_and_transform_hull, is_valid_hull,
     hash_hull_data, compare_hull_data,
-    create_cylinder, create_cone, create_rock,
+    create_cylinder, create_cone, create_rock, create_complex_hull,
     make_box_hull, make_cube_hull, make_offset_box_hull, make_transformed_box_hull,
     make_scaled_box_hull, scale_box,
     find_hull_support_vertex, find_hull_support_face,
@@ -1707,6 +1707,45 @@ create_rock(radius: float) -> *HullData {
     return h
 }
 
+// PEEL's BasicRandom: a 32-bit linear congruential generator.
+struct PeelRandom {
+    state: long
+}
+
+// The next value in [-0.5, 0.5], in single precision as PEEL's.
+peel_random_float(g: *PeelRandom) -> f32 {
+    g.state = (g.state * (2147001325 as long) + (715136305 as long)) & (4294967295 as long)
+    return ((g.state & (65535 as long)) as f32) / (65535.0 as f32) - (0.5 as f32)
+}
+
+// PEEL's large convex (b3CreateComplexHull): the hull of 32 random
+// directions at the radius, from a fixed seed so it is the same hull every
+// run. The directions are drawn and normalised in single precision as the
+// reference draws them, so the points are the reference's points.
+create_complex_hull(radius: float) -> *HullData {
+    count = 32
+    block = basics.alloc(count * sizeof(Vec3))
+    points = block as Vec3[]
+    g = PeelRandom { state: 42 as long }
+    r = radius as f32
+    i = 0
+    while i < count {
+        // A point of the cube kept when it falls in the ball of radius one half.
+        p = math.vec3f_zero()
+        while true {
+            p.x = peel_random_float(&g)
+            p.y = peel_random_float(&g)
+            p.z = peel_random_float(&g)
+            if math.dotf(p, p) <= (0.25 as f32) { break }
+        }
+        points[i] = math.widen3(math.mul_svf(r, math.normalizef(p)))
+        i = i + 1
+    }
+    h = create_hull(points, count, count)
+    basics.free_bytes(block, count * sizeof(Vec3))
+    return h
+}
+
 // --- the box hull ----------------------------------------------------------------
 
 // The signs of a box corner's coordinates, in the order of the box's

aephysics/manifold/module.ae

@@ -450,7 +450,7 @@ reduce_manifold_points(m: *LocalManifold, capacity: int, points: LocalManifoldPo
     // 1: the extreme touching point.
     best_index = NULL_INDEX
     best_score = 0.0 - math.MAX_FLOAT
-    search = arbitrary_perp(normal)
+    search = math.arbitrary_perp(normal)
     i = 0
     while i < count {
         if points[i].separation <= SPECULATIVE_DISTANCE {
@@ -541,22 +541,6 @@ reduce_manifold_points(m: *LocalManifold, capacity: int, points: LocalManifoldPo
     }
 }
 
-// A fixed perpendicular to a unit vector, the reference's, so the
-// reduction's tie-break is the same on every platform.
-arbitrary_perp(v: Vec3) -> Vec3 {
-    a = 0.67
-    b = 0.0 - 0.42
-    p = math.vec3_zero()
-    if v.x < 0.0 - 0.5 || 0.5 < v.x {
-        p = math.vec3(a * v.y + b * v.z, 0.0 - a * v.x, 0.0 - b * v.x)
-    } else if v.y < 0.0 - 0.5 || 0.5 < v.y {
-        p = math.vec3(a * v.y, 0.0 - a * v.x + b * v.z, 0.0 - b * v.y)
-    } else {
-        p = math.vec3(a * v.z, b * v.z, 0.0 - a * v.x - b * v.y)
-    }
-    return math.normalize(p)
-}
-
 // --- spheres and capsules ----------------------------------------------------------------
 
 one_point(m: *LocalManifold, normal: Vec3, point: Vec3, separation: float, pair: FeaturePair) {

aephysics/math/module.ae

@@ -41,7 +41,7 @@ exports (
     min_int, max_int, clamp_int, abs_float, min_float, max_float, clamp_float, lerp_float,
     sqrt_float, atan2_approx, compute_cos_sin, sin_approx, cos_approx, unwind_angle,
     add, sub, mul, neg, dot, length, length_squared, distance, distance_squared,
-    normalize, get_length_and_normalize, perp, is_normalized, mul_add, mul_sub, mul_sv,
+    normalize, get_length_and_normalize, perp, arbitrary_perp, is_normalized, mul_add, mul_sub, mul_sv,
     cross, lerp, blend2, blend3, skew, delta_quat_to_rotation, integrate_rotation, modified_cross, abs_vec3, sign, min_vec3, max_vec3, clamp_vec3, safe_scale, scalar_triple_product,
     is_normalized_quat, rotate_vector, inv_rotate_vector, dot_quat, mul_quat, inv_mul_quat,
     conjugate, negate_quat, normalize_quat, make_quat_from_axis_angle, get_axis_angle,
@@ -348,6 +348,24 @@ perp(a: Vec3) -> Vec3 {
     return normalize(Vec3 { x: 0.0, y: a.z, z: 0.0 - a.y })
 }
 
+// A fixed perpendicular to a unit vector (b3ArbitraryPerp): no component
+// of the vector goes into its own slot, so the dot product cancels. The
+// manifold reduction searches along it, so its tie-break is the same on
+// every platform.
+arbitrary_perp(v: Vec3) -> Vec3 {
+    a = 0.67
+    b = 0.0 - 0.42
+    p = vec3_zero()
+    if v.x < 0.0 - 0.5 || 0.5 < v.x {
+        p = vec3(a * v.y + b * v.z, 0.0 - a * v.x, 0.0 - b * v.x)
+    } else if v.y < 0.0 - 0.5 || 0.5 < v.y {
+        p = vec3(a * v.y, 0.0 - a * v.x + b * v.z, 0.0 - b * v.y)
+    } else {
+        p = vec3(a * v.z, b * v.z, 0.0 - a * v.x - b * v.y)
+    }
+    return normalize(p)
+}
+
 is_normalized(a: Vec3) -> bool {
     aa = dot(a, a)
     return abs_float(1.0 - aa) < 100.0 * EPSILON

aephysics/native/aephysics_inline.h

@@ -54,6 +54,12 @@ static inline long long aephysics_load_word(const void *p, long long index)
     return word;
 }
 
+// The byte at `index` into `p`, unsigned (a hashed block's tail).
+static inline int aephysics_load_u8(const void *p, long long index)
+{
+    return ((const unsigned char *)p)[index];
+}
+
 // The unsigned 16-bit value at `index` halves into `p` (a height field's
 // quanta), one load where two byte reads and a shift were the field's
 // queries' extra cost.

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_u16
+    name_copy, name_text, name_free, prefetch, load_word, load_u8, load_u16
 )
 
 // The most workers a world can have: every module's scratch is a block
@@ -38,6 +38,7 @@ extern aephysics_semaphore_wait(s: ptr)
 extern aephysics_semaphore_signal(s: ptr, count: int)
 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_name_copy(s: string) -> ptr
 extern aephysics_name_text(name: ptr) -> string
@@ -50,6 +51,9 @@ prefetch(p: ptr) { aephysics_prefetch(p) }
 // written as (C's aliasing rule).
 load_word(p: ptr, index: long) -> long { return aephysics_load_word(p, index) }
 
+// The byte at `index` into `p`, unsigned.
+load_u8(p: ptr, index: long) -> int { return aephysics_load_u8(p, index) }
+
 // The unsigned 16-bit value at `index` halves into `p` (little-endian, as
 // every platform aephysics runs on).
 load_u16(p: ptr, index: long) -> int { return aephysics_load_u16(p, index) }

bench/RESULTS.md

@@ -1188,6 +1188,57 @@ at 9 m/s is a step's 15 cm in before its contact exists), the farthest
 body 6.8991 m against 6.8992, at rest 0.8-4.7 mm against our 0.6-5.0,
 every body asleep at the end in both.
 
+## suite (the reference's benchmark suite)
+
+The rest of the reference's `benchmark/main.c` (#141): its nine other
+scenes, built by its own functions (`shared/benchmarks.c`, `e77352c`, the
+release sizes) in `bench/suite_box3d.c` and by their port in
+`bench/suite.ae`, each run as the reference's runner runs it: the scene's
+own step before every step, the first step untimed. One thread, pinned to
+a core, the better of two runs, ms per step:
+
+| scene | aephysics | Box3D `e77352c`, AVX2 | | height sum, contacts (ours / reference) |
+|---|---|---|---|---|
+| convex_pile: 5,120 of PEEL's 32-point convexes | 11.89 | 9.84 | 1.21x | 37,935 / 37,950; 51,767 / 51,706 |
+| junkyard: 10,584 rocks swept by a kinematic cylinder | 22.69 | 21.80 | 1.04x | 85,925 / 90,716; 195,896 / 193,126 |
+| large_world: a million static boxes, 100 spheres dropped | 0.0109 | 0.009 | 1.21x | 76.189 / 76.189; 380 / 380 |
+| sleep: ten pyramids, one woken each step | 2.18 | 1.80 | 1.21x | 14,322.8 / 14,322.8; 5,890 / 5,890 |
+| spinner: 1,500 spheres churned in a cylinder mesh | 3.46 | 3.17 | 1.09x | 19,206 / 19,059; 7,496 / 7,564 |
+| trees25 | 1.45 | 1.25 | 1.16x | 43.99 / 44.48; 1,112 / 1,103 |
+| trees50 | 0.543 | 0.513 | 1.06x | 45.59 / 45.60; 1,106 / 1,100 |
+| trees100 | 0.310 | 0.299 | 1.04x | 45.20 / 45.33; 1,100 / 1,100 |
+| washer: 8,000 cubes in a kinematic drum | 16.87 | 15.14 | 1.11x | 91,422 / 85,532; 41,649 / 43,736 |
+
+Where nothing is chaotic the two agree: the large world and the sleep
+scene to the printed digit, the spinner's motor angle after 800 steps to
+five (-0.157365 against -0.157361). The piles that tumble (the junkyard,
+the washer, the spinner's spheres) part at the rounding of the first
+collisions and end in different heaps of the same size. The reference's
+trees25 warns that a mesh query filled its 256-triangle buffer; ours has
+no such cap, hence its nine more contacts.
+
+The SAT runs (`shared/sat_benchmark.c`): 1,000 pairs of each hull pairing
+placed half the speculative distance apart, the best of four runs of
+twenty repeats, microseconds a query, ours / the reference's:
+
+| pair | GJK cold | GJK warm | SAT | SAT, no inscribed-sphere bound |
+|---|---|---|---|---|
+| complex/complex | 0.654 / 0.618 | 0.212 / 0.196 | 1.78 / 0.954 | 16.3 / 6.08 |
+| complex/rock | 0.433 / 0.411 | 0.139 / 0.129 | 1.02 / 0.645 | 5.63 / 2.76 |
+| complex/cylinder | 1.13 / 1.09 | 0.380 / 0.361 | 2.78 / 1.07 | 14.9 / 6.00 |
+| rock/rock | 0.240 / 0.225 | 0.080 / 0.073 | 0.560 / 0.379 | 1.99 / 1.18 |
+| rock/cylinder | 0.944 / 0.842 | 0.296 / 0.273 | 1.58 / 0.738 | 5.18 / 3.02 |
+| cylinder/cylinder | 1.45 / 1.41 | 0.485 / 0.474 | 4.77 / 1.41 | 13.3 / 7.48 |
+
+GJK is within 3-12% of the reference's. The separating axis test is
+1.5-3.4x: the reference's is the wide one of `convex_manifold_wide.inl`
+(the hull's vertices as eight-lane rows, the edge candidates filtered and
+sorted before any is tested), ours the four-lane edge pass of an earlier
+reference. That is the convex pile's gap, and #150. Across the scenes the
+pairs phase is 1.4-1.6x the reference's (2.59 against 1.63 ms on the
+convex pile): the broad phase's boxes are double, where the reference
+keeps float boxes rounded outward even in its large-world build (#151).
+
 ## Build flags (inlining)
 
 The generated C of every Aether function is a plain `static` function,