Make bevy_math's libm feature use libm for all f32methods with unspecified precision (#14693)

# Objective

Closes #14474

Previously, the `libm` feature of bevy_math would just pass the same
feature flag down to glam. However, bevy_math itself had many uses of
floating-point arithmetic with unspecified precision. For example,
`f32::sin_cos` and `f32::powi` have unspecified precision, which means
that the exact details of their output are not guaranteed to be stable
across different systems and/or versions of Rust. This means that users
of bevy_math could observe slightly different behavior on different
systems if these methods were used.

The goal of this PR is to make it so that the `libm` feature flag
actually guarantees some degree of determinacy within bevy_math itself
by switching to the libm versions of these functions when the `libm`
feature is enabled.

## Solution

bevy_math now has an internal module `bevy_math::ops`, which re-exports
either the standard versions of the operations or the libm versions
depending on whether the `libm` feature is enabled. For example,
`ops::sin` compiles to `f32::sin` without the `libm` feature and to
`libm::sinf` with it.

This approach has a small shortfall, which is that `f32::powi` (integer
powers of floating point numbers) does not have an equivalent in `libm`.
On the other hand, this method is only used for squaring and cubing
numbers in bevy_math. Accordingly, this deficit is covered by the
introduction of a trait `ops::FloatPow`:
```rust
pub(crate) trait FloatPow {
    fn squared(self) -> Self;
    fn cubed(self) -> Self;
}
```

Next, each current usage of the unspecified-precision methods has been
replaced by its equivalent in `ops`, so that when `libm` is enabled, the
libm version is used instead. The exception, of course, is that
`.powi(2)`/`.powi(3)` have been replaced with `.squared()`/`.cubed()`.

Finally, the usage of the plain `f32` methods with unspecified precision
is now linted out of bevy_math (and hence disallowed in CI). For
example, using `f32::sin` within bevy_math produces a warning that tells
the user to use the `ops::sin` version instead.

## Testing

Ran existing tests. It would be nice to check some benchmarks on NURBS
things once #14677 merges. I'm happy to wait until then if the rest of
this PR is fine.

---

## Discussion

In the future, it might make sense to actually expose `bevy_math::ops`
as public if any downstream Bevy crates want to provide similar
determinacy guarantees. For now, it's all just `pub(crate)`.

This PR also only covers `f32`. If we find ourselves using `f64`
internally in parts of bevy_math for better robustness, we could extend
the module and lints to cover the `f64` versions easily enough.

I don't know how feasible it is, but it would also be nice if we could
standardize the bevy_math tests with the `libm` feature in CI, since
their success is currently platform-dependent (e.g. 8 of them fail on my
machine when run locally).

---------

Co-authored-by: IQuick 143 <IQuick143cz@gmail.com>
This commit is contained in:
Matty 2024-08-12 12:13:36 -04:00 committed by GitHub
parent c8d30edf1a
commit 61a1530c56
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
18 changed files with 450 additions and 119 deletions

View file

@ -0,0 +1,29 @@
disallowed-methods = [
{ path = "f32::powi", reason = "use ops::FloatPow::squared, ops::FloatPow::cubed, or ops::powf instead for libm determinism" },
{ path = "f32::log", reason = "use ops::ln, ops::log2, or ops::log10 instead for libm determinism" },
{ path = "f32::abs_sub", reason = "deprecated and deeply confusing method" },
{ path = "f32::powf", reason = "use ops::powf instead for libm determinism" },
{ path = "f32::exp", reason = "use ops::exp instead for libm determinism" },
{ path = "f32::exp2", reason = "use ops::exp2 instead for libm determinism" },
{ path = "f32::ln", reason = "use ops::ln instead for libm determinism" },
{ path = "f32::log2", reason = "use ops::log2 instead for libm determinism" },
{ path = "f32::log10", reason = "use ops::log10 instead for libm determinism" },
{ path = "f32::cbrt", reason = "use ops::cbrt instead for libm determinism" },
{ path = "f32::hypot", reason = "use ops::hypot instead for libm determinism" },
{ path = "f32::sin", reason = "use ops::sin instead for libm determinism" },
{ path = "f32::cos", reason = "use ops::cos instead for libm determinism" },
{ path = "f32::tan", reason = "use ops::tan instead for libm determinism" },
{ path = "f32::asin", reason = "use ops::asin instead for libm determinism" },
{ path = "f32::acos", reason = "use ops::acos instead for libm determinism" },
{ path = "f32::atan", reason = "use ops::atan instead for libm determinism" },
{ path = "f32::atan2", reason = "use ops::atan2 instead for libm determinism" },
{ path = "f32::sin_cos", reason = "use ops::sin_cos instead for libm determinism" },
{ path = "f32::exp_m1", reason = "use ops::exp_m1 instead for libm determinism" },
{ path = "f32::ln_1p", reason = "use ops::ln_1p instead for libm determinism" },
{ path = "f32::sinh", reason = "use ops::sinh instead for libm determinism" },
{ path = "f32::cosh", reason = "use ops::cosh instead for libm determinism" },
{ path = "f32::tanh", reason = "use ops::tanh instead for libm determinism" },
{ path = "f32::asinh", reason = "use ops::asinh instead for libm determinism" },
{ path = "f32::acosh", reason = "use ops::acosh instead for libm determinism" },
{ path = "f32::atanh", reason = "use ops::atanh instead for libm determinism" },
]

View file

