diff --git a/src/pyrecest/distributions/hypertorus/hypertoroidal_uniform_distribution.py b/src/pyrecest/distributions/hypertorus/hypertoroidal_uniform_distribution.py
index cc83343fc..7d59dd22a 100644
--- a/src/pyrecest/distributions/hypertorus/hypertoroidal_uniform_distribution.py
+++ b/src/pyrecest/distributions/hypertorus/hypertoroidal_uniform_distribution.py
@@ -3,6 +3,7 @@
 # pylint: disable=no-name-in-module,no-member
 import numpy as np
 from pyrecest.backend import (
+    all as backend_all,
     asarray,
     int32,
     int64,
@@ -113,6 +114,11 @@ def _validate_boundary(name: str, value, dim: int):
     return value
 
 
+def _validate_boundary_order(left, right) -> None:
+    if not bool(backend_all(right >= left)):
+        raise ValueError("integration boundaries must be increasing in every dimension")
+
+
 class HypertoroidalUniformDistribution(
     AbstractUniformDistribution, AbstractHypertoroidalDistribution
 ):
@@ -198,6 +204,7 @@ def integrate(self, integration_boundaries=None) -> float:
             left, right = integration_boundaries
         left = _validate_boundary("left", left, self.dim)
         right = _validate_boundary("right", right, self.dim)
+        _validate_boundary_order(left, right)
 
         volume = prod(right - left)
         return 1.0 / (2.0 * pi) ** self.dim * volume
diff --git a/tests/distributions/test_hypertoroidal_uniform_distribution.py b/tests/distributions/test_hypertoroidal_uniform_distribution.py
index 2f446723b..57da29e05 100644
--- a/tests/distributions/test_hypertoroidal_uniform_distribution.py
+++ b/tests/distributions/test_hypertoroidal_uniform_distribution.py
@@ -56,6 +56,13 @@ def test_integrate_validates_boundary_shapes():
         dist.integrate((zeros((2,)), ones((1,))))
 
 
+def test_integrate_rejects_reversed_boundaries():
+    dist = HypertoroidalUniformDistribution(2)
+
+    with pytest.raises(ValueError, match="increasing"):
+        dist.integrate((array([0.0, 1.0]), array([1.0, 0.5])))
+
+
 def test_integrate_accepts_scalar_boundaries_for_one_dimension():
     dist = HypertoroidalUniformDistribution(1)