From 91e07261bb9b68016fd4af01422d15320eb03ef8 Mon Sep 17 00:00:00 2001
From: Rouven Spreckels <rs@qu1x.dev>
Date: Mon, 25 Mar 2024 20:12:46 +0100
Subject: [PATCH] Recover from numerical instability.

---
 RELEASES.md      |  5 +++++
 src/ball.rs      | 13 ++++++++-----
 src/enclosing.rs | 19 ++++++++++++-------
 3 files changed, 25 insertions(+), 12 deletions(-)

diff --git a/RELEASES.md b/RELEASES.md
index fb0d811..15de415 100644
--- a/RELEASES.md
+++ b/RELEASES.md
@@ -1,3 +1,8 @@
+# Version 0.4.1 (2024-03-26)
+
+  * Leverage move-to-front heuristic to prevent panic due to numerical instability.
+  * Let `Enclosing::contains` panic for infinite points.
+
 # Version 0.4.0 (2024-03-24)
 
   * Lower trait bounds on `Ball` and `Enclosing`.
diff --git a/src/ball.rs b/src/ball.rs
index 72ebe19..b18fcf0 100644
--- a/src/ball.rs
+++ b/src/ball.rs
@@ -34,7 +34,9 @@ where
 {
 	#[inline]
 	fn contains(&self, point: &OPoint<T, D>) -> bool {
-		(point - &self.center).norm_squared() <= self.radius_squared
+		let norm_squared = (point - &self.center).norm_squared();
+		assert!(norm_squared.is_finite());
+		norm_squared <= self.radius_squared
 	}
 	fn with_bounds(bounds: &[OPoint<T, D>]) -> Option<Self>
 	where
@@ -65,16 +67,17 @@ where
 			}
 		});
 		let vector = vector.view((0, 0), (length, 1));
-		matrix.try_inverse().map(|matrix| {
+		matrix.try_inverse().and_then(|matrix| {
 			let vector = matrix * vector;
 			let mut center = OVector::<T, D>::zeros();
 			for point in 0..length {
 				center += points.column(point) * vector[point].clone();
 			}
-			Self {
+			let radius_squared = center.norm_squared();
+			radius_squared.is_finite().then(|| Self {
 				center: &bounds[0] + &center,
-				radius_squared: center.norm_squared(),
-			}
+				radius_squared,
+			})
 		})
 	}
 }
diff --git a/src/enclosing.rs b/src/enclosing.rs
index a442b6a..9133eb9 100644
--- a/src/enclosing.rs
+++ b/src/enclosing.rs
@@ -147,13 +147,18 @@ where
 		DefaultAllocator: Allocator<T, D, D> + Allocator<OPoint<T, D>, DimNameSum<D, U1>>,
 		<DefaultAllocator as Allocator<OPoint<T, D>, DimNameSum<D, U1>>>::Buffer: Default,
 	{
-		maybe_grow(Self::RED_ZONE, Self::STACK_SIZE, || {
-			Self::enclosing_points_with_bounds(
-				points,
-				&mut OVec::<OPoint<T, D>, DimNameSum<D, U1>>::new(),
-			)
-			.expect("Empty point set")
-		})
+		assert!(!points.is_empty(), "empty point set");
+		for _ in 0..points.len() {
+			if let Some(ball) = maybe_grow(Self::RED_ZONE, Self::STACK_SIZE, || {
+				Self::enclosing_points_with_bounds(
+					points,
+					&mut OVec::<OPoint<T, D>, DimNameSum<D, U1>>::new(),
+				)
+			}) {
+				return ball;
+			}
+		}
+		unreachable!("numerical instability");
 	}
 	/// Returns minimum ball enclosing `points` with `bounds`.
 	///