Improve Bevy's double-precision story for third-party crates (#19194)

# Objective

Certain classes of games, usually those with enormous worlds, require
some amount of support for double-precision. Libraries like `big_space`
exist to allow for large worlds while integrating cleanly with Bevy's
primarily single-precision ecosystem, but even then, games will often
still work directly in double-precision throughout the part of the
pipeline that feeds into the Bevy interface.

Currently, working with double-precision types in Bevy is a pain. `glam`
provides types like `DVec3`, but Bevy doesn't provide double-precision
analogs for `glam` wrappers like `Dir3`. This is mostly because doing so
involves one of:

- code duplication
- generics
- templates (like `glam` uses)
- macros

Each of these has issues that are enough to be deal-breakers as far as
maintainability, usability or readability. To work around this, I'm
putting together `bevy_dmath`, a crate that duplicates `bevy_math` types
and functionality to allow downstream users to enjoy the ergonomics and
power of `bevy_math` in double-precision. For the most part, it's a
smooth process, but in order to fully integrate, there are some
necessary changes that can only be made in `bevy_math`.

## Solution

This PR addresses the first and easiest issue with downstream
double-precision math support: `VectorSpace` currently can only
represent vector spaces over `f32`. This automatically closes the door
to double-precision curves, among other things. This restriction can be
easily lifted by allowing vector spaces to specify the underlying scalar
field. This PR adds a new trait `ScalarField` that satisfies the
properties of a scalar field (the ones that can be upheld statically)
and adds a new associated type `type Scalar: ScalarField` to
`VectorSpace`. It's mostly an unintrusive change. The biggest annoyances
are:

- it touches a lot of curve code
- `bevy_math::ops` doesn't support `f64`, so there are some annoying
workarounds

As far as curves code, I wanted to make this change unintrusive and
bite-sized, so I'm trying to touch as little code as possible. To prove
to myself it can be done, I went ahead and (*not* in this PR) migrated
most of the curves API to support different `ScalarField`s and it went
really smoothly! The ugliest thing was adding `P::Scalar: From<usize>`
in several places. There's an argument to be made here that we should be
using `num-traits`, but that's not immediately relevant. The point is
that for now, the smallest change I could make was to go into every
curve impl and make them generic over `VectorSpace<Scalar = f32>`.
Curves work exactly like before and don't change the user API at all.

# Follow-up

- **Extend `bevy_math::ops` to work with `f64`.** `bevy_math::ops` is
used all over, and if curves are ever going to support different
`ScalarField` types, we'll need to be able to use the correct `std` or
`libm` ops for `f64` types as well. Adding an `ops64` mod turned out to
be really ugly, but I'll point out the maintenance burden is low because
we're not going to be adding new floating-point ops anytime soon.
Another solution is to build a floating-point trait that calls the right
op variant and impl it for `f32` and `f64`. This reduces maintenance
burden because on the off chance we ever *do* want to go modify it, it's
all tied together: you can't change the interface on one without
changing the trait, which forces you to update the other. A third option
is to use `num-traits`, which is basically option 2 but someone else did
the work for us. They already support `no_std` using `libm`, so it would
be more or less a drop-in replacement. They're missing a couple
floating-point ops like `floor` and `ceil`, but we could make our own
floating-point traits for those (there's even the potential for
upstreaming them into `num-traits`).
- **Tweak curves to accept vector spaces over any `ScalarField`.**
Curves are ready to support custom scalar types as soon as the bullet
above is addressed. I will admit that the code is not as fun to look at:
`P::Scalar` instead of `f32` everywhere. We could consider an alternate
design where we use `f32` even to interpolate something like a `DVec3`,
but personally I think that's a worse solution than parameterizing
curves over the vector space's scalar type. At the end of the day, it's
not really bad to deal with in my opinion... `ScalarType` supports
enough operations that working with them is almost like working with raw
float types, and it unlocks a whole ecosystem for games that want to use
double-precision.

benches/benches/bevy_math/bezier.rs

