Repository: sigma-py/quadpy (language: python)

Issue:
ndim
Some helper functions have moved over to https://github.com/nschloe/ndim, use them from there.


Discussion / hints:

Fix patch:
diff --git a/src/quadpy/tn/_helpers.py b/src/quadpy/tn/_helpers.py
index db607b47..337ae2c4 100644
--- a/src/quadpy/tn/_helpers.py
+++ b/src/quadpy/tn/_helpers.py
@@ -2,7 +2,6 @@
 from typing import Callable, List, Optional
 
 import numpy as np
-import sympy
 
 from .._exception import QuadpyError
 from ..helpers import QuadratureScheme
@@ -86,41 +85,3 @@ def get_vol(simplex):
 
     vol = np.sqrt((-1.0) ** (j + 1) / 2 ** j / math.factorial(j) ** 2 * det)
     return vol
-
-
-def integrate_monomial_over_unit_simplex(k, symbolic=False):
-    """The integrals of monomials over the standard triangle and tetrahedron are
-    given by
-
-    \\int_T x_0^k0 * x1^k1 = (k0!*k1!) / (2+k0+k1)!,
-    \\int_T x_0^k0 * x1^k1 * x2^k2 = (k0!*k1!*k2!) / (3+k0+k1+k2)!,
-
-    see, e.g.,
-    A set of symmetric quadrature rules on triangles and tetrahedra,
-    Linbo Zhang, Tao Cui and Hui Liu,
-    Journal of Computational Mathematics,
-    Vol. 27, No. 1 (January 2009), pp. 89-96,
-    <https://www.jstor.org/stable/43693493>.
-
-    See, e.g., <https://math.stackexchange.com/q/207073/36678> for a formula in
-    all dimensions.
-
-    To cope with the huge terms in numerator and denominator, it might make sense to use
-    exp(lgamma()). It's even easier though to represent the above expression by the
-    recurrence with the simple factor k_i / (sum(k) + n) which will never overflow. It's
-    also well-suited for symbolic computation.
-    """
-    frac = sympy.Rational if symbolic else lambda a, b: a / b
-
-    assert all(kk >= 0 for kk in k)
-
-    n = len(k)
-    if all(kk == 0 for kk in k):
-        return frac(1, math.factorial(n))
-
-    # find first nonzero
-    idx = next(i for i, j in enumerate(k) if j > 0)
-    alpha = frac(k[idx], sum(k) + n)
-    k2 = k.copy()
-    k2[idx] -= 1
-    return integrate_monomial_over_unit_simplex(k2, symbolic) * alpha
diff --git a/src/quadpy/tn/_stroud_1969.py b/src/quadpy/tn/_stroud_1969.py
index d12206a3..0cee1e0f 100644
--- a/src/quadpy/tn/_stroud_1969.py
+++ b/src/quadpy/tn/_stroud_1969.py
@@ -1,9 +1,10 @@
+import ndim
 import numpy as np
 from sympy import Rational as frac
 from sympy import sqrt
 
 from ..helpers import article, rd, untangle
-from ._helpers import TnScheme, integrate_monomial_over_unit_simplex
+from ._helpers import TnScheme
 
 source = article(
     authors=["A.H. Stroud"],
@@ -45,11 +46,9 @@ def stroud_1969(n):
         pts.append(rd(n + 1, [(u2, n - 1), (v2, 2)]))
         k_range.append(5)
 
-    b0 = integrate_monomial_over_unit_simplex(n * [0], symbolic=True)
+    b0 = ndim.nsimplex.integrate_monomial(n * [0], symbolic=True)
     b = [
-        integrate_monomial_over_unit_simplex(
-            np.array([k] + (n - 1) * [0]), symbolic=True
-        )
+        ndim.nsimplex.integrate_monomial(np.array([k] + (n - 1) * [0]), symbolic=True)
         / b0
         for k in k_range
     ]


Test patch:
diff --git a/tests/test_t3.py b/tests/test_t3.py
index d5c92393..01d9d688 100644
--- a/tests/test_t3.py
+++ b/tests/test_t3.py
@@ -1,10 +1,10 @@
+import ndim
 import numpy as np
 import pytest
 import sympy
 from helpers import check_degree
 
 import quadpy
-from quadpy.tn._helpers import integrate_monomial_over_unit_simplex
 
 
 def _integrate_exact(f, t3):
@@ -57,7 +57,7 @@ def test_scheme(scheme):
 
     degree, err = check_degree(
         lambda poly: scheme.integrate(poly, t3),
-        integrate_monomial_over_unit_simplex,
+        ndim.nsimplex.integrate_monomial,
         3,
         scheme.degree + 1,
         scheme.test_tolerance * 1.1,
diff --git a/tests/test_tn.py b/tests/test_tn.py
index a9e5a522..20a003e5 100644
--- a/tests/test_tn.py
+++ b/tests/test_tn.py
@@ -1,9 +1,9 @@
+import ndim
 import numpy as np
 import pytest
 from helpers import check_degree
 
 import quadpy
-from quadpy.tn._helpers import integrate_monomial_over_unit_simplex
 
 
 @pytest.mark.parametrize(
@@ -56,7 +56,7 @@ def test_scheme(scheme):
         simplex[k + 1, k] = 1
     degree, err = check_degree(
         lambda poly: scheme.integrate(poly, simplex),
-        integrate_monomial_over_unit_simplex,
+        ndim.nsimplex.integrate_monomial,
         n,
         scheme.degree + 1,
         tol=scheme.test_tolerance,