Renormalize rotation in 2D integration and clean up renormalize (#590)

# Objective

Currently, only 3D rotation is renormalized after position integration. However, 2D rotation can also become denormalized over time, sometimes causing crashes with `debug_assertions`.

## Solution

Add `fast_renormalize` for 2D `Rotation` and renormalize 2D rotations after position integration.

I also renamed the 3D `Rotation::renormalize` method to `fast_renormalize`, and made it return the value instead of mutating `self`, to match Bevy 0.15.

---

## Migration Guide

The 3D `Rotation::renormalize` method has been renamed to `fast_renormalize`, and it returns the value instead of mutating `self`. This was changed to match Bevy 0.15.

src/collision/collider/backend.rs

@@ -600,9 +600,9 @@ fn update_aabb<C: AnyCollider>(
             }
             #[cfg(feature = "3d")]
             {
-                let mut end_rot =
-                    Rotation(Quaternion::from_scaled_axis(ang_vel.0 * delta_secs) * rot.0);
-                end_rot.renormalize();
+                let end_rot =
+                    Rotation(Quaternion::from_scaled_axis(ang_vel.0 * delta_secs) * rot.0)
+                        .fast_renormalize();
                 (
                     pos.0
                         + (lin_vel.0 * delta_secs)

src/dynamics/ccd/mod.rs

@@ -646,7 +646,7 @@ fn solve_swept_ccd(
             {
                 let delta_rot = Quaternion::from_scaled_axis(ang_vel1.0 * min_toi);
                 rot1.0 = delta_rot * prev_rot.0 .0;
-                rot1.renormalize();
+                *rot1 = rot1.fast_renormalize();
             }
 
             if let Some(mut collider_translation) = body2.translation {
@@ -663,7 +663,7 @@ fn solve_swept_ccd(
                     let delta_rot = Quaternion::from_scaled_axis(collider_ang_vel * min_toi);
 
                     body2.rot.0 = delta_rot * collider_prev_rot.0 .0;
-                    body2.rot.renormalize();
+                    *body2.rot = body2.rot.fast_renormalize();
                 }
             }
         }

src/dynamics/integrator/semi_implicit_euler.rs

@@ -123,6 +123,7 @@ pub fn integrate_position(
         let delta_rot = Rotation::radians(ang_vel * delta_seconds);
         if delta_rot != Rotation::IDENTITY && delta_rot.is_finite() {
             *rot *= delta_rot;
+            *rot = rot.fast_renormalize();
         }
     }
     #[cfg(feature = "3d")]
@@ -133,7 +134,7 @@ pub fn integrate_position(
         if scaled_axis != AngularVelocity::ZERO.0 && scaled_axis.is_finite() {
             let delta_rot = Quaternion::from_scaled_axis(scaled_axis);
             rot.0 = delta_rot * rot.0;
-            rot.renormalize();
+            *rot = rot.fast_renormalize();
         }
     }
 }

src/dynamics/solver/xpbd/angular_constraint.rs

@@ -91,13 +91,13 @@ pub trait AngularConstraint: XpbdConstraint<2> {
             //       Maybe the math above can be done in a way that keeps rotations normalized?
             let delta_quat = Self::get_delta_rot(inv_inertia1, impulse);
             body1.rotation.0 = delta_quat * body1.rotation.0;
-            body1.rotation.renormalize();
+            *body1.rotation = body1.rotation.fast_renormalize();
         }
         if body2.rb.is_dynamic() {
             // See comments for `body1` above.
             let delta_quat = Self::get_delta_rot(inv_inertia2, -impulse);
             body2.rotation.0 = delta_quat * body2.rotation.0;
-            body2.rotation.renormalize();
+            *body2.rotation = body2.rotation.fast_renormalize();
         }
 
         impulse
@@ -238,13 +238,13 @@ pub trait AngularConstraint: XpbdConstraint<2> {
             //       Maybe the math above can be done in a way that keeps rotations normalized?
             let delta_quat = Self::get_delta_rot(inv_inertia1, p);
             body1.rotation.0 = delta_quat * body1.rotation.0;
-            body1.rotation.renormalize();
+            *body1.rotation = body1.rotation.fast_renormalize();
         }
         if body2.rb.is_dynamic() {
             // See comments for `body1` above.
             let delta_quat = Self::get_delta_rot(inv_inertia2, -p);
             body2.rotation.0 = delta_quat * body2.rotation.0;
-            body2.rotation.renormalize();
+            *body2.rotation = body2.rotation.fast_renormalize();
         }
 
         p

src/dynamics/solver/xpbd/positional_constraint.rs

@@ -59,7 +59,7 @@ pub trait PositionConstraint: XpbdConstraint<2> {
                 //       Maybe the math above can be done in a way that keeps rotations normalized?
                 let delta_quat = Self::get_delta_rot(inv_inertia1, r1, impulse);
                 body1.rotation.0 = delta_quat * body1.rotation.0;
-                body1.rotation.renormalize();
+                *body1.rotation = body1.rotation.fast_renormalize();
             }
         }
         if body2.rb.is_dynamic() && body2.dominance() <= body1.dominance() {
@@ -75,7 +75,7 @@ pub trait PositionConstraint: XpbdConstraint<2> {
                 // See comments for `body1` above.
                 let delta_quat = Self::get_delta_rot(inv_inertia2, r2, -impulse);
                 body2.rotation.0 = delta_quat * body2.rotation.0;
-                body2.rotation.renormalize();
+                *body2.rotation = body2.rotation.fast_renormalize();
             }
         }
 

src/position.rs

@@ -332,6 +332,7 @@ impl Rotation {
     /// accumulated floating point error, or if the rotation was constructed
     /// with invalid values.
     #[inline]
+    #[must_use]
     pub fn try_normalize(self) -> Option<Self> {
         let recip = self.length_recip();
         if recip.is_finite() && recip > 0.0 {
@@ -352,11 +353,25 @@ impl Rotation {
     ///
     /// Panics if `self` has a length of zero, NaN, or infinity when debug assertions are enabled.
     #[inline]
+    #[must_use]
     pub fn normalize(self) -> Self {
         let length_recip = self.length_recip();
         Self::from_sin_cos(self.sin * length_recip, self.cos * length_recip)
     }
 
+    /// Returns `self` after an approximate normalization,
+    /// assuming the value is already nearly normalized.
+    /// Useful for preventing numerical error accumulation.
+    #[inline]
+    #[must_use]
+    pub fn fast_renormalize(self) -> Self {
+        // First-order Tayor approximation
+        // 1/L = (L^2)^(-1/2) ≈ 1 - (L^2 - 1) / 2 = (3 - L^2) / 2
+        let length_squared = self.length_squared();
+        let approx_inv_length = 0.5 * (3.0 - length_squared);
+        Self::from_sin_cos(self.sin * approx_inv_length, self.cos * approx_inv_length)
+    }
+
     /// Returns `true` if the rotation is neither infinite nor NaN.
     #[inline]
     pub fn is_finite(self) -> bool {
@@ -687,6 +702,8 @@ pub struct Rotation(pub Quaternion);
 #[cfg(feature = "3d")]
 impl Rotation {
     /// Inverts the rotation.
+    #[inline]
+    #[must_use]
     pub fn inverse(&self) -> Self {
         Self(self.0.inverse())
     }
@@ -702,6 +719,7 @@ impl Rotation {
     /// the result resembles a kind of ease-in-out effect.
     ///
     /// If you would like the angular velocity to remain constant, consider using [`slerp`](Self::slerp) instead.
+    #[inline]
     pub fn nlerp(self, end: Self, t: Scalar) -> Self {
         Self(self.0.lerp(end.0, t))
     }
@@ -716,17 +734,22 @@ impl Rotation {
     ///
     /// If you would like the rotation to have a kind of ease-in-out effect, consider
     /// using the slightly more efficient [`nlerp`](Self::nlerp) instead.
+    #[inline]
     pub fn slerp(self, end: Self, t: Scalar) -> Self {
         Self(self.0.slerp(end.0, t))
     }
 
-    /// Performs a renormalization of the contained quaternion
+    /// Returns `self` after an approximate normalization,
+    /// assuming the value is already nearly normalized.
+    /// Useful for preventing numerical error accumulation.
     #[inline]
-    pub fn renormalize(&mut self) {
-        let length_squared = self.0.length_squared();
-        // 1/L = (L^2)^-(1/2) = ~= 1 - (L^2 - 1) / 2 = (3 - L^2) / 2
+    #[must_use]
+    pub fn fast_renormalize(self) -> Self {
+        // First-order Tayor approximation
+        // 1/L = (L^2)^(-1/2) ≈ 1 - (L^2 - 1) / 2 = (3 - L^2) / 2
+        let length_squared = self.length_squared();
         let approx_inv_length = 0.5 * (3.0 - length_squared);
-        self.0 = self.0 * approx_inv_length;
+        Self(self.0 * approx_inv_length)
     }
 }
 

src/tests/determinism_2d.rs

@@ -58,7 +58,7 @@ fn cross_platform_determinism_2d() {
     let hash = compute_hash(app.world(), query);
 
     // Update this value if simulation behavior changes.
-    let expected = 0x47306cfb;
+    let expected = 0x5b90b194;
 
     assert!(
         hash == expected,