@@ -32,7 +32,7 @@ fn segment_ease(c: &mut Criterion) {
 
 fn curve_position(c: &mut Criterion) {
     /// A helper function that benchmarks calling [`CubicCurve::position()`] over a generic [`VectorSpace`].
-    fn bench_curve<M: Measurement, P: VectorSpace>(
+    fn bench_curve<M: Measurement, P: VectorSpace<Scalar = f32>>(
         group: &mut BenchmarkGroup<M>,
         name: &str,
         curve: CubicCurve<P>,

crates/bevy_animation/src/gltf_curves.rs

@@ -55,7 +55,7 @@ pub struct CubicKeyframeCurve<T> {
 
 impl<V> Curve<V> for CubicKeyframeCurve<V>
 where
-    V: VectorSpace,
+    V: VectorSpace<Scalar = f32>,
 {
     #[inline]
     fn domain(&self) -> Interval {
@@ -179,7 +179,7 @@ pub struct WideLinearKeyframeCurve<T> {
 
 impl<T> IterableCurve<T> for WideLinearKeyframeCurve<T>
 where
-    T: VectorSpace,
+    T: VectorSpace<Scalar = f32>,
 {
     #[inline]
     fn domain(&self) -> Interval {
@@ -289,7 +289,7 @@ pub struct WideCubicKeyframeCurve<T> {
 
 impl<T> IterableCurve<T> for WideCubicKeyframeCurve<T>
 where
-    T: VectorSpace,
+    T: VectorSpace<Scalar = f32>,
 {
     #[inline]
     fn domain(&self) -> Interval {
@@ -406,7 +406,7 @@ fn cubic_spline_interpolation<T>(
     step_duration: f32,
 ) -> T
 where
-    T: VectorSpace,
+    T: VectorSpace<Scalar = f32>,
 {
     let coeffs = (vec4(2.0, 1.0, -2.0, 1.0) * lerp + vec4(-3.0, -2.0, 3.0, -1.0)) * lerp;
     value_start * (coeffs.x * lerp + 1.0)
@@ -415,7 +415,7 @@ where
         + tangent_in_end * step_duration * lerp * coeffs.w
 }
 
-fn cubic_spline_interpolate_slices<'a, T: VectorSpace>(
+fn cubic_spline_interpolate_slices<'a, T: VectorSpace<Scalar = f32>>(
     width: usize,
     first: &'a [T],
     second: &'a [T],

crates/bevy_color/src/lib.rs

@@ -262,6 +262,7 @@ macro_rules! impl_componentwise_vector_space {
         }
 
         impl bevy_math::VectorSpace for $ty {
+            type Scalar = f32;
             const ZERO: Self = Self {
                 $($element: 0.0,)+
             };

crates/bevy_math/src/common_traits.rs

@@ -1,6 +1,6 @@
 //! This module contains abstract mathematical traits shared by types used in `bevy_math`.
 
-use crate::{ops, Dir2, Dir3, Dir3A, Quat, Rot2, Vec2, Vec3, Vec3A, Vec4};
+use crate::{ops, DVec2, DVec3, DVec4, Dir2, Dir3, Dir3A, Quat, Rot2, Vec2, Vec3, Vec3A, Vec4};
 use core::{
     fmt::Debug,
     ops::{Add, Div, Mul, Neg, Sub},
@@ -9,7 +9,7 @@ use variadics_please::all_tuples_enumerated;
 
 /// A type that supports the mathematical operations of a real vector space, irrespective of dimension.
 /// In particular, this means that the implementing type supports:
-/// - Scalar multiplication and division on the right by elements of `f32`
+/// - Scalar multiplication and division on the right by elements of `Self::Scalar`
 /// - Negation
 /// - Addition and subtraction
 /// - Zero
@@ -19,16 +19,16 @@ use variadics_please::all_tuples_enumerated;
 /// - (Commutativity of addition) For all `u, v: Self`, `u + v == v + u`.
 /// - (Additive identity) For all `v: Self`, `v + Self::ZERO == v`.
 /// - (Additive inverse) For all `v: Self`, `v - v == v + (-v) == Self::ZERO`.
-/// - (Compatibility of multiplication) For all `a, b: f32`, `v: Self`, `v * (a * b) == (v * a) * b`.
+/// - (Compatibility of multiplication) For all `a, b: Self::Scalar`, `v: Self`, `v * (a * b) == (v * a) * b`.
 /// - (Multiplicative identity) For all `v: Self`, `v * 1.0 == v`.
-/// - (Distributivity for vector addition) For all `a: f32`, `u, v: Self`, `(u + v) * a == u * a + v * a`.
+/// - (Distributivity for vector addition) For all `a: Self::Scalar`, `u, v: Self`, `(u + v) * a == u * a + v * a`.
-/// - (Distributivity for scalar addition) For all `a, b: f32`, `v: Self`, `v * (a + b) == v * a + v * b`.
+/// - (Distributivity for scalar addition) For all `a, b: Self::Scalar`, `v: Self`, `v * (a + b) == v * a + v * b`.
 ///
 /// Note that, because implementing types use floating point arithmetic, they are not required to actually
 /// implement `PartialEq` or `Eq`.
 pub trait VectorSpace:
-    Mul<f32, Output = Self>
+    Mul<Self::Scalar, Output = Self>
-    + Div<f32, Output = Self>
+    + Div<Self::Scalar, Output = Self>
     + Add<Self, Output = Self>
     + Sub<Self, Output = Self>
     + Neg<Output = Self>
@@ -37,6 +37,9 @@ pub trait VectorSpace:
     + Clone
     + Copy
 {
+    /// The scalar type of this vector space.
+    type Scalar: ScalarField;
+
     /// The zero vector, which is the identity of addition for the vector space type.
     const ZERO: Self;
 
@@ -47,29 +50,99 @@ pub trait VectorSpace:
     /// Note that the value of `t` is not clamped by this function, so extrapolating outside
     /// of the interval `[0,1]` is allowed.
     #[inline]
-    fn lerp(self, rhs: Self, t: f32) -> Self {
+    fn lerp(self, rhs: Self, t: Self::Scalar) -> Self {
-        self * (1. - t) + rhs * t
+        self * (Self::Scalar::ONE - t) + rhs * t
     }
 }
 
 impl VectorSpace for Vec4 {
+    type Scalar = f32;
     const ZERO: Self = Vec4::ZERO;
 }
 
 impl VectorSpace for Vec3 {
+    type Scalar = f32;
     const ZERO: Self = Vec3::ZERO;
 }
 
 impl VectorSpace for Vec3A {
+    type Scalar = f32;
     const ZERO: Self = Vec3A::ZERO;
 }
 
 impl VectorSpace for Vec2 {
+    type Scalar = f32;
     const ZERO: Self = Vec2::ZERO;
 }
 
-impl VectorSpace for f32 {
+impl VectorSpace for DVec4 {
+    type Scalar = f64;
+    const ZERO: Self = DVec4::ZERO;
+}
+
+impl VectorSpace for DVec3 {
+    type Scalar = f64;
+    const ZERO: Self = DVec3::ZERO;
+}
+
+impl VectorSpace for DVec2 {
+    type Scalar = f64;
+    const ZERO: Self = DVec2::ZERO;
+}
+
+// Every scalar field is a 1-dimensional vector space over itself.
+impl<T: ScalarField> VectorSpace for T {
+    type Scalar = Self;
+    const ZERO: Self = Self::ZERO;
+}
+
+/// A type that supports the operations of a scalar field. An implementation should support:
+/// - Addition and subtraction
+/// - Multiplication and division
+/// - Negation
+/// - Zero (additive identity)
+/// - One (multiplicative identity)
+///
+/// Within the limitations of floating point arithmetic, all the following are required to hold:
+/// - (Associativity of addition) For all `u, v, w: Self`, `(u + v) + w == u + (v + w)`.
+/// - (Commutativity of addition) For all `u, v: Self`, `u + v == v + u`.
+/// - (Additive identity) For all `v: Self`, `v + Self::ZERO == v`.
+/// - (Additive inverse) For all `v: Self`, `v - v == v + (-v) == Self::ZERO`.
+/// - (Associativity of multiplication) For all `u, v, w: Self`, `(u * v) * w == u * (v * w)`.
+/// - (Commutativity of multiplication) For all `u, v: Self`, `u * v == v * u`.
+/// - (Multiplicative identity) For all `v: Self`, `v * Self::ONE == v`.
+/// - (Multiplicative inverse) For all `v: Self`, `v / v == v * v.inverse() == Self::ONE`.
+/// - (Distributivity over addition) For all `a, b: Self`, `u, v: Self`, `(u + v) * a == u * a + v * a`.
+pub trait ScalarField:
+    Mul<Self, Output = Self>
+    + Div<Self, Output = Self>
+    + Add<Self, Output = Self>
+    + Sub<Self, Output = Self>
+    + Neg<Output = Self>
+    + Default
+    + Debug
+    + Clone
+    + Copy
+{
+    /// The additive identity.
+    const ZERO: Self;
+    /// The multiplicative identity.
+    const ONE: Self;
+
+    /// The multiplicative inverse of this element. This is equivalent to `1.0 / self`.
+    fn recip(self) -> Self {
+        Self::ONE / self
+    }
+}
+
+impl ScalarField for f32 {
     const ZERO: Self = 0.0;
+    const ONE: Self = 1.0;
+}
+
+impl ScalarField for f64 {
+    const ZERO: Self = 0.0;
+    const ONE: Self = 1.0;
 }
 
 /// A type consisting of formal sums of elements from `V` and `W`. That is,
@@ -84,24 +157,24 @@ impl VectorSpace for f32 {
 #[cfg_attr(feature = "bevy_reflect", derive(bevy_reflect::Reflect))]
 pub struct Sum<V, W>(pub V, pub W);
 
-impl<V, W> Mul<f32> for Sum<V, W>
+impl<F: ScalarField, V, W> Mul<F> for Sum<V, W>
 where
-    V: VectorSpace,
+    V: VectorSpace<Scalar = F>,
-    W: VectorSpace,
+    W: VectorSpace<Scalar = F>,
 {
     type Output = Self;
-    fn mul(self, rhs: f32) -> Self::Output {
+    fn mul(self, rhs: F) -> Self::Output {
         Sum(self.0 * rhs, self.1 * rhs)
     }
 }
 
-impl<V, W> Div<f32> for Sum<V, W>
+impl<F: ScalarField, V, W> Div<F> for Sum<V, W>
 where
-    V: VectorSpace,
+    V: VectorSpace<Scalar = F>,
-    W: VectorSpace,
+    W: VectorSpace<Scalar = F>,
 {
     type Output = Self;
-    fn div(self, rhs: f32) -> Self::Output {
+    fn div(self, rhs: F) -> Self::Output {
         Sum(self.0 / rhs, self.1 / rhs)
     }
 }
@@ -149,11 +222,12 @@ where
     }
 }
 
-impl<V, W> VectorSpace for Sum<V, W>
+impl<F: ScalarField, V, W> VectorSpace for Sum<V, W>
 where
-    V: VectorSpace,
+    V: VectorSpace<Scalar = F>,
-    W: VectorSpace,
+    W: VectorSpace<Scalar = F>,
 {
+    type Scalar = F;
     const ZERO: Self = Sum(V::ZERO, W::ZERO);
 }
 
@@ -162,32 +236,32 @@ where
 /// relationships hold, within the limitations of floating point arithmetic:
 /// - (Nonnegativity) For all `v: Self`, `v.norm() >= 0.0`.
 /// - (Positive definiteness) For all `v: Self`, `v.norm() == 0.0` implies `v == Self::ZERO`.
-/// - (Absolute homogeneity) For all `c: f32`, `v: Self`, `(v * c).norm() == v.norm() * c.abs()`.
+/// - (Absolute homogeneity) For all `c: Self::Scalar`, `v: Self`, `(v * c).norm() == v.norm() * c.abs()`.
 /// - (Triangle inequality) For all `v, w: Self`, `(v + w).norm() <= v.norm() + w.norm()`.
 ///
 /// Note that, because implementing types use floating point arithmetic, they are not required to actually
 /// implement `PartialEq` or `Eq`.
 pub trait NormedVectorSpace: VectorSpace {
     /// The size of this element. The return value should always be nonnegative.
-    fn norm(self) -> f32;
+    fn norm(self) -> Self::Scalar;
 
     /// The squared norm of this element. Computing this is often faster than computing
     /// [`NormedVectorSpace::norm`].
     #[inline]
-    fn norm_squared(self) -> f32 {
+    fn norm_squared(self) -> Self::Scalar {
         self.norm() * self.norm()
     }
 
     /// The distance between this element and another, as determined by the norm.
     #[inline]
-    fn distance(self, rhs: Self) -> f32 {
+    fn distance(self, rhs: Self) -> Self::Scalar {
         (rhs - self).norm()
     }
 
     /// The squared distance between this element and another, as determined by the norm. Note that
     /// this is often faster to compute in practice than [`NormedVectorSpace::distance`].
     #[inline]
-    fn distance_squared(self, rhs: Self) -> f32 {
+    fn distance_squared(self, rhs: Self) -> Self::Scalar {
         (rhs - self).norm_squared()
     }
 }
@@ -245,10 +319,55 @@ impl NormedVectorSpace for f32 {
     fn norm(self) -> f32 {
         ops::abs(self)
     }
+}
+
+impl NormedVectorSpace for DVec4 {
+    #[inline]
+    fn norm(self) -> f64 {
+        self.length()
+    }
 
     #[inline]
-    fn norm_squared(self) -> f32 {
+    fn norm_squared(self) -> f64 {
-        self * self
+        self.length_squared()
+    }
+}
+
+impl NormedVectorSpace for DVec3 {
+    #[inline]
+    fn norm(self) -> f64 {
+        self.length()
+    }
+
+    #[inline]
+    fn norm_squared(self) -> f64 {
+        self.length_squared()
+    }
+}
+
+impl NormedVectorSpace for DVec2 {
+    #[inline]
+    fn norm(self) -> f64 {
+        self.length()
+    }
+
+    #[inline]
+    fn norm_squared(self) -> f64 {
+        self.length_squared()
+    }
+}
+
+impl NormedVectorSpace for f64 {
+    #[inline]
+    #[cfg(feature = "std")]
+    fn norm(self) -> f64 {
+        f64::abs(self)
+    }
+
+    #[inline]
+    #[cfg(all(any(feature = "libm", feature = "nostd-libm"), not(feature = "std")))]
+    fn norm(self) -> f64 {
+        libm::fabs(self)
     }
 }
 
@@ -353,7 +472,7 @@ pub trait StableInterpolate: Clone {
 // VectorSpace type, but the "natural from the semantics" part is less clear in general.
 impl<V> StableInterpolate for V
 where
-    V: NormedVectorSpace,
+    V: NormedVectorSpace<Scalar = f32>,
 {
     #[inline]
     fn interpolate_stable(&self, other: &Self, t: f32) -> Self {
@@ -462,10 +581,13 @@ impl<V: VectorSpace> HasTangent for V {
     type Tangent = V;
 }
 
-impl<M, N> HasTangent for (M, N)
+impl<F, U, V, M, N> HasTangent for (M, N)
 where
-    M: HasTangent,
+    F: ScalarField,
-    N: HasTangent,
+    U: VectorSpace<Scalar = F>,
+    V: VectorSpace<Scalar = F>,
+    M: HasTangent<Tangent = U>,
+    N: HasTangent<Tangent = V>,
 {
     type Tangent = Sum<M::Tangent, N::Tangent>;
 }

crates/bevy_math/src/cubic_splines/curve_impls.rs

@@ -10,7 +10,7 @@ use super::{CubicCurve, RationalCurve};
 
 // -- CubicSegment
 
-impl<P: VectorSpace> Curve<P> for CubicSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> Curve<P> for CubicSegment<P> {
     #[inline]
     fn domain(&self) -> Interval {
         Interval::UNIT
@@ -22,7 +22,7 @@ impl<P: VectorSpace> Curve<P> for CubicSegment<P> {
     }
 }
 
-impl<P: VectorSpace> SampleDerivative<P> for CubicSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleDerivative<P> for CubicSegment<P> {
     #[inline]
     fn sample_with_derivative_unchecked(&self, t: f32) -> WithDerivative<P> {
         WithDerivative {
@@ -32,7 +32,7 @@ impl<P: VectorSpace> SampleDerivative<P> for CubicSegment<P> {
     }
 }
 
-impl<P: VectorSpace> SampleTwoDerivatives<P> for CubicSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleTwoDerivatives<P> for CubicSegment<P> {
     #[inline]
     fn sample_with_two_derivatives_unchecked(&self, t: f32) -> WithTwoDerivatives<P> {
         WithTwoDerivatives {
@@ -46,7 +46,7 @@ impl<P: VectorSpace> SampleTwoDerivatives<P> for CubicSegment<P> {
 // -- CubicCurve
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> Curve<P> for CubicCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> Curve<P> for CubicCurve<P> {
     #[inline]
     fn domain(&self) -> Interval {
         // The non-emptiness invariant guarantees that this succeeds.
@@ -61,7 +61,7 @@ impl<P: VectorSpace> Curve<P> for CubicCurve<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> SampleDerivative<P> for CubicCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleDerivative<P> for CubicCurve<P> {
     #[inline]
     fn sample_with_derivative_unchecked(&self, t: f32) -> WithDerivative<P> {
         WithDerivative {
@@ -72,7 +72,7 @@ impl<P: VectorSpace> SampleDerivative<P> for CubicCurve<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> SampleTwoDerivatives<P> for CubicCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleTwoDerivatives<P> for CubicCurve<P> {
     #[inline]
     fn sample_with_two_derivatives_unchecked(&self, t: f32) -> WithTwoDerivatives<P> {
         WithTwoDerivatives {
@@ -85,7 +85,7 @@ impl<P: VectorSpace> SampleTwoDerivatives<P> for CubicCurve<P> {
 
 // -- RationalSegment
 
-impl<P: VectorSpace> Curve<P> for RationalSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> Curve<P> for RationalSegment<P> {
     #[inline]
     fn domain(&self) -> Interval {
         Interval::UNIT
@@ -97,7 +97,7 @@ impl<P: VectorSpace> Curve<P> for RationalSegment<P> {
     }
 }
 
-impl<P: VectorSpace> SampleDerivative<P> for RationalSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleDerivative<P> for RationalSegment<P> {
     #[inline]
     fn sample_with_derivative_unchecked(&self, t: f32) -> WithDerivative<P> {
         WithDerivative {
@@ -107,7 +107,7 @@ impl<P: VectorSpace> SampleDerivative<P> for RationalSegment<P> {
     }
 }
 
-impl<P: VectorSpace> SampleTwoDerivatives<P> for RationalSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleTwoDerivatives<P> for RationalSegment<P> {
     #[inline]
     fn sample_with_two_derivatives_unchecked(&self, t: f32) -> WithTwoDerivatives<P> {
         WithTwoDerivatives {
@@ -121,7 +121,7 @@ impl<P: VectorSpace> SampleTwoDerivatives<P> for RationalSegment<P> {
 // -- RationalCurve
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> Curve<P> for RationalCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> Curve<P> for RationalCurve<P> {
     #[inline]
     fn domain(&self) -> Interval {
         // The non-emptiness invariant guarantees the success of this.
@@ -136,7 +136,7 @@ impl<P: VectorSpace> Curve<P> for RationalCurve<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> SampleDerivative<P> for RationalCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleDerivative<P> for RationalCurve<P> {
     #[inline]
     fn sample_with_derivative_unchecked(&self, t: f32) -> WithDerivative<P> {
         WithDerivative {
@@ -147,7 +147,7 @@ impl<P: VectorSpace> SampleDerivative<P> for RationalCurve<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> SampleTwoDerivatives<P> for RationalCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> SampleTwoDerivatives<P> for RationalCurve<P> {
     #[inline]
     fn sample_with_two_derivatives_unchecked(&self, t: f32) -> WithTwoDerivatives<P> {
         WithTwoDerivatives {

crates/bevy_math/src/cubic_splines/mod.rs

@@ -68,7 +68,7 @@ impl<P: VectorSpace> CubicBezier<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CubicGenerator<P> for CubicBezier<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicGenerator<P> for CubicBezier<P> {
     type Error = CubicBezierError;
 
     #[inline]
@@ -176,7 +176,7 @@ impl<P: VectorSpace> CubicHermite<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CubicGenerator<P> for CubicHermite<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicGenerator<P> for CubicHermite<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -202,7 +202,7 @@ impl<P: VectorSpace> CubicGenerator<P> for CubicHermite<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CyclicCubicGenerator<P> for CubicHermite<P> {
+impl<P: VectorSpace<Scalar = f32>> CyclicCubicGenerator<P> for CubicHermite<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -313,7 +313,7 @@ impl<P: VectorSpace> CubicCardinalSpline<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CubicGenerator<P> for CubicCardinalSpline<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicGenerator<P> for CubicCardinalSpline<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -351,7 +351,7 @@ impl<P: VectorSpace> CubicGenerator<P> for CubicCardinalSpline<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CyclicCubicGenerator<P> for CubicCardinalSpline<P> {
+impl<P: VectorSpace<Scalar = f32>> CyclicCubicGenerator<P> for CubicCardinalSpline<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -471,7 +471,7 @@ impl<P: VectorSpace> CubicBSpline<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CubicGenerator<P> for CubicBSpline<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicGenerator<P> for CubicBSpline<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -494,7 +494,7 @@ impl<P: VectorSpace> CubicGenerator<P> for CubicBSpline<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CyclicCubicGenerator<P> for CubicBSpline<P> {
+impl<P: VectorSpace<Scalar = f32>> CyclicCubicGenerator<P> for CubicBSpline<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -620,7 +620,7 @@ pub struct CubicNurbs<P: VectorSpace> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CubicNurbs<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicNurbs<P> {
     /// Build a Non-Uniform Rational B-Spline.
     ///
     /// If provided, weights must be the same length as the control points. Defaults to equal weights.
@@ -781,7 +781,7 @@ impl<P: VectorSpace> CubicNurbs<P> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> RationalGenerator<P> for CubicNurbs<P> {
+impl<P: VectorSpace<Scalar = f32>> RationalGenerator<P> for CubicNurbs<P> {
     type Error = InsufficientDataError;
 
     #[inline]
@@ -962,7 +962,7 @@ pub struct CubicSegment<P: VectorSpace> {
     pub coeff: [P; 4],
 }
 
-impl<P: VectorSpace> CubicSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicSegment<P> {
     /// Instantaneous position of a point at parametric value `t`.
     #[inline]
     pub fn position(&self, t: f32) -> P {
@@ -1184,7 +1184,7 @@ pub struct CubicCurve<P: VectorSpace> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> CubicCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> CubicCurve<P> {
     /// Create a new curve from a collection of segments. If the collection of segments is empty,
     /// a curve cannot be built and `None` will be returned instead.
     pub fn from_segments(segments: impl IntoIterator<Item = CubicSegment<P>>) -> Option<Self> {
@@ -1347,7 +1347,7 @@ pub struct RationalSegment<P: VectorSpace> {
     /// The width of the domain of this segment.
     pub knot_span: f32,
 }
-impl<P: VectorSpace> RationalSegment<P> {
+impl<P: VectorSpace<Scalar = f32>> RationalSegment<P> {
     /// Instantaneous position of a point at parametric value `t` in `[0, 1]`.
     #[inline]
     pub fn position(&self, t: f32) -> P {
@@ -1484,7 +1484,7 @@ pub struct RationalCurve<P: VectorSpace> {
 }
 
 #[cfg(feature = "alloc")]
-impl<P: VectorSpace> RationalCurve<P> {
+impl<P: VectorSpace<Scalar = f32>> RationalCurve<P> {
     /// Create a new curve from a collection of segments. If the collection of segments is empty,
     /// a curve cannot be built and `None` will be returned instead.
     pub fn from_segments(segments: impl IntoIterator<Item = RationalSegment<P>>) -> Option<Self> {

crates/bevy_math/src/curve/derivatives/adaptor_impls.rs

@@ -208,10 +208,12 @@ where
 
 // -- ZipCurve
 
-impl<S, T, C, D> SampleDerivative<(S, T)> for ZipCurve<S, T, C, D>
+impl<U, V, S, T, C, D> SampleDerivative<(S, T)> for ZipCurve<S, T, C, D>
 where
-    S: HasTangent,
+    U: VectorSpace<Scalar = f32>,
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    S: HasTangent<Tangent = U>,
+    T: HasTangent<Tangent = V>,
     C: SampleDerivative<S>,
     D: SampleDerivative<T>,
 {
@@ -225,10 +227,12 @@ where
     }
 }
 
-impl<S, T, C, D> SampleTwoDerivatives<(S, T)> for ZipCurve<S, T, C, D>
+impl<U, V, S, T, C, D> SampleTwoDerivatives<(S, T)> for ZipCurve<S, T, C, D>
 where
-    S: HasTangent,
+    U: VectorSpace<Scalar = f32>,
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    S: HasTangent<Tangent = U>,
+    T: HasTangent<Tangent = V>,
     C: SampleTwoDerivatives<S>,
     D: SampleTwoDerivatives<T>,
 {
@@ -248,9 +252,10 @@ where
 
 // -- GraphCurve
 
-impl<T, C> SampleDerivative<(f32, T)> for GraphCurve<T, C>
+impl<V, T, C> SampleDerivative<(f32, T)> for GraphCurve<T, C>
 where
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    T: HasTangent<Tangent = V>,
     C: SampleDerivative<T>,
 {
     fn sample_with_derivative_unchecked(&self, t: f32) -> WithDerivative<(f32, T)> {
@@ -262,9 +267,10 @@ where
     }
 }
 
-impl<T, C> SampleTwoDerivatives<(f32, T)> for GraphCurve<T, C>
+impl<V, T, C> SampleTwoDerivatives<(f32, T)> for GraphCurve<T, C>
 where
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    T: HasTangent<Tangent = V>,
     C: SampleTwoDerivatives<T>,
 {
     fn sample_with_two_derivatives_unchecked(&self, t: f32) -> WithTwoDerivatives<(f32, T)> {
@@ -321,9 +327,10 @@ where
 
 // -- CurveReparamCurve
 
-impl<T, C, D> SampleDerivative<T> for CurveReparamCurve<T, C, D>
+impl<V, T, C, D> SampleDerivative<T> for CurveReparamCurve<T, C, D>
 where
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    T: HasTangent<Tangent = V>,
     C: SampleDerivative<T>,
     D: SampleDerivative<f32>,
 {
@@ -349,9 +356,10 @@ where
     }
 }
 
-impl<T, C, D> SampleTwoDerivatives<T> for CurveReparamCurve<T, C, D>
+impl<V, T, C, D> SampleTwoDerivatives<T> for CurveReparamCurve<T, C, D>
 where
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    T: HasTangent<Tangent = V>,
     C: SampleTwoDerivatives<T>,
     D: SampleTwoDerivatives<f32>,
 {
@@ -386,9 +394,10 @@ where
 
 // -- LinearReparamCurve
 
-impl<T, C> SampleDerivative<T> for LinearReparamCurve<T, C>
+impl<V, T, C> SampleDerivative<T> for LinearReparamCurve<T, C>
 where
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    T: HasTangent<Tangent = V>,
     C: SampleDerivative<T>,
 {
     fn sample_with_derivative_unchecked(&self, t: f32) -> WithDerivative<T> {
@@ -413,9 +422,10 @@ where
     }
 }
 
-impl<T, C> SampleTwoDerivatives<T> for LinearReparamCurve<T, C>
+impl<V, T, C> SampleTwoDerivatives<T> for LinearReparamCurve<T, C>
 where
-    T: HasTangent,
+    V: VectorSpace<Scalar = f32>,
+    T: HasTangent<Tangent = V>,
     C: SampleTwoDerivatives<T>,
 {
     fn sample_with_two_derivatives_unchecked(&self, t: f32) -> WithTwoDerivatives<T> {

crates/bevy_math/src/curve/easing.rs

@@ -32,7 +32,7 @@ pub trait Ease: Sized {
     fn interpolating_curve_unbounded(start: Self, end: Self) -> impl Curve<Self>;
 }
 
-impl<V: VectorSpace> Ease for V {
+impl<V: VectorSpace<Scalar = f32>> Ease for V {
     fn interpolating_curve_unbounded(start: Self, end: Self) -> impl Curve<Self> {
         FunctionCurve::new(Interval::EVERYWHERE, move |t| V::lerp(start, end, t))
     }

crates/bevy_math/src/sampling/shape_sampling.rs

@@ -40,11 +40,12 @@
 
 use core::f32::consts::{PI, TAU};
 
-use crate::{ops, primitives::*, NormedVectorSpace, Vec2, Vec3};
+use crate::{ops, primitives::*, NormedVectorSpace, ScalarField, Vec2, Vec3};
 use rand::{
     distributions::{Distribution, WeightedIndex},
     Rng,
 };
+use rand_distr::uniform::SampleUniform;
 
 /// Exposes methods to uniformly sample a variety of primitive shapes.
 pub trait ShapeSample {
@@ -281,22 +282,24 @@ impl ShapeSample for Cuboid {
 }
 
 /// Interior sampling for triangles which doesn't depend on the ambient dimension.
-fn sample_triangle_interior<P: NormedVectorSpace, R: Rng + ?Sized>(
+fn sample_triangle_interior<P, R>(vertices: [P; 3], rng: &mut R) -> P
-    vertices: [P; 3],
+where
-    rng: &mut R,
+    P: NormedVectorSpace,
-) -> P {
+    P::Scalar: SampleUniform + PartialOrd,
+    R: Rng + ?Sized,
+{
     let [a, b, c] = vertices;
     let ab = b - a;
     let ac = c - a;
 
     // Generate random points on a parallelepiped and reflect so that
     // we can use the points that lie outside the triangle
-    let u = rng.gen_range(0.0..=1.0);
+    let u = rng.gen_range(P::Scalar::ZERO..=P::Scalar::ONE);
-    let v = rng.gen_range(0.0..=1.0);
+    let v = rng.gen_range(P::Scalar::ZERO..=P::Scalar::ONE);
 
-    if u + v > 1. {
+    if u + v > P::Scalar::ONE {
-        let u1 = 1. - v;
+        let u1 = P::Scalar::ONE - v;
-        let v1 = 1. - u;
+        let v1 = P::Scalar::ONE - u;
         a + (ab * u1 + ac * v1)
     } else {
         a + (ab * u + ac * v)
@@ -304,16 +307,18 @@ fn sample_triangle_interior<P: NormedVectorSpace, R: Rng + ?Sized>(
 }
 
 /// Boundary sampling for triangles which doesn't depend on the ambient dimension.
-fn sample_triangle_boundary<P: NormedVectorSpace, R: Rng + ?Sized>(
+fn sample_triangle_boundary<P, R>(vertices: [P; 3], rng: &mut R) -> P
-    vertices: [P; 3],
+where
-    rng: &mut R,
+    P: NormedVectorSpace,
-) -> P {
+    P::Scalar: SampleUniform + PartialOrd + for<'a> ::core::ops::AddAssign<&'a P::Scalar>,
+    R: Rng + ?Sized,
+{
     let [a, b, c] = vertices;
     let ab = b - a;
     let ac = c - a;
     let bc = c - b;
 
-    let t = rng.gen_range(0.0..=1.0);
+    let t = rng.gen_range(P::Scalar::ZERO..=P::Scalar::ONE);
 
     if let Ok(dist) = WeightedIndex::new([ab.norm(), ac.norm(), bc.norm()]) {
         match dist.sample(rng) {

examples/animation/color_animation.rs

@@ -3,8 +3,8 @@
 use bevy::{math::VectorSpace, prelude::*};
 
 // We define this trait so we can reuse the same code for multiple color types that may be implemented using curves.
-trait CurveColor: VectorSpace + Into<Color> + Send + Sync + 'static {}
+trait CurveColor: VectorSpace<Scalar = f32> + Into<Color> + Send + Sync + 'static {}
-impl<T: VectorSpace + Into<Color> + Send + Sync + 'static> CurveColor for T {}
+impl<T: VectorSpace<Scalar = f32> + Into<Color> + Send + Sync + 'static> CurveColor for T {}
 
 // We define this trait so we can reuse the same code for multiple color types that may be implemented using mixing.
 trait MixedColor: Mix + Into<Color> + Send + Sync + 'static {}

release-content/migration-guides/scalar-field-on-vector-space.md

@@ -0,0 +1,8 @@
+---
+title: `VectorSpace` implementations
+pull_requests: [19194]
+---
+
+Previously, implementing `VectorSpace` for a type required your type to use or at least interface with `f32`. This made implementing `VectorSpace` for double-precision types (like `DVec3`) less meaningful and useful, requiring lots of casting. `VectorSpace` has a new required associated type `Scalar` that's bounded by a new trait `ScalarField`. `bevy_math` implements this trait for `f64` and `f32` out of the box, and `VectorSpace` is now implemented for `DVec[N]` types.
+
+If you manually implemented `VectorSpace` for any type, you'll need to implement `Scalar` for it. If you were working with single-precision floating-point types and you want the exact behavior from before, set it to `f32`.