From 1ac4855541e0c75f2c33ee719c94065f879fdd82 Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Tue, 25 Aug 2026 14:05:53 +0800
Subject: [PATCH 1/2] Fix high-dimensional ellipsoid volume overflow

---
 .../abstract_ellipsoidal_ball_distribution.py    | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/pyrecest/distributions/abstract_ellipsoidal_ball_distribution.py b/src/pyrecest/distributions/abstract_ellipsoidal_ball_distribution.py
index 928b7f741..8e34d530e 100644
--- a/src/pyrecest/distributions/abstract_ellipsoidal_ball_distribution.py
+++ b/src/pyrecest/distributions/abstract_ellipsoidal_ball_distribution.py
@@ -1,3 +1,5 @@
+from math import lgamma, log as scalar_log, pi as scalar_pi
+
 # pylint: disable=no-name-in-module,no-member
 from pyrecest.backend import all as backend_all
 from pyrecest.backend import (
@@ -14,7 +16,6 @@
     transpose,
 )
 from pyrecest.exceptions import ShapeError, ValidationError
-from scipy.special import gamma
 
 from .abstract_bounded_nonperiodic_distribution import (
     AbstractBoundedNonPeriodicDistribution,
@@ -97,8 +98,17 @@ def get_manifold_size(self):
         elif self.dim == 4:
             c = 0.5 * pi**2
         else:
-            c = (pi ** (self.dim / 2)) / gamma((self.dim / 2) + 1)
+            half_dim = self.dim / 2
+            # Keep the unit-ball volume in log space. gamma(half_dim + 1)
+            # overflows starting around dim=342 even while the volume itself
+            # is still representable in float64.
+            log_unit_ball_volume = half_dim * scalar_log(scalar_pi) - lgamma(
+                half_dim + 1
+            )
+
+        if self.dim <= 4:
+            log_unit_ball_volume = log(array(c))
 
         cholesky_factor = linalg.cholesky(self.shape_matrix)
-        log_volume = log(array(c)) + log(diagonal(cholesky_factor)).sum()
+        log_volume = log_unit_ball_volume + log(diagonal(cholesky_factor)).sum()
         return exp(log_volume)

From b4aed678a9544aa163fe515b28650e8db279218e Mon Sep 17 00:00:00 2001
From: Florian Pfaff <6773539+FlorianPfaff@users.noreply.github.com>
Date: Tue, 25 Aug 2026 14:06:14 +0800
Subject: [PATCH 2/2] Add regression test for ellipsoid volume overflow

---
 ...llipsoidal_ball_high_dimensional_volume.py | 35 +++++++++++++++++++
 1 file changed, 35 insertions(+)
 create mode 100644 tests/distributions/test_ellipsoidal_ball_high_dimensional_volume.py

diff --git a/tests/distributions/test_ellipsoidal_ball_high_dimensional_volume.py b/tests/distributions/test_ellipsoidal_ball_high_dimensional_volume.py
new file mode 100644
index 000000000..34648e63f
--- /dev/null
+++ b/tests/distributions/test_ellipsoidal_ball_high_dimensional_volume.py
@@ -0,0 +1,35 @@
+import unittest
+
+import mpmath
+import numpy as np
+import numpy.testing as npt
+
+# pylint: disable=no-name-in-module,no-member
+from pyrecest.backend import array
+from pyrecest.distributions import EllipsoidalBallUniformDistribution
+
+
+class TestEllipsoidalBallHighDimensionalVolume(unittest.TestCase):
+    def test_volume_and_pdf_remain_finite_after_gamma_overflow_threshold(self):
+        dim = 342
+        dist = EllipsoidalBallUniformDistribution(
+            array(np.zeros(dim)), array(np.eye(dim))
+        )
+
+        with mpmath.workdps(80):
+            expected_volume = float(
+                mpmath.pi ** (dim / 2) / mpmath.gamma(dim / 2 + 1)
+            )
+
+        volume = dist.get_manifold_size()
+        npt.assert_allclose(volume, expected_volume, rtol=1e-12, atol=0.0)
+        npt.assert_allclose(
+            dist.pdf(array(np.zeros(dim))),
+            1.0 / expected_volume,
+            rtol=1e-12,
+            atol=0.0,
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()