@ -2,6 +2,7 @@ mod primitive_impls;
use super::{BoundingVolume, IntersectsVolume};
use crate::{
ops::FloatPow,
prelude::{Mat2, Rot2, Vec2},
Isometry2d,
};
@ -251,7 +252,7 @@ impl IntersectsVolume<BoundingCircle> for Aabb2d {
fn intersects(&self, circle: &BoundingCircle) -> bool {
let closest_point = self.closest_point(circle.center);
let distance_squared = circle.center.distance_squared(closest_point);
let radius_squared = circle.radius().powi(2);
let radius_squared = circle.radius().squared();
distance_squared <= radius_squared
}
}
@ -261,7 +262,7 @@ mod aabb2d_tests {
use super::Aabb2d;
use crate::{
bounding::{BoundingCircle, BoundingVolume, IntersectsVolume},
Vec2,
ops, Vec2,
};
#[test]
@ -385,7 +386,7 @@ mod aabb2d_tests {
max: Vec2::new(2.0, 2.0),
};
let transformed = a.transformed_by(Vec2::new(2.0, -2.0), std::f32::consts::FRAC_PI_4);
let half_length = 2_f32.hypot(2.0);
let half_length = ops::hypot(2.0, 2.0);
assert_eq!(
transformed.min,
Vec2::new(2.0 - half_length, -half_length - 2.0)
@ -535,7 +536,7 @@ impl BoundingVolume for BoundingCircle {
#[inline(always)]
fn contains(&self, other: &Self) -> bool {
let diff = self.radius() - other.radius();
self.center.distance_squared(other.center) <= diff.powi(2).copysign(diff)
self.center.distance_squared(other.center) <= diff.squared().copysign(diff)
}
#[inline(always)]
@ -593,7 +594,7 @@ impl IntersectsVolume<Self> for BoundingCircle {
#[inline(always)]
fn intersects(&self, other: &Self) -> bool {
let center_distance_squared = self.center.distance_squared(other.center);
let radius_sum_squared = (self.radius() + other.radius()).powi(2);
let radius_sum_squared = (self.radius() + other.radius()).squared();
center_distance_squared <= radius_sum_squared
}
}

View file

@ -1,6 +1,7 @@
//! Contains [`Bounded2d`] implementations for [geometric primitives](crate::primitives).
use crate::{
ops,
primitives::{
Annulus, Arc2d, BoxedPolygon, BoxedPolyline2d, Capsule2d, Circle, CircularSector,
CircularSegment, Ellipse, Line2d, Plane2d, Polygon, Polyline2d, Rectangle, RegularPolygon,
@ -154,7 +155,7 @@ impl Bounded2d for Ellipse {
let (ux, uy) = (hw * alpha_cos, hw * alpha_sin);
let (vx, vy) = (hh * beta_cos, hh * beta_sin);
let half_size = Vec2::new(ux.hypot(vx), uy.hypot(vy));
let half_size = Vec2::new(ops::hypot(ux, vx), ops::hypot(uy, vy));
Aabb2d::new(isometry.translation, half_size)
}
@ -403,6 +404,7 @@ mod tests {
use crate::{
bounding::Bounded2d,
ops::{self, FloatPow},
primitives::{
Annulus, Arc2d, Capsule2d, Circle, CircularSector, CircularSegment, Ellipse, Line2d,
Plane2d, Polygon, Polyline2d, Rectangle, RegularPolygon, Rhombus, Segment2d,
@ -505,7 +507,7 @@ mod tests {
// The exact coordinates here are not obvious, but can be computed by constructing
// an altitude from the midpoint of the chord to the y-axis and using the right triangle
// similarity theorem.
bounding_circle_center: Vec2::new(-apothem / 2.0, apothem.powi(2)),
bounding_circle_center: Vec2::new(-apothem / 2.0, apothem.squared()),
bounding_circle_radius: 0.5,
},
// Test case: handling of axis-aligned extrema
@ -844,7 +846,7 @@ mod tests {
let bounding_circle = segment.bounding_circle(isometry);
assert_eq!(bounding_circle.center, translation);
assert_eq!(bounding_circle.radius(), 1.0_f32.hypot(0.5));
assert_eq!(bounding_circle.radius(), ops::hypot(1.0, 0.5));
}
#[test]
@ -920,7 +922,7 @@ mod tests {
let bounding_circle = rectangle.bounding_circle(Isometry2d::from_translation(translation));
assert_eq!(bounding_circle.center, translation);
assert_eq!(bounding_circle.radius(), 1.0_f32.hypot(0.5));
assert_eq!(bounding_circle.radius(), ops::hypot(1.0, 0.5));
}
#[test]

View file

@ -7,7 +7,7 @@ use crate::primitives::{
BoxedPolygon, BoxedPolyline2d, Capsule2d, Cuboid, Cylinder, Ellipse, Extrusion, Line2d,
Polygon, Polyline2d, Primitive2d, Rectangle, RegularPolygon, Segment2d, Triangle2d,
};
use crate::{Isometry2d, Isometry3d, Quat, Rot2};
use crate::{ops, Isometry2d, Isometry3d, Quat, Rot2};
use crate::{bounding::Bounded2d, primitives::Circle};
@ -231,8 +231,8 @@ pub trait BoundedExtrusion: Primitive2d + Bounded2d {
center,
circle: Circle { radius },
} = self.bounding_circle(Isometry2d::IDENTITY);
let radius = radius.hypot(half_depth);
let center = isometry.translation + isometry.rotation * Vec3A::from(center.extend(0.));
let radius = ops::hypot(radius, half_depth);
let center = isometry * Vec3A::from(center.extend(0.));
BoundingSphere::new(center, radius)
}
@ -246,6 +246,7 @@ mod tests {
use crate::{
bounding::{Bounded3d, BoundingVolume},
ops,
primitives::{
Capsule2d, Circle, Ellipse, Extrusion, Line2d, Polygon, Polyline2d, Rectangle,
RegularPolygon, Segment2d, Triangle2d,
@ -265,7 +266,7 @@ mod tests {
let bounding_sphere = cylinder.bounding_sphere(isometry);
assert_eq!(bounding_sphere.center, translation.into());
assert_eq!(bounding_sphere.radius(), 1f32.hypot(0.5));
assert_eq!(bounding_sphere.radius(), ops::hypot(1.0, 0.5));
}
#[test]

View file

@ -4,7 +4,7 @@ mod primitive_impls;
use glam::Mat3;
use super::{BoundingVolume, IntersectsVolume};
use crate::{Isometry3d, Quat, Vec3A};
use crate::{ops::FloatPow, Isometry3d, Quat, Vec3A};
#[cfg(feature = "bevy_reflect")]
use bevy_reflect::Reflect;
@ -253,7 +253,7 @@ impl IntersectsVolume<BoundingSphere> for Aabb3d {
fn intersects(&self, sphere: &BoundingSphere) -> bool {
let closest_point = self.closest_point(sphere.center);
let distance_squared = sphere.center.distance_squared(closest_point);
let radius_squared = sphere.radius().powi(2);
let radius_squared = sphere.radius().squared();
distance_squared <= radius_squared
}
}
@ -263,7 +263,7 @@ mod aabb3d_tests {
use super::Aabb3d;
use crate::{
bounding::{BoundingSphere, BoundingVolume, IntersectsVolume},
Quat, Vec3, Vec3A,
ops, Quat, Vec3, Vec3A,
};
#[test]
@ -389,7 +389,7 @@ mod aabb3d_tests {
Vec3A::new(2.0, -2.0, 4.0),
Quat::from_rotation_z(std::f32::consts::FRAC_PI_4),
);
let half_length = 2_f32.hypot(2.0);
let half_length = ops::hypot(2.0, 2.0);
assert_eq!(
transformed.min,
Vec3A::new(2.0 - half_length, -half_length - 2.0, 2.0)
@ -518,7 +518,7 @@ impl BoundingSphere {
let radius = self.radius();
let distance_squared = (point - self.center).length_squared();
if distance_squared <= radius.powi(2) {
if distance_squared <= radius.squared() {
// The point is inside the sphere.
point
} else {
@ -553,7 +553,7 @@ impl BoundingVolume for BoundingSphere {
#[inline(always)]
fn contains(&self, other: &Self) -> bool {
let diff = self.radius() - other.radius();
self.center.distance_squared(other.center) <= diff.powi(2).copysign(diff)
self.center.distance_squared(other.center) <= diff.squared().copysign(diff)
}
#[inline(always)]
@ -621,7 +621,7 @@ impl IntersectsVolume<Self> for BoundingSphere {
#[inline(always)]
fn intersects(&self, other: &Self) -> bool {
let center_distance_squared = self.center.distance_squared(other.center);
let radius_sum_squared = (self.radius() + other.radius()).powi(2);
let radius_sum_squared = (self.radius() + other.radius()).squared();
center_distance_squared <= radius_sum_squared
}
}

View file

@ -4,6 +4,7 @@ use glam::Vec3A;
use crate::{
bounding::{Bounded2d, BoundingCircle},
ops,
primitives::{
BoxedPolyline3d, Capsule3d, Cone, ConicalFrustum, Cuboid, Cylinder, InfinitePlane3d,
Line3d, Polyline3d, Segment3d, Sphere, Torus, Triangle2d, Triangle3d,
@ -137,7 +138,7 @@ impl Bounded3d for Cylinder {
}
fn bounding_sphere(&self, isometry: Isometry3d) -> BoundingSphere {
let radius = self.radius.hypot(self.half_height);
let radius = ops::hypot(self.radius, self.half_height);
BoundingSphere::new(isometry.translation, radius)
}
}
@ -345,7 +346,7 @@ impl Bounded3d for Triangle3d {
#[cfg(test)]
mod tests {
use crate::{bounding::BoundingVolume, Isometry3d};
use crate::{bounding::BoundingVolume, ops, Isometry3d};
use glam::{Quat, Vec3, Vec3A};
use crate::{
@ -441,7 +442,7 @@ mod tests {
let bounding_sphere = segment.bounding_sphere(isometry);
assert_eq!(bounding_sphere.center, translation.into());
assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(0.5));
assert_eq!(bounding_sphere.radius(), ops::hypot(1.0, 0.5));
}
#[test]
@ -461,7 +462,10 @@ mod tests {
let bounding_sphere = polyline.bounding_sphere(isometry);
assert_eq!(bounding_sphere.center, translation.into());
assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(1.0).hypot(1.0));
assert_eq!(
bounding_sphere.radius(),
ops::hypot(ops::hypot(1.0, 1.0), 1.0)
);
}
#[test]
@ -479,7 +483,10 @@ mod tests {
let bounding_sphere = cuboid.bounding_sphere(Isometry3d::from_translation(translation));
assert_eq!(bounding_sphere.center, translation.into());
assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(0.5).hypot(0.5));
assert_eq!(
bounding_sphere.radius(),
ops::hypot(ops::hypot(1.0, 0.5), 0.5)
);
}
#[test]
@ -500,7 +507,7 @@ mod tests {
let bounding_sphere = cylinder.bounding_sphere(isometry);
assert_eq!(bounding_sphere.center, translation.into());
assert_eq!(bounding_sphere.radius(), 1.0_f32.hypot(0.5));
assert_eq!(bounding_sphere.radius(), ops::hypot(1.0, 0.5));
}
#[test]

View file

@ -1,5 +1,5 @@
use super::{Aabb2d, BoundingCircle, IntersectsVolume};
use crate::{Dir2, Ray2d, Vec2};
use crate::{ops::FloatPow, Dir2, Ray2d, Vec2};
#[cfg(feature = "bevy_reflect")]
use bevy_reflect::Reflect;
@ -76,8 +76,8 @@ impl RayCast2d {
let offset = self.ray.origin - circle.center;
let projected = offset.dot(*self.ray.direction);
let closest_point = offset - projected * *self.ray.direction;
let distance_squared = circle.radius().powi(2) - closest_point.length_squared();
if distance_squared < 0. || projected.powi(2).copysign(-projected) < -distance_squared {
let distance_squared = circle.radius().squared() - closest_point.length_squared();
if distance_squared < 0. || projected.squared().copysign(-projected) < -distance_squared {
None
} else {
let toi = -projected - distance_squared.sqrt();

View file

@ -1,5 +1,5 @@
use super::{Aabb3d, BoundingSphere, IntersectsVolume};
use crate::{Dir3A, Ray3d, Vec3A};
use crate::{ops::FloatPow, Dir3A, Ray3d, Vec3A};
#[cfg(feature = "bevy_reflect")]
use bevy_reflect::Reflect;
@ -73,8 +73,8 @@ impl RayCast3d {
let offset = self.origin - sphere.center;
let projected = offset.dot(*self.direction);
let closest_point = offset - projected * *self.direction;
let distance_squared = sphere.radius().powi(2) - closest_point.length_squared();
if distance_squared < 0. || projected.powi(2).copysign(-projected) < -distance_squared {
let distance_squared = sphere.radius().squared() - closest_point.length_squared();
if distance_squared < 0. || projected.squared().copysign(-projected) < -distance_squared {
None
} else {
let toi = -projected - distance_squared.sqrt();

View file

@ -1,6 +1,6 @@
//! This module contains abstract mathematical traits shared by types used in `bevy_math`.
use crate::{Dir2, Dir3, Dir3A, Quat, Rot2, Vec2, Vec3, Vec3A, Vec4};
use crate::{ops, Dir2, Dir3, Dir3A, Quat, Rot2, Vec2, Vec3, Vec3A, Vec4};
use std::fmt::Debug;
use std::ops::{Add, Div, Mul, Neg, Sub};
@ -256,7 +256,7 @@ pub trait StableInterpolate: Clone {
/// object_position.smooth_nudge(&target_position, decay_rate, delta_time);
/// ```
fn smooth_nudge(&mut self, target: &Self, decay_rate: f32, delta: f32) {
self.interpolate_stable_assign(target, 1.0 - f32::exp(-decay_rate * delta));
self.interpolate_stable_assign(target, 1.0 - ops::exp(-decay_rate * delta));
}
}

View file

@ -2,7 +2,7 @@
use std::{fmt::Debug, iter::once};
use crate::{Vec2, VectorSpace};
use crate::{ops::FloatPow, Vec2, VectorSpace};
use itertools::Itertools;
use thiserror::Error;
@ -731,12 +731,12 @@ impl<P: VectorSpace> CubicNurbs<P> {
// t[5] := t_i+2
// t[6] := t_i+3
let m00 = (t[4] - t[3]).powi(2) / ((t[4] - t[2]) * (t[4] - t[1]));
let m02 = (t[3] - t[2]).powi(2) / ((t[5] - t[2]) * (t[4] - t[2]));
let m00 = (t[4] - t[3]).squared() / ((t[4] - t[2]) * (t[4] - t[1]));
let m02 = (t[3] - t[2]).squared() / ((t[5] - t[2]) * (t[4] - t[2]));
let m12 = (3.0 * (t[4] - t[3]) * (t[3] - t[2])) / ((t[5] - t[2]) * (t[4] - t[2]));
let m22 = 3.0 * (t[4] - t[3]).powi(2) / ((t[5] - t[2]) * (t[4] - t[2]));
let m33 = (t[4] - t[3]).powi(2) / ((t[6] - t[3]) * (t[5] - t[3]));
let m32 = -m22 / 3.0 - m33 - (t[4] - t[3]).powi(2) / ((t[5] - t[3]) * (t[5] - t[2]));
let m22 = 3.0 * (t[4] - t[3]).squared() / ((t[5] - t[2]) * (t[4] - t[2]));
let m33 = (t[4] - t[3]).squared() / ((t[6] - t[3]) * (t[5] - t[3]));
let m32 = -m22 / 3.0 - m33 - (t[4] - t[3]).squared() / ((t[5] - t[3]) * (t[5] - t[2]));
[
[m00, 1.0 - m00 - m02, m02, 0.0],
[-3.0 * m00, 3.0 * m00 - m12, m12, 0.0],
@ -1254,7 +1254,7 @@ impl<P: VectorSpace> RationalSegment<P> {
// Position = N/D therefore
// Velocity = (N/D)' = N'/D - N * D'/D^2 = (N' * D - N * D')/D^2
numerator_derivative / denominator
- numerator * (denominator_derivative / denominator.powi(2))
- numerator * (denominator_derivative / denominator.squared())
}
/// Instantaneous acceleration of a point at parametric value `t` in `[0, knot_span)`.
@ -1288,10 +1288,10 @@ impl<P: VectorSpace> RationalSegment<P> {
// Velocity = (N/D)' = N'/D - N * D'/D^2 = (N' * D - N * D')/D^2
// Acceleration = (N/D)'' = ((N' * D - N * D')/D^2)' = N''/D + N' * (-2D'/D^2) + N * (-D''/D^2 + 2D'^2/D^3)
numerator_second_derivative / denominator
+ numerator_derivative * (-2.0 * denominator_derivative / denominator.powi(2))
+ numerator_derivative * (-2.0 * denominator_derivative / denominator.squared())
+ numerator
* (-denominator_second_derivative / denominator.powi(2)
+ 2.0 * denominator_derivative.powi(2) / denominator.powi(3))
* (-denominator_second_derivative / denominator.squared()
+ 2.0 * denominator_derivative.squared() / denominator.cubed())
}
/// Calculate polynomial coefficients for the cubic polynomials using a characteristic matrix.
@ -1512,9 +1512,12 @@ impl<P: VectorSpace> From<CubicCurve<P>> for RationalCurve<P> {
mod tests {
use glam::{vec2, Vec2};
use crate::cubic_splines::{
CubicBSpline, CubicBezier, CubicGenerator, CubicNurbs, CubicSegment, RationalCurve,
RationalGenerator,
use crate::{
cubic_splines::{
CubicBSpline, CubicBezier, CubicGenerator, CubicNurbs, CubicSegment, RationalCurve,
RationalGenerator,
},
ops::{self, FloatPow},
};
/// How close two floats can be and still be considered equal
@ -1540,10 +1543,10 @@ mod tests {
/// Manual, hardcoded function for computing the position along a cubic bezier.
fn cubic_manual(t: f32, points: [Vec2; 4]) -> Vec2 {
let p = points;
p[0] * (1.0 - t).powi(3)
+ 3.0 * p[1] * t * (1.0 - t).powi(2)
+ 3.0 * p[2] * t.powi(2) * (1.0 - t)
+ p[3] * t.powi(3)
p[0] * (1.0 - t).cubed()
+ 3.0 * p[1] * t * (1.0 - t).squared()
+ 3.0 * p[2] * t.squared() * (1.0 - t)
+ p[3] * t.cubed()
}
/// Basic cubic Bezier easing test to verify the shape of the curve.
@ -1674,8 +1677,8 @@ mod tests {
// subjecting ones self to a lot of tedious matrix algebra.
let alpha = FRAC_PI_2;
let leg = 2.0 * f32::sin(alpha / 2.0) / (1.0 + 2.0 * f32::cos(alpha / 2.0));
let weight = (1.0 + 2.0 * f32::cos(alpha / 2.0)) / 3.0;
let leg = 2.0 * ops::sin(alpha / 2.0) / (1.0 + 2.0 * ops::cos(alpha / 2.0));
let weight = (1.0 + 2.0 * ops::cos(alpha / 2.0)) / 3.0;
let points = [
vec2(1.0, 0.0),
vec2(1.0, leg),

View file

@ -596,7 +596,7 @@ where
#[cfg(test)]
mod tests {
use super::*;
use crate::Quat;
use crate::{ops, Quat};
use approx::{assert_abs_diff_eq, AbsDiffEq};
use std::f32::consts::TAU;
@ -616,8 +616,8 @@ mod tests {
assert!(curve.sample_unchecked(2.0).abs_diff_eq(&4.0, f32::EPSILON));
assert!(curve.sample_unchecked(-3.0).abs_diff_eq(&9.0, f32::EPSILON));
let curve = function_curve(interval(0.0, f32::INFINITY).unwrap(), f32::log2);
assert_eq!(curve.sample_unchecked(3.5), f32::log2(3.5));
let curve = function_curve(interval(0.0, f32::INFINITY).unwrap(), ops::log2);
assert_eq!(curve.sample_unchecked(3.5), ops::log2(3.5));
assert!(curve.sample_unchecked(-1.0).is_nan());
assert!(curve.sample(-1.0).is_none());
}
@ -642,10 +642,10 @@ mod tests {
#[test]
fn reparametrization() {
let curve = function_curve(interval(1.0, f32::INFINITY).unwrap(), f32::log2);
let curve = function_curve(interval(1.0, f32::INFINITY).unwrap(), ops::log2);
let reparametrized_curve = curve
.by_ref()
.reparametrize(interval(0.0, f32::INFINITY).unwrap(), f32::exp2);
.reparametrize(interval(0.0, f32::INFINITY).unwrap(), ops::exp2);
assert_abs_diff_eq!(reparametrized_curve.sample_unchecked(3.5), 3.5);
assert_abs_diff_eq!(reparametrized_curve.sample_unchecked(100.0), 100.0);
assert_eq!(
@ -664,8 +664,8 @@ mod tests {
#[test]
fn multiple_maps() {
// Make sure these actually happen in the right order.
let curve = function_curve(interval(0.0, 1.0).unwrap(), f32::exp2);
let first_mapped = curve.map(f32::log2);
let curve = function_curve(interval(0.0, 1.0).unwrap(), ops::exp2);
let first_mapped = curve.map(ops::log2);
let second_mapped = first_mapped.map(|x| x * -2.0);
assert_abs_diff_eq!(second_mapped.sample_unchecked(0.0), 0.0);
assert_abs_diff_eq!(second_mapped.sample_unchecked(0.5), -1.0);
@ -675,8 +675,8 @@ mod tests {
#[test]
fn multiple_reparams() {
// Make sure these happen in the right order too.
let curve = function_curve(interval(0.0, 1.0).unwrap(), f32::exp2);
let first_reparam = curve.reparametrize(interval(1.0, 2.0).unwrap(), f32::log2);
let curve = function_curve(interval(0.0, 1.0).unwrap(), ops::exp2);
let first_reparam = curve.reparametrize(interval(1.0, 2.0).unwrap(), ops::log2);
let second_reparam = first_reparam.reparametrize(interval(0.0, 1.0).unwrap(), |t| t + 1.0);
assert_abs_diff_eq!(second_reparam.sample_unchecked(0.0), 1.0);
assert_abs_diff_eq!(second_reparam.sample_unchecked(0.5), 1.5);

View file

@ -860,6 +860,8 @@ impl approx::UlpsEq for Dir3A {
#[cfg(test)]
mod tests {
use crate::ops;
use super::*;
use approx::assert_relative_eq;
@ -916,7 +918,7 @@ mod tests {
#[test]
fn dir2_renorm() {
// Evil denormalized Rot2
let (sin, cos) = 1.0_f32.sin_cos();
let (sin, cos) = ops::sin_cos(1.0_f32);
let rot2 = Rot2::from_sin_cos(sin * (1.0 + 1e-5), cos * (1.0 + 1e-5));
let mut dir_a = Dir2::X;
let mut dir_b = Dir2::X;

View file

@ -21,6 +21,7 @@ pub mod curve;
mod direction;
mod float_ord;
mod isometry;
mod ops;
pub mod primitives;
mod ray;
mod rects;

290
crates/bevy_math/src/ops.rs Normal file
View file

@ -0,0 +1,290 @@
//! This mod re-exports the correct versions of floating-point operations with
//! unspecified precision in the standard library depending on whether the `libm`
//! crate feature is enabled.
//!
//! All the functions here are named according to their versions in the standard
//! library.
#![allow(dead_code)]
#![allow(clippy::disallowed_methods)]
// Note: There are some Rust methods with unspecified precision without a `libm`
// equivalent:
// - `f32::powi` (integer powers)
// - `f32::log` (logarithm with specified base)
// - `f32::abs_sub` (actually unsure if `libm` has this, but don't use it regardless)
//
// Additionally, the following nightly API functions are not presently integrated
// into this, but they would be candidates once standardized:
// - `f32::gamma`
// - `f32::ln_gamma`
#[cfg(not(feature = "libm"))]
mod std_ops {
#[inline(always)]
pub(crate) fn powf(x: f32, y: f32) -> f32 {
f32::powf(x, y)
}
#[inline(always)]
pub(crate) fn exp(x: f32) -> f32 {
f32::exp(x)
}
#[inline(always)]
pub(crate) fn exp2(x: f32) -> f32 {
f32::exp2(x)
}
#[inline(always)]
pub(crate) fn ln(x: f32) -> f32 {
f32::ln(x)
}
#[inline(always)]
pub(crate) fn log2(x: f32) -> f32 {
f32::log2(x)
}
#[inline(always)]
pub(crate) fn log10(x: f32) -> f32 {
f32::log10(x)
}
#[inline(always)]
pub(crate) fn cbrt(x: f32) -> f32 {
f32::cbrt(x)
}
#[inline(always)]
pub(crate) fn hypot(x: f32, y: f32) -> f32 {
f32::hypot(x, y)
}
#[inline(always)]
pub(crate) fn sin(x: f32) -> f32 {
f32::sin(x)
}
#[inline(always)]
pub(crate) fn cos(x: f32) -> f32 {
f32::cos(x)
}
#[inline(always)]
pub(crate) fn tan(x: f32) -> f32 {
f32::tan(x)
}
#[inline(always)]
pub(crate) fn asin(x: f32) -> f32 {
f32::asin(x)
}
#[inline(always)]
pub(crate) fn acos(x: f32) -> f32 {
f32::acos(x)
}
#[inline(always)]
pub(crate) fn atan(x: f32) -> f32 {
f32::atan(x)
}
#[inline(always)]
pub(crate) fn atan2(x: f32, y: f32) -> f32 {
f32::atan2(x, y)
}
#[inline(always)]
pub(crate) fn sin_cos(x: f32) -> (f32, f32) {
f32::sin_cos(x)
}
#[inline(always)]
pub(crate) fn exp_m1(x: f32) -> f32 {
f32::exp_m1(x)
}
#[inline(always)]
pub(crate) fn ln_1p(x: f32) -> f32 {
f32::ln_1p(x)
}
#[inline(always)]
pub(crate) fn sinh(x: f32) -> f32 {
f32::sinh(x)
}
#[inline(always)]
pub(crate) fn cosh(x: f32) -> f32 {
f32::cosh(x)
}
#[inline(always)]
pub(crate) fn tanh(x: f32) -> f32 {
f32::tanh(x)
}
#[inline(always)]
pub(crate) fn asinh(x: f32) -> f32 {
f32::asinh(x)
}
#[inline(always)]
pub(crate) fn acosh(x: f32) -> f32 {
f32::acosh(x)
}
#[inline(always)]
pub(crate) fn atanh(x: f32) -> f32 {
f32::atanh(x)
}
}
#[cfg(feature = "libm")]
mod libm_ops {
#[inline(always)]
pub(crate) fn powf(x: f32, y: f32) -> f32 {
libm::powf(x, y)
}
#[inline(always)]
pub(crate) fn exp(x: f32) -> f32 {
libm::expf(x)
}
#[inline(always)]
pub(crate) fn exp2(x: f32) -> f32 {
libm::exp2f(x)
}
#[inline(always)]
pub(crate) fn ln(x: f32) -> f32 {
// This isn't documented in `libm` but this is actually the base e logarithm.
libm::logf(x)
}
#[inline(always)]
pub(crate) fn log2(x: f32) -> f32 {
libm::log2f(x)
}
#[inline(always)]
pub(crate) fn log10(x: f32) -> f32 {
libm::log10f(x)
}
#[inline(always)]
pub(crate) fn cbrt(x: f32) -> f32 {
libm::cbrtf(x)
}
#[inline(always)]
pub(crate) fn hypot(x: f32, y: f32) -> f32 {
libm::hypotf(x, y)
}
#[inline(always)]
pub(crate) fn sin(x: f32) -> f32 {
libm::sinf(x)
}
#[inline(always)]
pub(crate) fn cos(x: f32) -> f32 {
libm::cosf(x)
}
#[inline(always)]
pub(crate) fn tan(x: f32) -> f32 {
libm::tanf(x)
}
#[inline(always)]
pub(crate) fn asin(x: f32) -> f32 {
libm::asinf(x)
}
#[inline(always)]
pub(crate) fn acos(x: f32) -> f32 {
libm::acosf(x)
}
#[inline(always)]
pub(crate) fn atan(x: f32) -> f32 {
libm::atanf(x)
}
#[inline(always)]
pub(crate) fn atan2(x: f32, y: f32) -> f32 {
libm::atan2f(x, y)
}
#[inline(always)]
pub(crate) fn sin_cos(x: f32) -> (f32, f32) {
libm::sincosf(x)
}
#[inline(always)]
pub(crate) fn exp_m1(x: f32) -> f32 {
libm::expm1f(x)
}
#[inline(always)]
pub(crate) fn ln_1p(x: f32) -> f32 {
libm::log1pf(x)
}
#[inline(always)]
pub(crate) fn sinh(x: f32) -> f32 {
libm::sinhf(x)
}
#[inline(always)]
pub(crate) fn cosh(x: f32) -> f32 {
libm::coshf(x)
}
#[inline(always)]
pub(crate) fn tanh(x: f32) -> f32 {
libm::tanhf(x)
}
#[inline(always)]
pub(crate) fn asinh(x: f32) -> f32 {
libm::asinhf(x)
}
#[inline(always)]
pub(crate) fn acosh(x: f32) -> f32 {
libm::acoshf(x)
}
#[inline(always)]
pub(crate) fn atanh(x: f32) -> f32 {
libm::atanhf(x)
}
}
#[cfg(feature = "libm")]
pub(crate) use libm_ops::*;
#[cfg(not(feature = "libm"))]
pub(crate) use std_ops::*;
/// This extension trait covers shortfall in determinacy from the lack of a `libm` counterpart
/// to `f32::powi`. Use this for the common small exponents.
pub(crate) trait FloatPow {
fn squared(self) -> Self;
fn cubed(self) -> Self;
}
impl FloatPow for f32 {
#[inline]
fn squared(self) -> Self {
self * self
}
#[inline]
fn cubed(self) -> Self {
self * self * self
}
}

View file

@ -1,7 +1,10 @@
use std::f32::consts::{FRAC_PI_2, FRAC_PI_3, PI};
use std::f32::consts::{FRAC_1_SQRT_2, FRAC_PI_2, FRAC_PI_3, PI};
use super::{Measured2d, Primitive2d, WindingOrder};
use crate::{Dir2, Vec2};
use crate::{
ops::{self, FloatPow},
Dir2, Vec2,
};
#[cfg(feature = "bevy_reflect")]
use bevy_reflect::{std_traits::ReflectDefault, Reflect};
@ -54,7 +57,7 @@ impl Circle {
pub fn closest_point(&self, point: Vec2) -> Vec2 {
let distance_squared = point.length_squared();
if distance_squared <= self.radius.powi(2) {
if distance_squared <= self.radius.squared() {
// The point is inside the circle.
point
} else {
@ -70,7 +73,7 @@ impl Measured2d for Circle {
/// Get the area of the circle
#[inline(always)]
fn area(&self) -> f32 {
PI * self.radius.powi(2)
PI * self.radius.squared()
}
/// Get the perimeter or circumference of the circle
@ -200,7 +203,7 @@ impl Arc2d {
/// Get half the distance between the endpoints (half the length of the chord)
#[inline(always)]
pub fn half_chord_length(&self) -> f32 {
self.radius * f32::sin(self.half_angle)
self.radius * ops::sin(self.half_angle)
}
/// Get the distance between the endpoints (the length of the chord)
@ -226,7 +229,7 @@ impl Arc2d {
// used by Wolfram MathWorld, which is the distance rather than the segment.
pub fn apothem(&self) -> f32 {
let sign = if self.is_minor() { 1.0 } else { -1.0 };
sign * f32::sqrt(self.radius.powi(2) - self.half_chord_length().powi(2))
sign * f32::sqrt(self.radius.squared() - self.half_chord_length().squared())
}
/// Get the length of the sagitta of this arc, that is,
@ -388,7 +391,7 @@ impl CircularSector {
/// Returns the area of this sector
#[inline(always)]
pub fn area(&self) -> f32 {
self.arc.radius.powi(2) * self.arc.half_angle
self.arc.radius.squared() * self.arc.half_angle
}
}
@ -527,7 +530,7 @@ impl CircularSegment {
/// Returns the area of this segment
#[inline(always)]
pub fn area(&self) -> f32 {
0.5 * self.arc.radius.powi(2) * (self.arc.angle() - self.arc.angle().sin())
0.5 * self.arc.radius.squared() * (self.arc.angle() - ops::sin(self.arc.angle()))
}
}
@ -882,11 +885,11 @@ impl Measured2d for Ellipse {
// The algorithm used here is the Gauss-Kummer infinite series expansion of the elliptic integral expression for the perimeter of ellipses
// For more information see https://www.wolframalpha.com/input/?i=gauss-kummer+series
// We only use the terms up to `i == 20` for this approximation
let h = ((a - b) / (a + b)).powi(2);
let h = ((a - b) / (a + b)).squared();
PI * (a + b)
* (0..=20)
.map(|i| BINOMIAL_COEFFICIENTS[i] * h.powi(i as i32))
.map(|i| BINOMIAL_COEFFICIENTS[i] * ops::powf(h, i as f32))
.sum::<f32>()
}
}
@ -953,8 +956,8 @@ impl Annulus {
pub fn closest_point(&self, point: Vec2) -> Vec2 {
let distance_squared = point.length_squared();
if self.inner_circle.radius.powi(2) <= distance_squared {
if distance_squared <= self.outer_circle.radius.powi(2) {
if self.inner_circle.radius.squared() <= distance_squared {
if distance_squared <= self.outer_circle.radius.squared() {
// The point is inside the annulus.
point
} else {
@ -976,7 +979,7 @@ impl Measured2d for Annulus {
/// Get the area of the annulus
#[inline(always)]
fn area(&self) -> f32 {
PI * (self.outer_circle.radius.powi(2) - self.inner_circle.radius.powi(2))
PI * (self.outer_circle.radius.squared() - self.inner_circle.radius.squared())
}
/// Get the perimeter or circumference of the annulus,
@ -1031,7 +1034,7 @@ impl Rhombus {
#[inline(always)]
pub fn from_side(side: f32) -> Self {
Self {
half_diagonals: Vec2::splat(side.hypot(side) / 2.0),
half_diagonals: Vec2::splat(side * FRAC_1_SQRT_2),
}
}
@ -1694,13 +1697,13 @@ impl RegularPolygon {
#[inline(always)]
#[doc(alias = "apothem")]
pub fn inradius(&self) -> f32 {
self.circumradius() * (PI / self.sides as f32).cos()
self.circumradius() * ops::cos(PI / self.sides as f32)
}
/// Get the length of one side of the regular polygon
#[inline(always)]
pub fn side_length(&self) -> f32 {
2.0 * self.circumradius() * (PI / self.sides as f32).sin()
2.0 * self.circumradius() * ops::sin(PI / self.sides as f32)
}
/// Get the internal angle of the regular polygon in degrees.
@ -1750,7 +1753,7 @@ impl RegularPolygon {
(0..self.sides).map(move |i| {
let theta = start_angle + i as f32 * step;
let (sin, cos) = theta.sin_cos();
let (sin, cos) = ops::sin_cos(theta);
Vec2::new(cos, sin) * self.circumcircle.radius
})
}
@ -1761,7 +1764,7 @@ impl Measured2d for RegularPolygon {
#[inline(always)]
fn area(&self) -> f32 {
let angle: f32 = 2.0 * PI / (self.sides as f32);
(self.sides as f32) * self.circumradius().powi(2) * angle.sin() / 2.0
(self.sides as f32) * self.circumradius().squared() * ops::sin(angle) / 2.0
}
/// Get the perimeter of the regular polygon.
@ -1821,7 +1824,7 @@ mod tests {
// Reference values were computed by hand and/or with external tools
use super::*;
use approx::assert_relative_eq;
use approx::{assert_abs_diff_eq, assert_relative_eq};
#[test]
fn rectangle_closest_point() {
@ -1913,10 +1916,10 @@ mod tests {
assert_eq!(rhombus.inradius(), 0.0, "incorrect inradius");
assert_eq!(rhombus.circumradius(), 0.0, "incorrect circumradius");
let rhombus = Rhombus::from_side(std::f32::consts::SQRT_2);
assert_eq!(rhombus, Rhombus::new(2.0, 2.0));
assert_eq!(
rhombus,
Rhombus::from_inradius(std::f32::consts::FRAC_1_SQRT_2)
assert_abs_diff_eq!(rhombus.half_diagonals, Vec2::new(1.0, 1.0));
assert_abs_diff_eq!(
rhombus.half_diagonals,
Rhombus::from_inradius(std::f32::consts::FRAC_1_SQRT_2).half_diagonals
);
}

View file

@ -1,7 +1,7 @@
use std::f32::consts::{FRAC_PI_3, PI};
use super::{Circle, Measured2d, Measured3d, Primitive2d, Primitive3d};
use crate::{Dir3, InvalidDirectionError, Mat3, Vec2, Vec3};
use crate::{ops, ops::FloatPow, Dir3, InvalidDirectionError, Mat3, Vec2, Vec3};
#[cfg(feature = "bevy_reflect")]
use bevy_reflect::{std_traits::ReflectDefault, Reflect};
@ -54,7 +54,7 @@ impl Sphere {
pub fn closest_point(&self, point: Vec3) -> Vec3 {
let distance_squared = point.length_squared();
if distance_squared <= self.radius.powi(2) {
if distance_squared <= self.radius.squared() {
// The point is inside the sphere.
point
} else {
@ -70,13 +70,13 @@ impl Measured3d for Sphere {
/// Get the surface area of the sphere
#[inline(always)]
fn area(&self) -> f32 {
4.0 * PI * self.radius.powi(2)
4.0 * PI * self.radius.squared()
}
/// Get the volume of the sphere
#[inline(always)]
fn volume(&self) -> f32 {
4.0 * FRAC_PI_3 * self.radius.powi(3)
4.0 * FRAC_PI_3 * self.radius.cubed()
}
}
@ -505,7 +505,7 @@ impl Cylinder {
/// Get the surface area of one base of the cylinder
#[inline(always)]
pub fn base_area(&self) -> f32 {
PI * self.radius.powi(2)
PI * self.radius.squared()
}
}
@ -642,7 +642,7 @@ impl Cone {
#[inline(always)]
#[doc(alias = "side_length")]
pub fn slant_height(&self) -> f32 {
self.radius.hypot(self.height)
ops::hypot(self.radius, self.height)
}
/// Get the surface area of the side of the cone,
@ -656,7 +656,7 @@ impl Cone {
/// Get the surface area of the base of the cone
#[inline(always)]
pub fn base_area(&self) -> f32 {
PI * self.radius.powi(2)
PI * self.radius.squared()
}
}
@ -828,14 +828,14 @@ impl Measured3d for Torus {
/// the expected result when the torus has a ring and isn't self-intersecting
#[inline(always)]
fn area(&self) -> f32 {
4.0 * PI.powi(2) * self.major_radius * self.minor_radius
4.0 * PI.squared() * self.major_radius * self.minor_radius
}
/// Get the volume of the torus. Note that this only produces
/// the expected result when the torus has a ring and isn't self-intersecting
#[inline(always)]
fn volume(&self) -> f32 {
2.0 * PI.powi(2) * self.major_radius * self.minor_radius.powi(2)
2.0 * PI.squared() * self.major_radius * self.minor_radius.squared()
}
}

View file

@ -1,6 +1,9 @@
use glam::FloatExt;
use crate::prelude::{Mat2, Vec2};
use crate::{
ops,
prelude::{Mat2, Vec2},
};
#[cfg(feature = "bevy_reflect")]
use bevy_reflect::{std_traits::ReflectDefault, Reflect};
@ -99,14 +102,7 @@ impl Rot2 {
/// Creates a [`Rot2`] from a counterclockwise angle in radians.
#[inline]
pub fn radians(radians: f32) -> Self {
#[cfg(feature = "libm")]
let (sin, cos) = (
libm::sin(radians as f64) as f32,
libm::cos(radians as f64) as f32,
);
#[cfg(not(feature = "libm"))]
let (sin, cos) = radians.sin_cos();
let (sin, cos) = ops::sin_cos(radians);
Self::from_sin_cos(sin, cos)
}
@ -136,14 +132,7 @@ impl Rot2 {
/// Returns the rotation in radians in the `(-pi, pi]` range.
#[inline]
pub fn as_radians(self) -> f32 {
#[cfg(feature = "libm")]
{
libm::atan2(self.sin as f64, self.cos as f64) as f32
}
#[cfg(not(feature = "libm"))]
{
f32::atan2(self.sin, self.cos)
}
ops::atan2(self.sin, self.cos)
}
/// Returns the rotation in degrees in the `(-180, 180]` range.

View file

@ -40,7 +40,7 @@
use std::f32::consts::{PI, TAU};
use crate::{primitives::*, NormedVectorSpace, Vec2, Vec3};
use crate::{ops, primitives::*, NormedVectorSpace, Vec2, Vec3};
use rand::{
distributions::{Distribution, WeightedIndex},
Rng,
@ -155,12 +155,14 @@ impl ShapeSample for Circle {
let theta = rng.gen_range(0.0..TAU);
let r_squared = rng.gen_range(0.0..=(self.radius * self.radius));
let r = r_squared.sqrt();
Vec2::new(r * theta.cos(), r * theta.sin())
let (sin, cos) = ops::sin_cos(theta);
Vec2::new(r * cos, r * sin)
}
fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
let theta = rng.gen_range(0.0..TAU);
Vec2::new(self.radius * theta.cos(), self.radius * theta.sin())
let (sin, cos) = ops::sin_cos(theta);
Vec2::new(self.radius * cos, self.radius * sin)
}
}
@ -168,7 +170,7 @@ impl ShapeSample for Circle {
#[inline]
fn sample_unit_sphere_boundary<R: Rng + ?Sized>(rng: &mut R) -> Vec3 {
let z = rng.gen_range(-1f32..=1f32);
let (a_sin, a_cos) = rng.gen_range(-PI..=PI).sin_cos();
let (a_sin, a_cos) = ops::sin_cos(rng.gen_range(-PI..=PI));
let c = (1f32 - z * z).sqrt();
let x = a_sin * c;
let y = a_cos * c;
@ -181,7 +183,7 @@ impl ShapeSample for Sphere {
fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
let r_cubed = rng.gen_range(0.0..=(self.radius * self.radius * self.radius));
let r = r_cubed.cbrt();
let r = ops::cbrt(r_cubed);
r * sample_unit_sphere_boundary(rng)
}
@ -202,8 +204,9 @@ impl ShapeSample for Annulus {
let r_squared = rng.gen_range((inner_radius * inner_radius)..(outer_radius * outer_radius));
let r = r_squared.sqrt();
let theta = rng.gen_range(0.0..TAU);
let (sin, cos) = ops::sin_cos(theta);
Vec2::new(r * theta.cos(), r * theta.sin())
Vec2::new(r * cos, r * sin)
}
fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
@ -624,7 +627,7 @@ mod tests {
for _ in 0..5000 {
let point = circle.sample_boundary(&mut rng);
let angle = f32::atan(point.y / point.x) + PI / 2.0;
let angle = ops::atan(point.y / point.x) + PI / 2.0;
let wedge = (angle * 8.0 / PI).floor() as usize;
wedge_hits[wedge] += 1;
}