lib/linear/src/test.zig
daab053ee43316e1809a84551d573ddd1e5bf3d2
1 const std = @import("std");
2 const hypothesis = @import("hypothesis");
3 const linear = @import("root.zig");
4
5 const Allocator = std.mem.Allocator;
6 const Affine = linear.Affine;
7 const Mat4 = linear.Mat4;
8 const Quat = linear.Quat;
9 const Vec2 = linear.Vec2;
10 const Vec3 = linear.Vec3;
11 const Vec4 = linear.Vec4;
12
13 test "linear package namespace" {
14 std.testing.refAllDecls(linear);
15 }
16
17 test {
18 _ = @import("affine.zig");
19 _ = @import("matrix.zig");
20 _ = @import("quaternion.zig");
21 _ = @import("vector.zig");
22 }
23
24 test "dot rounds (1 + 2^-12)^2 to 1 + 2^-11 before the sum where a fused multiply-add keeps 2^-24" {
25 var a: f32 = 1.0 + 0x1p-12;
26 var c: f32 = -(1.0 + 0x1p-11);
27 std.mem.doNotOptimizeAway(&a);
28 std.mem.doNotOptimizeAway(&c);
29 try std.testing.expectEqual(@as(f32, 0), Vec2.init(a, 1).dot(Vec2.init(a, c)));
30 try std.testing.expectEqual(@as(f32, 0), Vec3.init(a, 1, 0).dot(Vec3.init(a, c, 0)));
31 try std.testing.expectEqual(@as(f32, 0), Vec4.init(a, 1, 0, 0).dot(Vec4.init(a, c, 0, 0)));
32 try std.testing.expectEqual(@as(f32, 0), Quat.init(a, 1, 0, 0).dot(Quat.init(a, c, 0, 0)));
33 }
34
35 /// Folds the raw bits of every correctly rounded operation over a fixed input
36 /// sequence. Transcendental functions stay out, because a linked libc may
37 /// supply its own `sinf` and `tanf`.
38 fn arithmeticDigest() u64 {
39 var state: u64 = 0x9e37_79b9_7f4a_7c15;
40 var digest = std.hash.Wyhash.init(0);
41 for (0..256) |_| {
42 var lanes: [16]f32 = undefined;
43 for (&lanes) |*lane| {
44 state ^= state << 13;
45 state ^= state >> 7;
46 state ^= state << 17;
47 const unit: f32 = @floatFromInt(state >> 40);
48 lane.* = unit / @as(f32, 1 << 24) * 8 - 4;
49 }
50 const a = Vec3.init(lanes[0], lanes[1], lanes[2]);
51 const b = Vec3.init(lanes[3], lanes[4], lanes[5]);
52 const q = Quat.init(lanes[6], lanes[7], lanes[8], lanes[9]).normalized(0).?;
53 const r = Quat.init(lanes[10], lanes[11], lanes[12], lanes[13]).normalized(0).?;
54 const m = Mat4.fromCols(
55 Vec4.fromVec3(a, lanes[14]),
56 Vec4.fromVec3(b, lanes[15]),
57 Vec4.fromVec3(a.cross(b), 0),
58 Vec4.init(lanes[15], lanes[14], lanes[0], 1),
59 );
60 const transform = Affine.fromScaleRotationTranslation(b.abs().add(Vec3.splat(0.5)), q, a);
61 const results = .{
62 a.add(b),
63 a.sub(b).mul(b),
64 a.dot(b),
65 a.cross(b),
66 a.normalized(0).?,
67 a.clampLength(1),
68 a.lerp(b, lanes[14]),
69 q.mul(r),
70 q.rotate(a),
71 q.slerp(r, 0.25).dot(q),
72 linear.Mat3.fromQuat(q).mulVec(b),
73 m.mul(m.transpose()),
74 m.inverse() orelse Mat4.identity,
75 m.determinant(),
76 transform.transformPoint(b),
77 transform.inverse().?.transformPoint(b),
78 transform.normalMatrix().?.mulVec(a),
79 Mat4.lookAt(a, b, Vec3.init(0, 1, 0)) orelse Mat4.identity,
80 };
81 inline for (results) |result| digest.update(std.mem.asBytes(&result));
82 }
83 return digest.final();
84 }
85
86 /// Recorded from x86-64 in Debug and ReleaseFast and from aarch64 under QEMU in
87 /// Debug, ReleaseSafe, and ReleaseFast. A change to any evaluation order
88 /// changes it, and so changes consumers' results.
89 const recorded_arithmetic_digest: u64 = 17181897460903490440;
90
91 test "arithmetic matches the digest recorded on x86-64 and aarch64" {
92 try std.testing.expectEqual(recorded_arithmetic_digest, arithmeticDigest());
93 }
94
95 fn settings() hypothesis.Settings {
96 return hypothesis.Settings.quick()
97 .withSeed(0x4c49_4e45_4152_0001)
98 .withSeedFromEnv()
99 .withDatabase("zig-out/hypothesis-failures/linear");
100 }
101
102 fn drawLane(data: *hypothesis.ConjectureData, bound: f64) !f32 {
103 return @floatCast(try data.drawFloat(-bound, bound));
104 }
105
106 fn drawVec3(data: *hypothesis.ConjectureData, bound: f64) !Vec3 {
107 return .{
108 .x = try drawLane(data, bound),
109 .y = try drawLane(data, bound),
110 .z = try drawLane(data, bound),
111 };
112 }
113
114 /// A unit quaternion from an axis drawn away from zero and any angle.
115 fn drawRotation(data: *hypothesis.ConjectureData) !Quat {
116 const axis = (try drawVec3(data, 1)).normalized(0.1) orelse Vec3.init(0, 1, 0);
117 return Quat.fromAxisAngle(axis, try drawLane(data, std.math.pi));
118 }
119
120 /// Scale lanes bounded away from zero, so the transform stays well conditioned.
121 fn drawScale(data: *hypothesis.ConjectureData) !Vec3 {
122 return .{
123 .x = @floatCast(try data.drawFloat(0.25, 4)),
124 .y = @floatCast(try data.drawFloat(0.25, 4)),
125 .z = @floatCast(try data.drawFloat(0.25, 4)),
126 };
127 }
128
129 fn expectClose(expected: Vec3, actual: Vec3, tolerance: f32) !void {
130 try std.testing.expectApproxEqAbs(expected.x, actual.x, tolerance);
131 try std.testing.expectApproxEqAbs(expected.y, actual.y, tolerance);
132 try std.testing.expectApproxEqAbs(expected.z, actual.z, tolerance);
133 }
134
135 const RotationPreservesGeometry = struct {
136 pub fn property(data: *hypothesis.ConjectureData, _: Allocator) !void {
137 const q = try drawRotation(data);
138 const a = try drawVec3(data, 10);
139 const b = try drawVec3(data, 10);
140 const ra = q.rotate(a);
141 const rb = q.rotate(b);
142 try std.testing.expectApproxEqAbs(a.length(), ra.length(), 1e-4);
143 try std.testing.expectApproxEqAbs(a.dot(b), ra.dot(rb), 1e-3);
144 try expectClose(ra.cross(rb), q.rotate(a.cross(b)), 1e-3);
145 try expectClose(a, q.conjugate().rotate(ra), 1e-4);
146 }
147 };
148
149 const AffineInverseRoundTrip = struct {
150 pub fn property(data: *hypothesis.ConjectureData, _: Allocator) !void {
151 const transform = Affine.fromScaleRotationTranslation(
152 try drawScale(data),
153 try drawRotation(data),
154 try drawVec3(data, 10),
155 );
156 const p = try drawVec3(data, 10);
157 try expectClose(p, transform.inverse().?.transformPoint(transform.transformPoint(p)), 1e-3);
158 const product = transform.toMat4().mul(transform.toMat4().inverse().?);
159 for (product.cols, Mat4.identity.cols) |actual, expected| {
160 try expectClose(expected.xyz(), actual.xyz(), 1e-4);
161 try std.testing.expectApproxEqAbs(expected.w, actual.w, 1e-4);
162 }
163 }
164 };
165
166 const SlerpStaysUnit = struct {
167 pub fn property(data: *hypothesis.ConjectureData, _: Allocator) !void {
168 const a = try drawRotation(data);
169 const b = try drawRotation(data);
170 const t: f32 = @floatCast(try data.drawFloat(0, 1));
171 try std.testing.expectApproxEqAbs(@as(f32, 1), a.slerp(b, t).length(), 1e-5);
172 }
173 };
174
175 test "property: rotation preserves length, dot, and cross" {
176 try hypothesis.checkNamed(RotationPreservesGeometry, "linear-rotation-isometry", settings());
177 }
178
179 test "property: affine and homogeneous inverses round-trip" {
180 try hypothesis.checkNamed(AffineInverseRoundTrip, "linear-affine-inverse", settings());
181 }
182
183 test "property: slerp returns unit quaternions" {
184 try hypothesis.checkNamed(SlerpStaysUnit, "linear-slerp-unit", settings());
185 }