Removed simplify

.github/CHANGELOG.md

@@ -21,6 +21,7 @@
 * The `README` has been improved to describe the expanding capabilities of `haarpy` [(#39)](https://github.com/polyquantique/haarpy/pull/39).
 * Added references and examples to the docstrings. Slight modification to the docstring format  [(#47)](https://github.com/polyquantique/haarpy/pull/47).
 * `get_conjugacy_class()` has been sped up [(#51)](https://github.com/polyquantique/haarpy/pull/51).
+* `weingarten_...()` and `haar_integral_...()` have been sped up by building logic to simplify rational functions in `utils._simplify`  [(#53)](https://github.com/polyquantique/haarpy/pull/53).
 
 ### Bug fixes
 

.gitignore

@@ -1,5 +1,6 @@
 __pycache__/
 _build/
+notebooks/
 
 .coverage*
 coverage.xml

haarpy/_utils.py

@@ -0,0 +1,62 @@
+# Copyright 2026 Polyquantique
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#    http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Utility Python interface
+"""
+
+from collections.abc import Iterable
+from fractions import Fraction
+from sympy import Add, Mul, together, fraction, factor, factor_list, Expr, PolificationFailed
+
+
+def _simplify(expr_iter: Iterable[Expr], constant: Fraction = Fraction(1, 1)) -> Expr:
+    """Factorizes a sum of rational fraction into a
+    single factorized, simplified fraction
+
+    Parameters
+    ----------
+        expr (Expr) : the expression to be simplified
+        constant (Fraction) : A constant fraction multiplying the simplification result
+
+    Returns
+    -------
+        Expr : the simplified fraction
+    """
+    equation = together(Add(*expr_iter))
+    #automatically factorises the denominator
+    num, denum = fraction(constant * equation)
+
+    # Error occurs if numerator is a constant
+    try:
+        num_factor_list = factor_list(num)
+    except PolificationFailed:
+        return factor(num) / denum
+
+    denum_factor_list = factor_list(denum)
+
+    # gets rid of common factors
+    num_factor_dict, denum_factor_dict = dict(num_factor_list[1]), dict(denum_factor_list[1])
+    for fact in num_factor_dict:
+        if fact in denum_factor_dict:
+            common = min(num_factor_dict[fact], denum_factor_dict[fact])
+            num_factor_dict[fact] -= common
+            denum_factor_dict[fact] -= common
+
+    num_simplified = Mul(*(fact**power for fact, power in num_factor_dict.items()))
+    denum_simplified = Mul(*(fact**power for fact, power in denum_factor_dict.items()))
+    constant_simplified = Fraction(num_factor_list[0], denum_factor_list[0])
+
+    return Mul(constant_simplified.numerator, num_simplified, evaluate=False) / Mul(
+        constant_simplified.denominator, denum_simplified, evaluate=False
+    )

haarpy/circular_ensembles.py

@@ -25,7 +25,7 @@
 from functools import lru_cache
 from typing import Union
 from collections import Counter
-from sympy import Symbol, Expr, fraction, factor, simplify
+from sympy import Symbol, Expr
 from sympy.combinatorics import Permutation, SymmetricGroup
 from sympy.core.numbers import Integer
 from haarpy import (
@@ -34,6 +34,7 @@
     coset_type,
     stabilizer_coset,
 )
+from ._utils import _simplify
 
 
 @lru_cache
@@ -158,16 +159,12 @@ def haar_integral_circular_orthogonal(
         coset_type(permutation) for permutation in stabilizer_coset(seq_i, seq_j)
     )
 
-    integral = sum(
+    integral_gen = (
         count * weingarten_circular_orthogonal(coset, group_dimension)
         for coset, count in coset_mapping.items()
     )
 
-    if isinstance(group_dimension, Expr):
-        numerator, denominator = fraction(simplify(integral))
-        integral = factor(numerator) / factor(denominator)
-
-    return integral
+    return sum(integral_gen) if isinstance(group_dimension, int) else _simplify(integral_gen)
 
 
 @lru_cache
@@ -298,13 +295,13 @@ def haar_integral_circular_symplectic(sequences: tuple[tuple[Expr]], half_dimens
         if permutation(shifted_i) == shifted_j
     )
 
-    integral = coefficient * sum(
+    integral_gen = (
         weingarten_circular_symplectic(permutation, half_dimension)
         for permutation in permutation_tuple
     )
 
-    if isinstance(half_dimension, Expr):
-        numerator, denominator = fraction(simplify(integral))
-        integral = factor(numerator) / factor(denominator)
-
-    return integral
+    return (
+        coefficient * sum(integral_gen)
+        if isinstance(half_dimension, int)
+        else _simplify(integral_gen, Fraction(coefficient, 1))
+    )

haarpy/orthogonal.py

@@ -29,7 +29,7 @@
 from functools import lru_cache
 from typing import Union
 from collections import Counter
-from sympy import Symbol, Expr, factorial, factor, fraction, simplify
+from sympy import Symbol, Expr, factorial
 from sympy.combinatorics import Permutation
 from sympy.utilities.iterables import partitions
 from haarpy import (
@@ -41,6 +41,7 @@
     coset_type,
     coset_type_representative,
 )
+from ._utils import _simplify
 
 
 @lru_cache
@@ -184,20 +185,16 @@ def weingarten_orthogonal(
             factorial(degree),
         )
     else:
-        weingarten = (
-            sum(
-                irrep_dim * zonal_spherical / coefficient
-                for irrep_dim, zonal_spherical, coefficient in zip(
-                    irrep_dimension_gen, zonal_spherical_gen, coefficient_gen
-                )
-                if coefficient
+        weingarten_gen = (
+            irrep_dim * zonal_spherical / coefficient
+            for irrep_dim, zonal_spherical, coefficient in zip(
+                irrep_dimension_gen, zonal_spherical_gen, coefficient_gen
             )
-            * 2**half_degree
-            * factorial(half_degree)
-            / factorial(degree)
+            if coefficient
+        )
+        weingarten = _simplify(
+            weingarten_gen, Fraction(2**half_degree * factorial(half_degree), factorial(degree))
         )
-        numerator, denominator = fraction(simplify(weingarten))
-        weingarten = simplify(factor(numerator) / factor(denominator))
 
     return weingarten
 
@@ -263,13 +260,9 @@ def haar_integral_orthogonal(sequences: tuple[tuple[int]], orthogonal_dimension:
         coset_type(cycle_j * ~cycle_i) for cycle_i, cycle_j in product(permutation_i, permutation_j)
     )
 
-    integral = sum(
+    integral_gen = (
         count * weingarten_orthogonal(coset, orthogonal_dimension)
         for coset, count in coset_mapping.items()
     )
 
-    if isinstance(orthogonal_dimension, Expr):
-        numerator, denominator = fraction(simplify(integral))
-        integral = factor(numerator) / factor(denominator)
-
-    return integral
+    return sum(integral_gen) if isinstance(orthogonal_dimension, int) else _simplify(integral_gen)

haarpy/permutation.py

@@ -25,8 +25,9 @@
 from itertools import product
 from collections.abc import Sequence
 from fractions import Fraction
-from sympy import Symbol, simplify, binomial, factor, fraction
+from sympy import Symbol, binomial
 from haarpy import set_partitions, meet_operation, join_operation, partial_order
+from ._utils import _simplify
 
 
 @lru_cache
@@ -117,14 +118,13 @@ def weingarten_permutation(
             for partition in inferieur_partition_tuple
         )
     else:
-        weingarten = sum(
+        weingarten_gen = (
             mobius_function(partition, first_partition)
             * mobius_function(partition, second_partition)
             / prod(dimension - i for i, _ in enumerate(partition))
             for partition in inferieur_partition_tuple
         )
-        numerator, denominator = fraction(simplify(weingarten))
-        weingarten = factor(numerator) / factor(denominator)
+        weingarten = _simplify(weingarten_gen)
 
     return weingarten
 
@@ -192,7 +192,7 @@ def weingarten_centered_permutation(
             for i in range(singleton_set_size + 1)
         )
     else:
-        weingarten = sum(
+        weingarten_gen = (
             (-1) ** i
             * binomial(singleton_set_size, i)
             * sum(
@@ -205,8 +205,7 @@ def weingarten_centered_permutation(
             for i in range(singleton_set_size + 1)
         )
 
-        num, denum = fraction(simplify(weingarten))
-        weingarten = factor(num) / factor(denum)
+        weingarten = _simplify(weingarten_gen)
 
     return weingarten
 
@@ -321,7 +320,7 @@ def haar_integral_centered_permutation(
         for unique in set(column_indices)
     )
 
-    integral = sum(
+    integral_gen = (
         weingarten_centered_permutation(
             partition_sigma,
             partition_tau,
@@ -335,8 +334,4 @@ def haar_integral_centered_permutation(
         and partial_order(partition_tau, column_partition)
     )
 
-    if isinstance(dimension, Symbol):
-        num, denum = fraction(simplify(integral))
-        integral = factor(num) / factor(denum)
-
-    return integral
+    return sum(integral_gen) if isinstance(dimension, int) else _simplify(integral_gen)

haarpy/symplectic.py

@@ -27,7 +27,7 @@
 from fractions import Fraction
 from functools import lru_cache
 from itertools import product
-from sympy import Symbol, factorial, factor, fraction, simplify, Expr
+from sympy import Symbol, factorial, Expr
 from sympy.combinatorics import Permutation
 from sympy.utilities.iterables import partitions
 from sympy.core.numbers import Integer
@@ -38,6 +38,7 @@
     irrep_dimension,
     hyperoctahedral_transversal,
 )
+from ._utils import _simplify
 
 
 @lru_cache
@@ -165,20 +166,16 @@ def weingarten_symplectic(permutation: Permutation, half_dimension: Symbol) -> E
             factorial(degree),
         )
     else:
-        weingarten = (
-            sum(
-                irrep_dim * zonal_spherical / coefficient
-                for irrep_dim, zonal_spherical, coefficient in zip(
-                    irrep_dimension_gen, twisted_spherical_gen, coefficient_gen
-                )
-                if coefficient
+        weingarten_gen = (
+            irrep_dim * zonal_spherical / coefficient
+            for irrep_dim, zonal_spherical, coefficient in zip(
+                irrep_dimension_gen, twisted_spherical_gen, coefficient_gen
             )
-            * 2**half_degree
-            * factorial(half_degree)
-            / factorial(degree)
+            if coefficient
+        )
+        weingarten = (
+            _simplify(weingarten_gen) * 2**half_degree * factorial(half_degree) / factorial(degree)
         )
-        numerator, denominator = fraction(simplify(weingarten))
-        weingarten = simplify(factor(numerator) / factor(denominator))
 
     return weingarten
 
@@ -296,14 +293,10 @@ def twisted_delta(seq_value, seq_pos, perm):
         for perm in hyperoctahedral_transversal(degree)
     )
 
-    integral = sum(
+    integral_gen = (
         perm_i[1] * perm_j[1] * weingarten_symplectic(perm_j[0] * ~perm_i[0], half_dimension)
         for perm_i, perm_j in product(permutation_i_tuple, permutation_j_tuple)
         if perm_i[1] * perm_j[1]
     )
 
-    if isinstance(half_dimension, Expr):
-        numerator, denominator = fraction(simplify(integral))
-        integral = factor(numerator) / factor(denominator)
-
-    return integral
+    return sum(integral_gen) if isinstance(half_dimension, int) else _simplify(integral_gen)

haarpy/tests/test_circular_ensembles.py

@@ -56,7 +56,7 @@ def test_weingarten_circular_orthogonal_hyperoctahedral_symbolic(half_degree):
         for permutation in SymmetricGroup(2 * half_degree).generate():
             hyperoctahedral = ap.HyperoctahedralGroup(half_degree)
             coefficient = 1 / (d * (d + 1) * (d + 3))
-            assert ap.weingarten_circular_orthogonal(permutation, d) == (
+            assert ap.weingarten_circular_orthogonal(permutation, d).equals(
                 simplify((d + 2) * coefficient) if permutation in hyperoctahedral else -coefficient
             )
 
@@ -94,7 +94,7 @@ def test_weingarten_circular_symplectic_hyperoctahedral_symbolic(half_degree):
         for permutation in SymmetricGroup(2 * half_degree).generate():
             hyperoctahedral = ap.HyperoctahedralGroup(half_degree)
             coefficient = permutation.signature() / (d * (2 * d - 1) * (2 * d - 3))
-            assert ap.weingarten_circular_symplectic(permutation, d) == (
+            assert ap.weingarten_circular_symplectic(permutation, d).equals(
                 simplify((d - 1) * coefficient)
                 if permutation in hyperoctahedral
                 else coefficient / 2

haarpy/tests/test_orthogonal.py

@@ -223,7 +223,7 @@ def test_zonal_spherical_type_error(permutation, partition):
         (Permutation(0, 1, 2, 3, 4, 5), 2, d * (d - 1) * (d - 2) * (d + 2) * (d + 4)),
         (
             Permutation(0, 1, 2, 3, 4, 5, 6, 7),
-            -5 * d - 6,
+            -(5 * d + 6),
             d * (d - 1) * (d - 2) * (d - 3) * (d + 1) * (d + 2) * (d + 4) * (d + 6),
         ),
         (
@@ -238,7 +238,7 @@ def test_zonal_spherical_type_error(permutation, partition):
         ),
         (
             Permutation(4, 5, 6, 7),
-            -(d**3) - 6 * d**2 - 3 * d + 6,
+            -(d**3 + 6 * d**2 + 3 * d - 6),
             d * (d - 1) * (d - 2) * (d - 3) * (d + 1) * (d + 2) * (d + 4) * (d + 6),
         ),
         (
@@ -253,7 +253,7 @@ def test_weingarten_orthogonal_literature(permutation, num, denum):
     `Collins et al. Integration with Respect to the Haar Measure on Unitary, Orthogonal
     and Symplectic Group: <https://link.springer.com/article/10.1007/s00220-006-1554-3>`_
     """
-    assert ap.weingarten_orthogonal(permutation, d) == num / denum
+    assert ap.weingarten_orthogonal(permutation, d).equals(num / denum)
 
 
 @pytest.mark.parametrize("degree", range(2, 10, 2))

haarpy/tests/test_permutation.py

@@ -136,7 +136,7 @@ def test_weingarten_centered_permutation_hand_calculated(partition1, partition2,
     dimension = randint(10, 99)
     assert ap.weingarten_centered_permutation(
         partition1, partition2, d
-    ) == hand_calculated_weingarten[result_key] and ap.weingarten_centered_permutation(
+    ).equals(hand_calculated_weingarten[result_key]) and ap.weingarten_centered_permutation(
         partition1, partition2, dimension
     ) == Fraction(
         hand_calculated_weingarten[result_key].subs(d, dimension)
@@ -271,9 +271,9 @@ def test_haar_integral_centered_permutation_weingarten(row_indices, column_indic
         sum(value * hand_calculated_weingarten[key] for key, value in result_dict.items())
     )
     num, denum = fraction(result)
-    assert ap.haar_integral_centered_permutation(row_indices, column_indices, d) == factor(
-        num
-    ) / factor(denum)
+    assert ap.haar_integral_centered_permutation(row_indices, column_indices, d).equals(
+        factor(num) / factor(denum)
+    )
 
 
 @pytest.mark.parametrize(