diff --git a/src/sgl/avx.cppm b/src/sgl/avx.cppm index 9fe0e27..f4ee7a2 100644 --- a/src/sgl/avx.cppm +++ b/src/sgl/avx.cppm @@ -251,45 +251,37 @@ inline __m128 load(const mat2& m) noexcept { return _mm_load_ps(&m.cols[0].x); } -} // namespace sgl::detail - -export namespace sgl { - -/* non-constexpr but simd and branchless */ -template bool nearly_equal(const Vec& lhs, const Vec& rhs, const float abs_tol = fp32_abs_tol, const float rel_tol = fp32_rel_tol) noexcept { - const __m128 a = detail::load(lhs); - const __m128 b = detail::load(rhs); - - /* abs(a-b) */ - const __m128 sign_mask = _mm_set1_ps(-0.0f); - - /* Exact equality (handles +/-inf) */ - const __m128 exact = _mm_cmpeq_ps(a, b); - - /* absolute check */ - const __m128 diff = _mm_andnot_ps(sign_mask, _mm_sub_ps(a, b)); - const __m128 abs_ok = _mm_cmple_ps(diff, _mm_set1_ps(abs_tol)); - - /* relative check */ - const __m128 abs_a = _mm_andnot_ps(sign_mask, a); - const __m128 abs_b = _mm_andnot_ps(sign_mask, b); - const __m128 largest = _mm_max_ps(abs_a, abs_b); - const __m128 rel_ok = _mm_cmple_ps(diff, _mm_mul_ps(largest, _mm_set1_ps(rel_tol))); +/* Per-axis slab interval of a ray. An axis the ray is parallel to has an infinite reciprocal, and + * (min - o) * inf is NaN when the origin lies on the slab: minps / maxps would resolve that by + * operand order and NEON by propagating it. Such an axis is a containment test instead: unbounded + * when the origin is within the closed slab, empty otherwise. */ +struct slab_interval_result { + __m128 t_near; + __m128 t_far; +}; - const __m128 ok = _mm_or_ps(exact, _mm_or_ps(abs_ok, rel_ok)); +inline slab_interval_result slab_interval(const __m128 orig, const __m128 inv_dir, const __m128 bmin, const __m128 bmax) noexcept { + const __m128 inf{_mm_set1_ps(std::numeric_limits::infinity())}; + const __m128 neg_inf{_mm_set1_ps(-std::numeric_limits::infinity())}; + const __m128 t1{_mm_mul_ps(_mm_sub_ps(bmin, orig), inv_dir)}; + const __m128 t2{_mm_mul_ps(_mm_sub_ps(bmax, orig), inv_dir)}; + const __m128 parallel{_mm_cmpeq_ps(_mm_andnot_ps(_mm_set1_ps(-0.0f), inv_dir), inf)}; + const __m128 inside{_mm_and_ps(_mm_cmple_ps(bmin, orig), _mm_cmple_ps(orig, bmax))}; + return {_mm_blendv_ps(_mm_min_ps(t1, t2), _mm_blendv_ps(inf, neg_inf, inside), parallel), + _mm_blendv_ps(_mm_max_ps(t1, t2), _mm_blendv_ps(neg_inf, inf, inside), parallel)}; +} - constexpr auto mask{detail::lane_mask()}; - return (_mm_movemask_ps(ok) & mask) == mask; +/* float -> int32 saturating to the int32 range, as NEON's vcvt does; cvtps alone gives INT_MIN for + * anything out of range. NaN gives INT_MIN (maxps returns its second operand), on NEON too. */ +inline __m128i saturating_cvt(const __m128 f) noexcept { + const __m128i clamped_low{_mm_cvtps_epi32(_mm_max_ps(f, _mm_set1_ps(-2147483648.0f)))}; + const __m128 over{_mm_cmpge_ps(f, _mm_set1_ps(2147483648.0f))}; + return _mm_blendv_epi8(clamped_low, _mm_set1_epi32(std::numeric_limits::max()), _mm_castps_si128(over)); } -template bool nearly_zero(const Vec& val, const float tol = fp32_abs_tol) noexcept { - const __m128 sign_mask = _mm_set1_ps(-0.0f); - const __m128 abs_val = _mm_andnot_ps(sign_mask, detail::load(val)); - const __m128 ok = _mm_cmplt_ps(abs_val, _mm_set1_ps(tol)); +} // namespace sgl::detail - constexpr auto mask{detail::lane_mask()}; - return (_mm_movemask_ps(ok) & mask) == mask; -} +export namespace sgl { template bool any_nearly_zero(const Vec& val, const float tol = fp32_abs_tol) noexcept { const __m128 sign_mask = _mm_set1_ps(-0.0f); @@ -596,21 +588,6 @@ template bool is_perpendicular(const Vec& a, const Vec& b, float to return std::abs(d) <= tol; } -template bool is_parallel(const Vec& a, const Vec& b, float tol = fp32_rel_tol) noexcept { - constexpr auto mask{detail::dp_scalar_mask()}; - - /* Cauchy-Schwarz equality - * dot(a,b)^2 ~= dot(a,a) * dot(b,b) */ - const auto va = detail::load(a); - const auto vb = detail::load(b); - - const auto ab{_mm_cvtss_f32(_mm_dp_ps(va, vb, mask))}; - const auto aa{_mm_cvtss_f32(_mm_dp_ps(va, va, mask))}; - const auto bb{_mm_cvtss_f32(_mm_dp_ps(vb, vb, mask))}; - - return std::abs(ab * ab - aa * bb) <= tol * aa * bb; -} - template bool all_greater_than(const Vec& a, const Vec& b) noexcept { constexpr auto mask{detail::lane_mask()}; return (_mm_movemask_ps(_mm_cmpgt_ps(detail::load(a), detail::load(b))) & mask) == mask; @@ -626,11 +603,6 @@ template bool any_greater_than(const Vec& a, const Vec& b) noexcept return static_cast(_mm_movemask_ps(_mm_cmpgt_ps(detail::load(a), detail::load(b))) & mask); } -/* AABB containment: point inside box */ -inline bool contains(const box3d& box, const vec3& point) { - return all_greater_than(point, box.min) && all_less_than(point, box.max); -} - /* vec2: returns scalar (the z-component of the implied 3D cross) */ inline float cross(const vec2& a, const vec2& b) noexcept { const __m128 va = detail::load(a); @@ -662,13 +634,23 @@ inline vec3 cross(const vec3& a, const vec3& b) noexcept { return out; } -/* Specializations, faster */ -template <> inline bool is_parallel(const vec2& a, const vec2& b, float tol) noexcept { - return std::abs(cross(a, b)) <= tol; -} - -template <> inline bool is_parallel(const vec3& a, const vec3& b, float tol) noexcept { - return nearly_zero(cross(a, b), tol); +/* Whether sin^2 of the angle between a and b is within tol: |a x b|^2 <= tol |a|^2 |b|^2, which + * does not depend on the lengths. vec4 has no cross product and uses Lagrange's identity + * |a|^2 |b|^2 - (a.b)^2 instead. A zero vector is parallel to everything. */ +template bool is_parallel(const Vec& a, const Vec& b, const float tol = fp32_rel_tol) noexcept { + const float aa{dot(a, a)}; + const float bb{dot(b, b)}; + float cross_sq{}; + if constexpr (std::is_same_v) { + const float c{cross(a, b)}; + cross_sq = c * c; + } else if constexpr (std::is_same_v) { + cross_sq = length_squared(cross(a, b)); + } else { + const float ab{dot(a, b)}; + cross_sq = std::abs(aa * bb - ab * ab); + } + return cross_sq <= tol * aa * bb; } /* vec2: {-y, x}, a 90 degree counter-clockwise rotation */ @@ -704,21 +686,12 @@ inline vec3 perpendicular(const vec3& v) noexcept { return cross(v, vec3{0.0f, 0.0f, 1.0f}); } +/* Kahan's formula, 2 atan2(|a |b| - b |a||, |a |b| + b |a||): accurate at every angle, where acos + * of the cosine loses the small ones. A zero vector gives 0. */ template float angle_between(const Vec& a, const Vec& b) noexcept { - const auto va = detail::load(a); - const auto vb = detail::load(b); - - /* cos(theta) = dot(a,b) / sqrt(dot(a,a) * dot(b,b)) */ - const auto ab = _mm_dp_ps(va, vb, detail::dp_scalar_mask()); - const auto aa = _mm_dp_ps(va, va, detail::dp_scalar_mask()); - const auto bb = _mm_dp_ps(vb, vb, detail::dp_scalar_mask()); - - const auto cos_theta = _mm_div_ss(ab, _mm_sqrt_ss(_mm_mul_ss(aa, bb))); - - /* clamp to [-1,1], guarding against fp overshoot before acos */ - const auto ct{_mm_cvtss_f32(_mm_min_ss(_mm_max_ss(cos_theta, _mm_set_ss(-1.0f)), _mm_set_ss(1.0f)))}; - - return std::acos(ct); + const Vec u{a * length(b)}; + const Vec v{b * length(a)}; + return 2.0f * std::atan2(length(u - v), length(u + v)); } template bool is_valid(const Box& b) noexcept { @@ -820,6 +793,15 @@ template bool contains_point(const Box& b, const detail::box_vec_ return (ok & mask) == mask; } +/* Closed, like contains_point: a point on a face is inside. */ +inline bool contains(const box2d& box, const vec2& point) noexcept { + return contains_point(box, point); +} + +inline bool contains(const box3d& box, const vec3& point) noexcept { + return contains_point(box, point); +} + template Box intersection(const Box& a, const Box& b) noexcept { __m128 a_lo, a_hi, b_lo, b_hi; detail::load_min_max(a, a_lo, a_hi); @@ -1526,32 +1508,22 @@ inline mat4 to_mat4(const quat& q, const vec3& translation) noexcept { } inline bool nearly_equal(const quat& a, const quat& b, const float abs_tol = fp32_abs_tol, const float rel_tol = fp32_rel_tol) noexcept { - /* quaternions q and -q represent the same rotation,c heck both orientations */ - const __m128 av = detail::load(a); - const __m128 bv = detail::load(b); - const __m128 neg_bv = _mm_xor_ps(bv, _mm_set1_ps(-0.0f)); - - /* |a - b| */ - const __m128 diff1 = _mm_andnot_ps(_mm_set1_ps(-0.0f), _mm_sub_ps(av, bv)); - /* |a + b| */ - const __m128 diff2 = _mm_andnot_ps(_mm_set1_ps(-0.0f), _mm_sub_ps(av, neg_bv)); - /* min(|a-b|, |a+b|) */ - const __m128 diff = _mm_min_ps(diff1, diff2); - - /* Combined tolerance: abs_tol + rel_tol * max(|a|, |b|) */ - const __m128 abs_a = _mm_andnot_ps(_mm_set1_ps(-0.0f), av); - const __m128 abs_b = _mm_andnot_ps(_mm_set1_ps(-0.0f), bv); - const __m128 tol = _mm_add_ps(_mm_set1_ps(abs_tol), _mm_mul_ps(_mm_set1_ps(rel_tol), _mm_max_ps(abs_a, abs_b))); - - const __m128 cmp = _mm_cmple_ps(diff, tol); - return _mm_movemask_ps(cmp) == 0xF; -} - -/* Returns the shortest rotation angle (radians) to go from orientation a to b: - * angle = 2 * acos(|dot(a, b)|) */ + /* q and -q are the same rotation: equal when every lane of a - b, or every lane of a + b, is within + * tolerance. Not a per-lane choice, which would accept a mix of the two. */ + const __m128 sign{_mm_set1_ps(-0.0f)}; + const __m128 av{detail::load(a)}; + const __m128 bv{detail::load(b)}; + const __m128 tol{_mm_add_ps(_mm_set1_ps(abs_tol), _mm_mul_ps(_mm_set1_ps(rel_tol), _mm_max_ps(_mm_andnot_ps(sign, av), _mm_andnot_ps(sign, bv))))}; + const __m128 same{_mm_cmple_ps(_mm_andnot_ps(sign, _mm_sub_ps(av, bv)), tol)}; + const __m128 flipped{_mm_cmple_ps(_mm_andnot_ps(sign, _mm_add_ps(av, bv)), tol)}; + return _mm_movemask_ps(same) == 0xF || _mm_movemask_ps(flipped) == 0xF; +} + +/* The shortest rotation angle (radians, in [0, pi]) from orientation a to b: 2 atan2(|xyz|, |w|) of + * conj(a) * b, accurate where 2 acos(|dot(a, b)|) loses the small angles. */ inline float angle_between(const quat& a, const quat& b) noexcept { - /* note that we clamp for acos safety */ - return 2.0f * std::acos(std::min(std::abs(dot(a, b)), 1.0f)); + const quat d{conjugate(a) * b}; + return 2.0f * std::atan2(std::sqrt(d.x * d.x + d.y * d.y + d.z * d.z), std::abs(d.w)); } /* Returns true when `direction` is within `tolerance` of a cardinal axis. */ @@ -1586,6 +1558,7 @@ template Vec closest_point_on_segment(const Vec& point, con /* t = dot(ap, ab) / dot(ab, ab), clamped to [0, 1] */ auto t{_mm_div_ps(ab_dot_ap, ab_dot_ab)}; + /* maxps returns its second operand for NaN, so a zero-length segment's 0 / 0 gives t = 0: the start. */ t = _mm_max_ps(t, _mm_setzero_ps()); t = _mm_min_ps(t, _mm_set1_ps(1.0f)); @@ -1827,8 +1800,8 @@ inline segment_intersection_2d intersect_segments_2d(const vec2& a0, const vec2& * Computes t_entry and t_exit for all slab pairs simultaneously. * * Uses reciprocal direction to turn division into multiplication. - * Handles infinities correctly when direction component is 0 - * (IEEE 754: 1/0 = +/-inf, then min/max propagate correctly). + * The box is closed: a ray grazing a face hits it. A zero direction component makes that axis a + * containment test (see detail::slab_interval). */ inline ray_box_hit intersect_ray_box(const vec2& origin, const vec2& direction, const box2d& box) noexcept { /* For vec2 we only care about lanes 0 and 1. @@ -1838,12 +1811,7 @@ inline ray_box_hit intersect_ray_box(const vec2& origin, const vec2& direction, __m128 bmin, bmax; detail::load_min_max(box, bmin, bmax); - /* t values for near and far slabs */ - const __m128 t1 = _mm_mul_ps(_mm_sub_ps(bmin, orig), inv_dir); - const __m128 t2 = _mm_mul_ps(_mm_sub_ps(bmax, orig), inv_dir); - - const __m128 t_near = _mm_min_ps(t1, t2); - const __m128 t_far = _mm_max_ps(t1, t2); + const auto [t_near, t_far]{detail::slab_interval(orig, inv_dir, bmin, bmax)}; /* t_entry = max of all near values, t_exit = min of all far values */ /* For vec2: max(t_near.x, t_near.y) and min(t_far.x, t_far.y) */ @@ -1867,11 +1835,7 @@ inline ray_box_hit intersect_ray_box(const vec3& origin, const vec3& direction, const __m128 bmin = _mm_load_ps(&box.min.x); const __m128 bmax = _mm_load_ps(&box.max.x); - const __m128 t1 = _mm_mul_ps(_mm_sub_ps(bmin, orig), inv_dir); - const __m128 t2 = _mm_mul_ps(_mm_sub_ps(bmax, orig), inv_dir); - - const __m128 t_near = _mm_min_ps(t1, t2); - const __m128 t_far = _mm_max_ps(t1, t2); + const auto [t_near, t_far]{detail::slab_interval(orig, inv_dir, bmin, bmax)}; /* Horizontal max of t_near lanes 0,1,2 and horizontal min of t_far lanes 0,1,2 */ /* Shuffle: {y, z, x, _} */ @@ -1968,8 +1932,7 @@ inline ivec2 grid_cell(const vec2& point, const vec2& grid_origin, const float c const __m128 inv_size = _mm_set1_ps(1.0f / cell_size); const __m128 cell_f = _mm_floor_ps(_mm_mul_ps(p, inv_size)); - /* Convert to int. _mm_cvttps_epi32 truncates, but we already floored */ - const __m128i cell_i = _mm_cvtps_epi32(cell_f); /* rounds to nearest, but value is already integer after floor */ + const __m128i cell_i{detail::saturating_cvt(cell_f)}; alignas(16) int cells[4]{}; _mm_store_si128(reinterpret_cast<__m128i*>(cells), cell_i); @@ -1982,7 +1945,7 @@ inline ivec3 grid_cell(const vec3& point, const vec3& grid_origin, const float c const __m128 inv_size = _mm_set1_ps(1.0f / cell_size); const __m128 cell_f = _mm_floor_ps(_mm_mul_ps(p, inv_size)); - const __m128i cell_i = _mm_cvtps_epi32(cell_f); + const __m128i cell_i{detail::saturating_cvt(cell_f)}; alignas(16) int cells[4]{}; _mm_store_si128(reinterpret_cast<__m128i*>(cells), cell_i); @@ -2043,18 +2006,11 @@ inline transform inverse_uniform(const transform& tf) noexcept { return {inv_pos, inv_rot, {inv_s, inv_s, inv_s}}; } -/* Inverse of a non-uniform-scale transform. - * - * inv_scale = { 1/sx, 1/sy, 1/sz } - * inv_rotation = conjugate(rotation) - * inv_position = inv_rotation * (-position / scale) - * - * Note: non-uniform scale + rotation doesn't compose cleanly. - * The inverse is correct for a single transform, but chaining - * non-uniform-scale transforms is not equivalent to multiplying - * these structs (requires a full mat4 for that). - */ +/* Inverse of a transform with uniform scale. With a non-uniform scale the inverse applies the scale + * after the rotation, which a scale-then-rotate transform cannot represent: use inverse(to_mat4(tf)). */ inline transform inverse(const transform& tf) noexcept { + assert( + nearly_equal(tf.scale.x, tf.scale.y) && nearly_equal(tf.scale.x, tf.scale.z) && "inverse(transform) needs a uniform scale; use inverse(to_mat4(tf))"); const __m128 one = _mm_set1_ps(1.0f); const __m128 s = _mm_load_ps(&tf.scale.x); const __m128 inv_s = _mm_div_ps(one, s); diff --git a/src/sgl/neon.cppm b/src/sgl/neon.cppm index b45bfad..821517e 100644 --- a/src/sgl/neon.cppm +++ b/src/sgl/neon.cppm @@ -165,6 +165,12 @@ inline uint32x2_t n_or(uint32x2_t a, uint32x2_t b) noexcept { inline uint32x4_t n_or(uint32x4_t a, uint32x4_t b) noexcept { return vorrq_u32(a, b); } +inline uint32x2_t n_and(uint32x2_t a, uint32x2_t b) noexcept { + return vand_u32(a, b); +} +inline uint32x4_t n_and(uint32x4_t a, uint32x4_t b) noexcept { + return vandq_u32(a, b); +} /* --- Blend (bitwise select) --- */ inline float32x2_t n_bsl(uint32x2_t m, float32x2_t a, float32x2_t b) noexcept { @@ -462,6 +468,31 @@ inline float32x4_t load(const mat2& m) noexcept { return vld1q_f32(&m.cols[0].x); } +/* @see the AVX `slab_interval`: no NaN reaches min/max, so the result matches it lane for lane. */ +template struct slab_interval_result { + R t_near; + R t_far; +}; + +template slab_interval_result slab_interval(const R orig, const R inv_dir, const R bmin, const R bmax) noexcept { + const R inf{n_dup(std::numeric_limits::infinity(), orig)}; + const R neg_inf{n_dup(-std::numeric_limits::infinity(), orig)}; + const R t1{n_mul(n_sub(bmin, orig), inv_dir)}; + const R t2{n_mul(n_sub(bmax, orig), inv_dir)}; + const auto parallel{n_ceq(n_abs(inv_dir), inf)}; + const auto inside{n_and(n_cle(bmin, orig), n_cle(orig, bmax))}; + return {n_bsl(parallel, n_bsl(inside, neg_inf, inf), n_min(t1, t2)), n_bsl(parallel, n_bsl(inside, inf, neg_inf), n_max(t1, t2))}; +} + +/* @see the AVX `saturating_cvt`. vcvt saturates already; NaN would give 0 and is mapped to INT_MIN. */ +inline int32x4_t saturating_cvt(const float32x4_t f) noexcept { + return vbslq_s32(vceqq_f32(f, f), vcvtq_s32_f32(f), vdupq_n_s32(std::numeric_limits::min())); +} + +inline int32x2_t saturating_cvt(const float32x2_t f) noexcept { + return vbsl_s32(vceq_f32(f, f), vcvt_s32_f32(f), vdup_n_s32(std::numeric_limits::min())); +} + } // namespace sgl::detail export namespace sgl { @@ -470,26 +501,6 @@ export namespace sgl { * Comparisons * ============================================================ */ -template bool nearly_equal(const Vec& lhs, const Vec& rhs, const float abs_tol = fp32_abs_tol, const float rel_tol = fp32_rel_tol) noexcept { - const auto a = detail::load(lhs); - const auto b = detail::load(rhs); - const auto diff = detail::n_abs(detail::n_sub(a, b)); - const auto exact = detail::n_ceq(a, b); - const auto abs_ok = detail::n_cle(diff, detail::n_dup(abs_tol, a)); - const auto largest = detail::n_max(detail::n_abs(a), detail::n_abs(b)); - const auto rel_ok = detail::n_cle(diff, detail::n_mul(largest, detail::n_dup(rel_tol, a))); - const auto ok = detail::n_or(exact, detail::n_or(abs_ok, rel_ok)); - constexpr auto mask{detail::lane_mask()}; - return (detail::n_movemask(ok) & mask) == mask; -} - -template bool nearly_zero(const Vec& val, const float tol = fp32_abs_tol) noexcept { - const auto v = detail::load(val); - const auto ok = detail::n_clt(detail::n_abs(v), detail::n_dup(tol, v)); - constexpr auto mask{detail::lane_mask()}; - return (detail::n_movemask(ok) & mask) == mask; -} - template bool any_nearly_zero(const Vec& val, const float tol = fp32_abs_tol) noexcept { const auto v = detail::load(val); const auto ok = detail::n_clt(detail::n_abs(v), detail::n_dup(tol, v)); @@ -771,15 +782,6 @@ template bool is_perpendicular(const Vec& a, const Vec& b, float to return std::abs(detail::dot_scalar(detail::load(a), detail::load(b))) <= tol; } -template bool is_parallel(const Vec& a, const Vec& b, float tol = fp32_rel_tol) noexcept { - const auto va = detail::load(a); - const auto vb = detail::load(b); - const auto ab = detail::dot_scalar(va, vb); - const auto aa = detail::dot_scalar(va, va); - const auto bb = detail::dot_scalar(vb, vb); - return std::abs(ab * ab - aa * bb) <= tol * aa * bb; -} - template bool all_greater_than(const Vec& a, const Vec& b) noexcept { constexpr auto mask{detail::lane_mask()}; return (detail::n_movemask(detail::n_cgt(detail::load(a), detail::load(b))) & mask) == mask; @@ -795,10 +797,6 @@ template bool any_greater_than(const Vec& a, const Vec& b) noexcept return static_cast(detail::n_movemask(detail::n_cgt(detail::load(a), detail::load(b))) & mask); } -inline bool contains(const box3d& box, const vec3& point) { - return all_greater_than(point, box.min) && all_less_than(point, box.max); -} - /* ============================================================ * Cross products * ============================================================ */ @@ -824,11 +822,21 @@ inline vec3 cross(const vec3& a, const vec3& b) noexcept { return out; } -template <> inline bool is_parallel(const vec2& a, const vec2& b, float tol) noexcept { - return std::abs(cross(a, b)) <= tol; -} -template <> inline bool is_parallel(const vec3& a, const vec3& b, float tol) noexcept { - return nearly_zero(cross(a, b), tol); +/* @see the AVX `is_parallel`. */ +template bool is_parallel(const Vec& a, const Vec& b, const float tol = fp32_rel_tol) noexcept { + const float aa{dot(a, a)}; + const float bb{dot(b, b)}; + float cross_sq{}; + if constexpr (std::is_same_v) { + const float c{cross(a, b)}; + cross_sq = c * c; + } else if constexpr (std::is_same_v) { + cross_sq = length_squared(cross(a, b)); + } else { + const float ab{dot(a, b)}; + cross_sq = std::abs(aa * bb - ab * ab); + } + return cross_sq <= tol * aa * bb; } /* vec2 perpendicular: {-y, x}, D-form */ @@ -854,11 +862,11 @@ inline vec3 perpendicular(const vec3& v) noexcept { return cross(v, vec3{0, 0, 1}); } +/* @see the AVX `angle_between`. */ template float angle_between(const Vec& a, const Vec& b) noexcept { - const auto va = detail::load(a); - const auto vb = detail::load(b); - auto ct = detail::dot_scalar(va, vb) / std::sqrt(detail::dot_scalar(va, va) * detail::dot_scalar(vb, vb)); - return std::acos(std::max(-1.0f, std::min(ct, 1.0f))); + const Vec u{a * length(b)}; + const Vec v{b * length(a)}; + return 2.0f * std::atan2(length(u - v), length(u + v)); } /* ============================================================ @@ -954,6 +962,15 @@ template bool contains_point(const Box& b, const detail::box_vec_ return ((detail::n_movemask(detail::n_cle(lo, pv)) & detail::n_movemask(detail::n_cle(pv, hi))) & mask) == mask; } +/* @see the AVX `contains`: closed, like contains_point. */ +inline bool contains(const box2d& box, const vec2& point) noexcept { + return contains_point(box, point); +} + +inline bool contains(const box3d& box, const vec3& point) noexcept { + return contains_point(box, point); +} + template Box intersection(const Box& a, const Box& b) noexcept { detail::box_reg_t a_lo, a_hi, b_lo, b_hi; detail::load_min_max(a, a_lo, a_hi); @@ -1377,15 +1394,20 @@ inline mat4 to_mat4(const quat& q, const vec3& t) noexcept { return to_mat4(to_mat3(q), t); } +/* @see the AVX `nearly_equal(quat)`. */ inline bool nearly_equal(const quat& a, const quat& b, float abs_tol = fp32_abs_tol, float rel_tol = fp32_rel_tol) noexcept { - float32x4_t av = detail::load(a), bv = detail::load(b); - float32x4_t diff = vminq_f32(vabsq_f32(vsubq_f32(av, bv)), vabsq_f32(vaddq_f32(av, bv))); - float32x4_t tol = vaddq_f32(vdupq_n_f32(abs_tol), vmulq_f32(vdupq_n_f32(rel_tol), vmaxq_f32(vabsq_f32(av), vabsq_f32(bv)))); - return detail::n_movemask(vcleq_f32(diff, tol)) == 0xF; + const float32x4_t av{detail::load(a)}; + const float32x4_t bv{detail::load(b)}; + const float32x4_t tol{vaddq_f32(vdupq_n_f32(abs_tol), vmulq_f32(vdupq_n_f32(rel_tol), vmaxq_f32(vabsq_f32(av), vabsq_f32(bv))))}; + const bool same{detail::n_movemask(vcleq_f32(vabsq_f32(vsubq_f32(av, bv)), tol)) == 0xF}; + const bool flipped{detail::n_movemask(vcleq_f32(vabsq_f32(vaddq_f32(av, bv)), tol)) == 0xF}; + return same || flipped; } +/* @see the AVX `angle_between(quat)`. */ inline float angle_between(const quat& a, const quat& b) noexcept { - return 2.0f * std::acos(std::min(std::abs(dot(a, b)), 1.0f)); + const quat d{conjugate(a) * b}; + return 2.0f * std::atan2(std::sqrt(d.x * d.x + d.y * d.y + d.z * d.z), std::abs(d.w)); } inline bool is_axis_aligned(const vec3& direction, float tolerance = fp32_rel_tol) noexcept { @@ -1400,7 +1422,12 @@ template Vec closest_point_on_segment(const Vec& point, con const auto p = detail::load(point), a = detail::load(seg_start), b = detail::load(seg_end); const auto ab = detail::n_sub(b, a), ap = detail::n_sub(p, a); auto t = detail::n_div(detail::dot_broadcast(ab, ap), detail::dot_broadcast(ab, ab)); - t = detail::n_min(detail::n_max(t, detail::n_dup(0.0f, t)), detail::n_dup(1.0f, t)); + /* t > 0 ? t : 0, the order maxps uses on x86, not FMAX, which propagates NaN: a zero-length + * segment's 0 / 0 then gives t = 0, the start, on both backends. */ + const auto zero{detail::n_dup(0.0f, t)}; + const auto one{detail::n_dup(1.0f, t)}; + t = detail::n_bsl(detail::n_cgt(t, zero), t, zero); + t = detail::n_bsl(detail::n_clt(t, one), t, one); Vec out{}; detail::store(detail::n_fma(a, t, ab), out); return out; @@ -1566,11 +1593,11 @@ inline segment_intersection_2d intersect_segments_2d(const vec2& a0, const vec2& } inline ray_box_hit intersect_ray_box(const vec2& origin, const vec2& direction, const box2d& box) noexcept { - float32x2_t orig = detail::load(origin), inv_dir = vdiv_f32(vdup_n_f32(1.0f), detail::load(direction)); + const float32x2_t orig{detail::load(origin)}; + const float32x2_t inv_dir{vdiv_f32(vdup_n_f32(1.0f), detail::load(direction))}; float32x2_t bmin, bmax; detail::load_min_max(box, bmin, bmax); - float32x2_t t1 = vmul_f32(vsub_f32(bmin, orig), inv_dir), t2 = vmul_f32(vsub_f32(bmax, orig), inv_dir); - float32x2_t t_near = vmin_f32(t1, t2), t_far = vmax_f32(t1, t2); + const auto [t_near, t_far]{detail::slab_interval(orig, inv_dir, bmin, bmax)}; const float te{std::max(vget_lane_f32(t_near, 0), vget_lane_f32(t_near, 1))}; const float tx{std::min(vget_lane_f32(t_far, 0), vget_lane_f32(t_far, 1))}; if (te > tx || tx < 0) { @@ -1580,9 +1607,9 @@ inline ray_box_hit intersect_ray_box(const vec2& origin, const vec2& direction, } inline ray_box_hit intersect_ray_box(const vec3& origin, const vec3& direction, const box3d& box) noexcept { - float32x4_t orig = vld1q_f32(&origin.x), inv_dir = vdivq_f32(vdupq_n_f32(1.0f), vld1q_f32(&direction.x)); - float32x4_t t1 = vmulq_f32(vsubq_f32(vld1q_f32(&box.min.x), orig), inv_dir), t2 = vmulq_f32(vsubq_f32(vld1q_f32(&box.max.x), orig), inv_dir); - float32x4_t tn = vminq_f32(t1, t2), tf = vmaxq_f32(t1, t2); + const float32x4_t orig{vld1q_f32(&origin.x)}; + const float32x4_t inv_dir{vdivq_f32(vdupq_n_f32(1.0f), vld1q_f32(&direction.x))}; + const auto [tn, tf]{detail::slab_interval(orig, inv_dir, vld1q_f32(&box.min.x), vld1q_f32(&box.max.x))}; const float te{std::max({vgetq_lane_f32(tn, 0), vgetq_lane_f32(tn, 1), vgetq_lane_f32(tn, 2)})}; const float tx{std::min({vgetq_lane_f32(tf, 0), vgetq_lane_f32(tf, 1), vgetq_lane_f32(tf, 2)})}; if (te > tx || tx < 0) { @@ -1639,13 +1666,13 @@ inline segment_polygon_hit intersect_segment_polygon(const vec2& a, const vec2& inline ivec2 grid_cell(const vec2& point, const vec2& grid_origin, float cell_size) noexcept { float32x2_t p = vmul_n_f32(vsub_f32(detail::load(point), detail::load(grid_origin)), 1.0f / cell_size); - int32x2_t ci = vcvt_s32_f32(vrndm_f32(p)); + const int32x2_t ci{detail::saturating_cvt(vrndm_f32(p))}; return {vget_lane_s32(ci, 0), vget_lane_s32(ci, 1)}; } inline ivec3 grid_cell(const vec3& point, const vec3& grid_origin, float cell_size) noexcept { float32x4_t p = vmulq_n_f32(vsubq_f32(vld1q_f32(&point.x), vld1q_f32(&grid_origin.x)), 1.0f / cell_size); - int32x4_t ci = vcvtq_s32_f32(vrndmq_f32(p)); + const int32x4_t ci{detail::saturating_cvt(vrndmq_f32(p))}; return {vgetq_lane_s32(ci, 0), vgetq_lane_s32(ci, 1), vgetq_lane_s32(ci, 2)}; } @@ -1675,7 +1702,10 @@ inline transform inverse_uniform(const transform& tf) noexcept { return {rotate(ir, {-tf.position.x * is, -tf.position.y * is, -tf.position.z * is}), ir, {is, is, is}}; } +/* @see the AVX `inverse(transform)`: uniform scale only. */ inline transform inverse(const transform& tf) noexcept { + assert( + nearly_equal(tf.scale.x, tf.scale.y) && nearly_equal(tf.scale.x, tf.scale.z) && "inverse(transform) needs a uniform scale; use inverse(to_mat4(tf))"); float32x4_t inv_s = vdivq_f32(vdupq_n_f32(1.0f), vld1q_f32(&tf.scale.x)); vec3 is{}; vst1q_f32(&is.x, inv_s); diff --git a/src/sgl/scalar.cppm b/src/sgl/scalar.cppm index 10896e1..b709a37 100644 --- a/src/sgl/scalar.cppm +++ b/src/sgl/scalar.cppm @@ -45,7 +45,8 @@ constexpr float fp32_rel_tol{128.0f * std::numeric_limits::epsilon()}; /** * Returns true when lhs and rhs agree within an absolute or relative tolerance. - * Bitwise equality is checked first to handle +/-inf and +/-0 correctly. + * Exact equality is checked first, so equal infinities and +0 / -0 compare equal; an infinity is + * otherwise equal to nothing. */ inline SGL_CMATH_CONSTEXPR bool nearly_equal( const float lhs, const float rhs, const float abs_tol = fp32_abs_tol, const float rel_tol = fp32_rel_tol) noexcept { @@ -53,6 +54,11 @@ inline SGL_CMATH_CONSTEXPR bool nearly_equal( return true; } + /* Unequal with an infinity on either side: |inf - x| = inf would pass the relative check against inf * rel_tol. */ + if (std::isinf(lhs) || std::isinf(rhs)) { + return false; + } + const auto diff{std::abs(lhs - rhs)}; /* absolute check handles values near zero */ diff --git a/src/sgl/types.cppm b/src/sgl/types.cppm index 5fe59cf..fb2419b 100644 --- a/src/sgl/types.cppm +++ b/src/sgl/types.cppm @@ -43,14 +43,17 @@ export struct alignas(8) vec2 { float x; ///< @brief x component. float y; ///< @brief y component. + /* Member selection, not (&x)[i]: indexing past a scalar member is UB and not a constant + * expression. A constant index folds to the same load; a runtime one costs a compare or two + * (cmov or branch) where (&x)[i] was one indexed load. */ constexpr float operator[](std::size_t i) const noexcept { assert(i < 2); - return (&x)[i]; + return i == 0 ? x : y; } constexpr float& operator[](std::size_t i) noexcept { assert(i < 2); - return (&x)[i]; + return i == 0 ? x : y; } }; @@ -70,12 +73,26 @@ export struct alignas(16) vec3 { constexpr float operator[](std::size_t i) const noexcept { assert(i < 3); - return (&x)[i]; + switch (i) { + case 0: + return x; + case 1: + return y; + default: + return z; + } } constexpr float& operator[](std::size_t i) noexcept { assert(i < 3); - return (&x)[i]; + switch (i) { + case 0: + return x; + case 1: + return y; + default: + return z; + } } }; @@ -88,12 +105,30 @@ export struct alignas(16) vec4 { constexpr float operator[](std::size_t i) const noexcept { assert(i < 4); - return (&x)[i]; + switch (i) { + case 0: + return x; + case 1: + return y; + case 2: + return z; + default: + return w; + } } constexpr float& operator[](std::size_t i) noexcept { assert(i < 4); - return (&x)[i]; + switch (i) { + case 0: + return x; + case 1: + return y; + case 2: + return z; + default: + return w; + } } }; diff --git a/tests/src/box_test.cpp b/tests/src/box_test.cpp index 837edbf..f94b6f5 100644 --- a/tests/src/box_test.cpp +++ b/tests/src/box_test.cpp @@ -67,6 +67,13 @@ TEST(Box2d, Contains_Box2d) { EXPECT_FALSE(sgl::contains(outer, outside)); } +TEST(Box2d, ContainsIsClosed) { + const sgl::box2d b{.min = {0, 0}, .max = {10, 10}}; + EXPECT_TRUE(sgl::contains(b, sgl::vec2{0, 5})); + EXPECT_TRUE(sgl::contains(b, sgl::vec2{10, 10})); + EXPECT_FALSE(sgl::contains(b, sgl::vec2{-0.001f, 5})); +} + TEST(Box2d, ContainsPoint) { sgl::box2d b{.min = {0, 0}, .max = {10, 10}}; EXPECT_TRUE(sgl::contains_point(b, sgl::vec2{5, 5})); @@ -141,6 +148,14 @@ TEST(Box3d, ContainsPoint_Box3d) { EXPECT_TRUE(sgl::contains(sgl::box3d{.min = {0, 0, 0}, .max = {10, 10, 10}}, sgl::vec3{5, 5, 5})); } +/* Regression: contains was strict while contains_point (and sgl::soa::contains) are closed. */ +TEST(Box3d, ContainsIsClosed) { + const sgl::box3d b{.min = {0, 0, 0}, .max = {10, 10, 10}}; + EXPECT_TRUE(sgl::contains(b, sgl::vec3{0, 5, 5})); + EXPECT_TRUE(sgl::contains(b, sgl::vec3{10, 10, 10})); + EXPECT_FALSE(sgl::contains(b, sgl::vec3{10.001f, 5, 5})); +} + TEST(Box3d, Volume) { EXPECT_FLOAT_EQ(sgl::volume(sgl::box3d{.min = {0, 0, 0}, .max = {2, 3, 4}}), 24.0f); } diff --git a/tests/src/geometry_test.cpp b/tests/src/geometry_test.cpp index 55c8b4b..d6f9618 100644 --- a/tests/src/geometry_test.cpp +++ b/tests/src/geometry_test.cpp @@ -2,7 +2,9 @@ * @file geometry_test.cpp * @brief Unit tests for geometric operations (segments, planes, intersections, grids, transforms). */ +#include #include +#include import sgl; @@ -186,6 +188,52 @@ TEST(Geometry, GridCellRange2d) { EXPECT_EQ(r.max.x, 3); } +/* Regression: out-of-range cells were INT_MIN on x86 and INT_MAX (saturated) on NEON. */ +TEST(Geometry, GridCellSaturates) { + constexpr std::int32_t lo{std::numeric_limits::min()}; + constexpr std::int32_t hi{std::numeric_limits::max()}; + const auto c2{sgl::grid_cell(sgl::vec2{3e9f, -3e9f}, sgl::vec2{0, 0}, 1.0f)}; + EXPECT_EQ(c2.x, hi); + EXPECT_EQ(c2.y, lo); + const auto c3{sgl::grid_cell(sgl::vec3{-1e20f, 2.5f, 1e20f}, sgl::vec3{0, 0, 0}, 1.0f)}; + EXPECT_EQ(c3.x, lo); + EXPECT_EQ(c3.y, 2); + EXPECT_EQ(c3.z, hi); + const auto nan_cell{sgl::grid_cell(sgl::vec2{std::numeric_limits::quiet_NaN(), 1.5f}, sgl::vec2{0, 0}, 1.0f)}; + EXPECT_EQ(nan_cell.x, lo) << "NaN maps to INT_MIN"; + EXPECT_EQ(nan_cell.y, 1); +} + +/* Regression: the backends disagreed on these (NaN handling of minps/maxps vs FMIN/FMAX). */ +TEST(Geometry, RayBoxGrazingEdgesHitOnEveryBackend) { + const sgl::box2d square{.min = {0, 0}, .max = {1, 1}}; + for (const float y : {0.0f, 1.0f}) { + const auto h{sgl::intersect_ray_box(sgl::vec2{-1, y}, sgl::vec2{1, 0}, square)}; + EXPECT_TRUE(h.hit) << "y=" << y; + EXPECT_EQ(h.t_entry, 1.0f) << "y=" << y; + EXPECT_EQ(h.t_exit, 2.0f) << "y=" << y; + } + EXPECT_FALSE(sgl::intersect_ray_box(sgl::vec2{-1, -0.0001f}, sgl::vec2{1, 0}, square).hit); + EXPECT_FALSE(sgl::intersect_ray_box(sgl::vec2{-1, 1.0001f}, sgl::vec2{1, -0.0f}, square).hit); + + const sgl::box3d cube{.min = {0, 0, 0}, .max = {1, 1, 1}}; + const auto edge{sgl::intersect_ray_box(sgl::vec3{-1, 0, 1}, sgl::vec3{1, 0, 0}, cube)}; + EXPECT_TRUE(edge.hit); + EXPECT_EQ(edge.t_entry, 1.0f); + EXPECT_EQ(edge.t_exit, 2.0f); + EXPECT_FALSE(sgl::intersect_ray_box(sgl::vec3{-1, 0, 1.0001f}, sgl::vec3{1, 0, 0}, cube).hit); +} + +TEST(Geometry, ClosestPointOnDegenerateSegmentIsItsStart) { + const auto p3{sgl::closest_point_on_segment(sgl::vec3{5, 5, 5}, sgl::vec3{1, 2, 3}, sgl::vec3{1, 2, 3})}; + EXPECT_EQ(p3.x, 1.0f); + EXPECT_EQ(p3.y, 2.0f); + EXPECT_EQ(p3.z, 3.0f); + const auto p2{sgl::closest_point_on_segment(sgl::vec2{5, 5}, sgl::vec2{1, 2}, sgl::vec2{1, 2})}; + EXPECT_EQ(p2.x, 1.0f); + EXPECT_EQ(p2.y, 2.0f); +} + /* ============================================================ * Transform operations * ============================================================ */ @@ -225,10 +273,24 @@ TEST(Geometry, InverseUniformTransform) { TEST(Geometry, InverseTransform) { sgl::transform tf{}; tf.position = {1, 2, 3}; - tf.scale = {2, 4, 8}; - auto inv = sgl::inverse(tf); + tf.rotation = sgl::quat_from_axis_angle(sgl::vec3{0, 0, 1}, 1.5707963f); + tf.scale = {2, 2, 2}; + const auto inv{sgl::inverse(tf)}; EXPECT_NEAR(inv.scale.x, 0.5f, 1e-5f); - EXPECT_NEAR(inv.scale.y, 0.25f, 1e-5f); + const sgl::vec3 p{1, -2, 0.5f}; + const auto back{sgl::transform_point(inv, sgl::transform_point(tf, p))}; + EXPECT_NEAR(back.x, p.x, 1e-5f); + EXPECT_NEAR(back.y, p.y, 1e-5f); + EXPECT_NEAR(back.z, p.z, 1e-5f); +} + +/* Regression: a non-uniform inverse is not a scale-then-rotate transform, so inverse() returned a + * wrong one; it now asserts uniform scale (release builds keep computing the uniform formula). */ +TEST(Geometry, InverseTransformRequiresUniformScale) { + sgl::transform tf{}; + tf.rotation = sgl::quat_from_axis_angle(sgl::vec3{0, 0, 1}, 1.5707963f); + tf.scale = {2, 1, 1}; + EXPECT_DEBUG_DEATH((void)sgl::inverse(tf), "uniform"); } TEST(Geometry, ComposeTransform) { diff --git a/tests/src/quaternion_test.cpp b/tests/src/quaternion_test.cpp index 069a763..f6e1b5a 100644 --- a/tests/src/quaternion_test.cpp +++ b/tests/src/quaternion_test.cpp @@ -138,6 +138,14 @@ TEST(Quat, NearlyEqual) { EXPECT_TRUE(sgl::nearly_equal(a, c)); /* same rotation, opposite sign */ } +/* Regression: the sign was chosen per lane, so two different rotations could compare equal. */ +TEST(Quat, NearlyEqualDoesNotMixSigns) { + EXPECT_FALSE(sgl::nearly_equal(sgl::quat{0.5f, 0.5f, 0.5f, 0.5f}, sgl::quat{0.5f, -0.5f, 0.5f, -0.5f})); + const auto q{sgl::quat_from_axis_angle(sgl::vec3{1, 2, 3}, 0.7f)}; + EXPECT_TRUE(sgl::nearly_equal(q, sgl::negate(q))); + EXPECT_FALSE(sgl::nearly_equal(q, sgl::quat_from_axis_angle(sgl::vec3{1, 2, 3}, 0.8f))); +} + TEST(Quat, AngleBetween) { sgl::quat q1 = sgl::quat_from_axis_angle(sgl::vec3{0, 0, 1}, 0.0f); sgl::quat q2 = sgl::quat_from_axis_angle(sgl::vec3{0, 0, 1}, std::acos(-1.0f) * 0.5f); /* 90 degrees */ @@ -145,6 +153,18 @@ TEST(Quat, AngleBetween) { EXPECT_NEAR(angle, std::acos(-1.0f) * 0.5f, 1e-4f); } +/* Regression: acos near 1 lost the small angles (1e-3 rad came back as 9.77e-4). */ +TEST(Quat, AngleBetweenSmallAndLarge) { + for (const float angle : {1e-4f, 1e-3f, 0.5f, 3.0f}) { + const auto q{sgl::quat_from_axis_angle(sgl::vec3{0, 0, 1}, angle)}; + EXPECT_NEAR(sgl::angle_between(sgl::quat_identity, q), angle, angle * 1e-4f) << "angle=" << angle; + EXPECT_NEAR(sgl::angle_between(q, sgl::quat_identity), angle, angle * 1e-4f) << "angle=" << angle; + } + /* q and -q are the same rotation */ + const auto q{sgl::quat_from_axis_angle(sgl::vec3{0, 1, 0}, 0.25f)}; + EXPECT_NEAR(sgl::angle_between(sgl::quat_identity, sgl::negate(q)), 0.25f, 1e-5f); +} + TEST(Quat, IsAxisAligned) { EXPECT_TRUE(sgl::is_axis_aligned(sgl::vec3{1, 0, 0})); EXPECT_TRUE(sgl::is_axis_aligned(sgl::vec3{0, 0, -1})); diff --git a/tests/src/vec_test.cpp b/tests/src/vec_test.cpp index 1bca5bc..6463752 100644 --- a/tests/src/vec_test.cpp +++ b/tests/src/vec_test.cpp @@ -4,6 +4,7 @@ */ #include #include +#include import sgl; @@ -27,6 +28,32 @@ TEST(Vec4, NearlyEqual) { EXPECT_TRUE(sgl::nearly_equal(sgl::vec4{1, 2, 3, 4}, sgl::vec4{1, 2, 3, 4})); } +/* Regression: |inf - x| = inf passed the relative check against inf * rel_tol. */ +TEST(Scalar, NearlyEqualInfinity) { + constexpr float inf{std::numeric_limits::infinity()}; + EXPECT_FALSE(sgl::nearly_equal(inf, 0.0f)); + EXPECT_FALSE(sgl::nearly_equal(inf, -inf)); + EXPECT_FALSE(sgl::nearly_equal(inf, 1e30f)); + EXPECT_FALSE(sgl::nearly_equal(-1e30f, -inf)); + EXPECT_TRUE(sgl::nearly_equal(inf, inf)); + EXPECT_TRUE(sgl::nearly_equal(-inf, -inf)); + EXPECT_FALSE(sgl::nearly_equal(sgl::vec3{inf, 0, 0}, sgl::vec3{0, 0, 0})); + EXPECT_TRUE(sgl::nearly_equal(sgl::vec3{inf, 1, 2}, sgl::vec3{inf, 1, 2})); +} + +/* Regression: (&x)[i] read past a scalar member, which also broke constant evaluation. */ +TEST(Vec, IndexingIsConstexpr) { + static_assert(sgl::vec2{1, 2}[1] == 2.0f); + static_assert(sgl::vec3{1, 2, 3}[1] == 2.0f && sgl::vec3{1, 2, 3}[2] == 3.0f); + static_assert(sgl::vec4{1, 2, 3, 4}[3] == 4.0f); + sgl::vec3 v{1, 2, 3}; + v[2] = 5.0f; + EXPECT_EQ(v.z, 5.0f); + sgl::vec4 w{}; + w[3] = 7.0f; + EXPECT_EQ(w.w, 7.0f); +} + TEST(Vec2, NearlyZero) { EXPECT_TRUE(sgl::nearly_zero(sgl::vec2{0, 0})); EXPECT_FALSE(sgl::nearly_zero(sgl::vec2{1, 0})); @@ -326,6 +353,19 @@ TEST(Vec2, IsParallel) { EXPECT_FALSE(sgl::is_parallel(sgl::vec2{1, 0}, sgl::vec2{0, 1})); } +/* Regression: vec2/vec3 compared |cross| to an absolute tolerance, so short perpendicular vectors + * read as parallel and long nearly parallel ones did not; vec4 was relative. */ +TEST(Vec, IsParallelIsScaleInvariant) { + EXPECT_FALSE(sgl::is_parallel(sgl::vec3{1e-3f, 0, 0}, sgl::vec3{0, 1e-3f, 0})); + EXPECT_TRUE(sgl::is_parallel(sgl::vec3{1000, 0, 0}, sgl::vec3{1000, 0.01f, 0})); + EXPECT_FALSE(sgl::is_parallel(sgl::vec2{1e-3f, 0}, sgl::vec2{0, 1e-3f})); + EXPECT_TRUE(sgl::is_parallel(sgl::vec2{1000, 0}, sgl::vec2{1000, 0.01f})); + /* the same answer in every dimension */ + EXPECT_FALSE(sgl::is_parallel(sgl::vec4{1e-3f, 0, 0, 0}, sgl::vec4{0, 1e-3f, 0, 0})); + EXPECT_TRUE(sgl::is_parallel(sgl::vec4{1000, 0, 0, 0}, sgl::vec4{1000, 0.01f, 0, 0})); + EXPECT_TRUE(sgl::is_parallel(sgl::vec3{1, 2, 3}, sgl::vec3{-2, -4, -6})); +} + /* --- Comparison --- */ TEST(Vec3, AllGreaterThan) { @@ -379,3 +419,12 @@ TEST(Vec2, AngleBetween) { auto a = sgl::angle_between(sgl::vec2{1, 0}, sgl::vec2{-1, 0}); EXPECT_NEAR(a, std::acos(-1.0f), 1e-5f); /* pi */ } + +/* Regression: acos near 1 returned 0 for a 1e-4 rad angle. */ +TEST(Vec, AngleBetweenSmallAngles) { + EXPECT_NEAR(sgl::angle_between(sgl::vec3{1, 0, 0}, sgl::vec3{1, 1e-4f, 0}), 1e-4f, 1e-8f); + EXPECT_NEAR(sgl::angle_between(sgl::vec2{1, 0}, sgl::vec2{1, 1e-4f}), 1e-4f, 1e-8f); + EXPECT_NEAR(sgl::angle_between(sgl::vec4{1, 0, 0, 0}, sgl::vec4{1, 0, 0, 1e-4f}), 1e-4f, 1e-8f); + EXPECT_NEAR(sgl::angle_between(sgl::vec3{3, 0, 0}, sgl::vec3{-2, 1e-4f, 0}), std::acos(-1.0f) - 5e-5f, 1e-6f); + EXPECT_EQ(sgl::angle_between(sgl::vec3{0, 0, 0}, sgl::vec3{1, 0, 0}), 0.0f) << "zero vector"; +}