fix: avoid SymPy 1.13 simplification slowdown

doc/changelog.rst

@@ -15,11 +15,12 @@ Changelog
   time out after 600 s on SymPy ≥ 1.13.  SymPy PR #26390 added an O(N·M)
   ``.replace()`` traversal inside ``TR3``/``futrig`` that is a no-op for
   galgebra's symbolic trig arguments but dominated each of the ~70
-  ``Simp.apply`` calls during ``Ga.build(norm=True)`` for curvilinear
-  coordinates.  The fix uses ``trigsimp(method='old')`` via ``Simp.profile``
-  for the affected example, cutting run time from > 600 s to < 6 s.
-  A notebook note documents the two cosmetic output differences from the
-  pre-1.13 form; a proper upstream fix is tracked in :issue:`576`.
+  ``Simp.apply`` calls for large curvilinear-coordinate expressions.
+  ``Simp`` now detects expressions likely to trigger that traversal and uses
+  ``trigsimp(method='old')`` for those expressions while retaining
+  ``simplify`` for smaller expressions.  Explicit ``Simp.profile`` settings
+  continue to replace the default simplifier.  This library-level fallback
+  also removes the need for an example-wide profile override.
 
 - :support:`589` Added Step 0 to the release-process runbook
   (``doc/dev/release-process.md``): open a release issue before preparing the

examples/LaTeX/curvi_linear_latex.py

@@ -182,26 +182,14 @@ def main():
     #Eprint()
     Format()
 
-    # SymPy >= 1.13 (PR #26390) added a slow O(N*M) traversal inside
-    # sympy.simplify.fu that causes timeouts on curvilinear coordinate
-    # expressions.  Use trigsimp(method='old') via Simp.profile to avoid
-    # that code path entirely for this example.
-    from sympy import trigsimp
-    from galgebra.metric import Simp
-
-    orig_modes = Simp.modes[:]
-    Simp.profile([lambda e: trigsimp(e, method='old')])
-    try:
-        derivatives_in_spherical_coordinates()
-        derivatives_in_paraboloidal_coordinates()
-        # FIXME This takes ~600 seconds
-        # derivatives_in_elliptic_cylindrical_coordinates()
-        derivatives_in_prolate_spheroidal_coordinates()
-        #derivatives_in_oblate_spheroidal_coordinates()
-        #derivatives_in_bipolar_coordinates()
-        #derivatives_in_toroidal_coordinates()
-    finally:
-        Simp.profile(orig_modes)
+    derivatives_in_spherical_coordinates()
+    derivatives_in_paraboloidal_coordinates()
+    # FIXME This takes ~600 seconds
+    # derivatives_in_elliptic_cylindrical_coordinates()
+    derivatives_in_prolate_spheroidal_coordinates()
+    #derivatives_in_oblate_spheroidal_coordinates()
+    #derivatives_in_bipolar_coordinates()
+    #derivatives_in_toroidal_coordinates()
 
     # xpdf()
     xpdf(pdfprog=None)

galgebra/_utils/simplify.py

@@ -0,0 +1,51 @@
+"""Compatibility helpers for simplification across SymPy releases."""
+
+import re
+
+import sympy
+from sympy import preorder_traversal, simplify, trigsimp
+from sympy.functions.elementary.hyperbolic import HyperbolicFunction
+from sympy.functions.elementary.trigonometric import TrigonometricFunction
+
+
+def _major_minor(version):
+    """Return the leading major and minor numbers from a version string."""
+    match = re.match(r'^(\d+)\.(\d+)', version)
+    if match is None:
+        return (0, 0)
+    return tuple(map(int, match.groups()))
+
+
+_SYMPY_MAJOR_MINOR = _major_minor(sympy.__version__)
+
+# SymPy 1.13's gh-26390 added a nested replace traversal to the FU
+# simplifier. Its cost is proportional to the expression tree size times the
+# number of trig and hyperbolic nodes. Values below this limit are cheap enough
+# to retain SymPy's canonical ``simplify`` output.
+_FU_TRAVERSAL_COST_LIMIT = 4096
+_TRIG_FUNCTIONS = (TrigonometricFunction, HyperbolicFunction)
+
+
+def _has_expensive_fu_traversal(expr):
+    """Whether ``simplify`` is likely to hit SymPy's slow FU traversal."""
+    if _SYMPY_MAJOR_MINOR < (1, 13):
+        return False
+
+    trig_nodes = 0
+    for node_count, node in enumerate(preorder_traversal(expr), 1):
+        if isinstance(node, _TRIG_FUNCTIONS):
+            trig_nodes += 1
+        if node_count * trig_nodes >= _FU_TRAVERSAL_COST_LIMIT:
+            return True
+    return False
+
+
+def simplify_compat(expr):
+    """Simplify while avoiding a known SymPy 1.13+ performance regression.
+
+    Remove this fallback after SymPy replaces the nested traversal introduced
+    by gh-26390 and galgebra's minimum supported SymPy includes that fix.
+    """
+    if _has_expensive_fu_traversal(expr):
+        return trigsimp(expr, method='old')
+    return simplify(expr)

galgebra/metric.py

@@ -14,6 +14,7 @@
 
 from . import printer
 from ._utils import cached_property as _cached_property
+from ._utils.simplify import simplify_compat
 from .atoms import (
     BasisVectorSymbol, DotProductSymbol, MatrixFunction, Determinant,
 )
@@ -299,7 +300,7 @@ def symbols_list(s, indices=None, sub=True, commutative=False):
 
 
 class Simp:
-    modes = [simplify]
+    modes = [simplify_compat]
 
     @staticmethod
     def profile(s):

scripts/validate_nb_refresh.py

@@ -49,7 +49,7 @@
    DeprecationWarnings from mpmath are environment-specific and ignored.
    Only ``display_data`` and ``execute_result`` outputs are compared.
 
-5. SymPy 1.13 ``trigsimp(method='old')`` algebraic form differences
+5. SymPy 1.13 compatibility-fallback algebraic form differences
    (curvilinear coordinates example, ``examples/ipython/LaTeX.ipynb``):
 
    a. Pythagorean identity: ``sin²(η)+sinh²(ξ)`` <-> ``-cos²(η)+cosh²(ξ)``
@@ -58,14 +58,13 @@
       form <-> distributed ``A + B/r + C/r²`` form.
    d. Whitespace inside ``\\frac{...}`` arguments.
    e. Outer ``\\left(…\\right)`` wrapper before a basis blade.
+   f. Spherical curl: the former notebook-wide ``trigsimp(method='old')``
+      profile wrote the radial coefficient as one fraction and kept a minus
+      sign outside the polar coefficient.  Scoped routing restores the
+      standard ``simplify`` factored/sign-distributed forms.
 
    **Known remaining differences (require symbolic algebra to verify):**
 
-   * Spherical curl ``e_r`` component: SymPy writes a parenthesised
-     ``(A/tan + B - C/sin²)|sin|`` form; trigsimp produces a single
-     fraction over ``tan|sin|``.  Both are mathematically equal but have
-     fundamentally different structure.
-
    * Prolate-spheroidal divergence: SymPy collects all terms into one
      large fraction; trigsimp distributes into seven separate fractions.
      The two forms are mathematically equal and can be verified with
@@ -121,6 +120,55 @@ def _norm_sum_sqrt_to_power32(text: str) -> str:
     return re.sub(pattern, replacement, text)
 
 
+def _norm_spherical_curl(text: str) -> str:
+    r"""Normalize the two known forms of the spherical-curl coefficients.
+
+    The spherical coordinates in this notebook are real, so
+    ``sin(theta)**2 == Abs(sin(theta))**2``.  Together with
+    ``sin(theta)*cos(theta)*tan(theta) == sin(theta)**2``, that converts the
+    old single fraction into the factored form below.  The exact replacement
+    deliberately does not accept nearby expressions: an unverified change
+    must still fail validation.
+    """
+    old_radial = (
+        r'\frac{A^{\phi }  {\sin{\left (\theta  \right )}}^{2} + '
+        r'A^{\phi }  \sin{\left (\theta  \right )} '
+        r'\cos{\left (\theta  \right )} \tan{\left (\theta  \right )} + '
+        r'{\sin{\left (\theta  \right )}}^{2} '
+        r'\tan{\left (\theta  \right )} \partial_{\theta } A^{\phi }  - '
+        r'\tan{\left (\theta  \right )} \partial_{\phi } A^{\theta } }'
+        r'{\tan{\left (\theta  \right )} '
+        r'\left|{\sin{\left (\theta  \right )}}\right|} '
+        r'\boldsymbol{e}_{r}'
+    )
+    factored_radial = (
+        r'\left(\frac{2 A^{\phi } }{\tan{\left (\theta  \right )}} + '
+        r'\partial_{\theta } A^{\phi }  - '
+        r'\frac{\partial_{\phi } A^{\theta } }'
+        r'{{\sin{\left (\theta  \right )}}^{2}}\right) '
+        r'\left|{\sin{\left (\theta  \right )}}\right| '
+        r'\boldsymbol{e}_{r}'
+    )
+    old_polar = (
+        r'- \frac{r^{2} {\sin{\left (\theta  \right )}}^{2} '
+        r'\partial_{r} A^{\phi }  + 2 r A^{\phi }  '
+        r'{\sin{\left (\theta  \right )}}^{2} - '
+        r'\partial_{\phi } A^{r} }'
+        r'{r^{2} \left|{\sin{\left (\theta  \right )}}\right|} '
+        r'\boldsymbol{e}_{\theta }'
+    )
+    sign_distributed_polar = (
+        r'+ \frac{- r^{2} {\sin{\left (\theta  \right )}}^{2} '
+        r'\partial_{r} A^{\phi }  - 2 r A^{\phi }  '
+        r'{\sin{\left (\theta  \right )}}^{2} + '
+        r'\partial_{\phi } A^{r} }'
+        r'{r^{2} \left|{\sin{\left (\theta  \right )}}\right|} '
+        r'\boldsymbol{e}_{\theta }'
+    )
+    text = text.replace(old_radial, factored_radial)
+    return text.replace(old_polar, sign_distributed_polar)
+
+
 # ---------------------------------------------------------------------------
 # Brace-counting helpers for structural LaTeX normalizers
 # ---------------------------------------------------------------------------
@@ -293,10 +341,11 @@ def _norm_collapse_spaces(text: str) -> str:
 LATEX_NORMALIZERS = [
     _norm_array_colspec,
     _norm_cdot,
-    # SymPy 1.13 (trigsimp method='old') uses different algebraic forms for
-    # some curvilinear-coordinate expressions.  The normalizers below bring
-    # both forms to a common representation so the validator can confirm the
-    # changes are cosmetic.
+    # The SymPy 1.13 compatibility fallback uses different algebraic forms
+    # for some curvilinear-coordinate expressions. The normalizers below
+    # bring both forms to a common representation so the validator can
+    # confirm the changes are cosmetic.
+    _norm_spherical_curl,           # factored/sign-distributed spherical-curl coefficients
     _norm_sin2_sinh2_identity,      # sin²+sinh²  ↔  -cos²+cosh²  (prolate spheroidal)
     _norm_sum_sqrt_to_power32,      # (X²+Y²)^{3/2}  ↔  X²√(…)+Y²√(…)  (paraboloidal)
     _norm_distribute_r2_denominator,  # \frac{r²A+rB+C}{r²}  ↔  A+B/r+C/r²  (spherical)

test/test_simplify.py

@@ -0,0 +1,92 @@
+from unittest import mock
+
+from sympy import Add, cos, cosh, sin, sinh, symbols
+
+from galgebra._utils import simplify as simplify_module
+from galgebra.ga import Ga
+from galgebra.metric import Simp
+
+
+x = symbols('x')
+
+
+def test_major_minor():
+    assert simplify_module._major_minor('1.13.3') == (1, 13)
+    assert simplify_module._major_minor('1.15.dev') == (1, 15)
+    assert simplify_module._major_minor('unknown') == (0, 0)
+
+
+def test_small_expression_uses_standard_simplify():
+    expr = sin(x)**2 + cos(x)**2
+
+    with (
+        mock.patch.object(simplify_module, 'simplify', return_value=1) as new,
+        mock.patch.object(simplify_module, 'trigsimp') as old,
+    ):
+        assert simplify_module.simplify_compat(expr) == 1
+
+    new.assert_called_once_with(expr)
+    old.assert_not_called()
+
+
+def test_large_mixed_expression_avoids_standard_simplify():
+    terms = [sin(x + i) + sinh(x + i) for i in range(40)]
+    expr = Add(*terms)
+
+    with (
+        mock.patch.object(simplify_module, 'simplify') as new,
+        mock.patch.object(simplify_module, 'trigsimp', return_value=expr) as old,
+    ):
+        assert simplify_module.simplify_compat(expr) == expr
+
+    new.assert_not_called()
+    old.assert_called_once_with(expr, method='old')
+
+
+def test_sympy_before_1_13_uses_standard_simplify():
+    terms = [sin(x + i) + sinh(x + i) for i in range(40)]
+    expr = Add(*terms)
+
+    with (
+        mock.patch.object(simplify_module, '_SYMPY_MAJOR_MINOR', (1, 12)),
+        mock.patch.object(simplify_module, 'simplify', return_value=1) as new,
+        mock.patch.object(simplify_module, 'trigsimp') as old,
+    ):
+        assert simplify_module.simplify_compat(expr) == 1
+
+    new.assert_called_once_with(expr)
+    old.assert_not_called()
+
+
+def test_custom_simp_profile_overrides_compatibility_default():
+    original_modes = Simp.modes[:]
+    custom = mock.Mock(return_value=x)
+    Simp.profile([custom])
+    try:
+        assert Simp.apply(sin(x)) == x
+    finally:
+        Simp.profile(original_modes)
+
+    custom.assert_called_once_with(sin(x))
+
+
+def test_prolate_spheroidal_divergence_renders():
+    a = symbols('a', real=True)
+    coords = xi, eta, phi = symbols('xi eta phi', real=True)
+    ps3d, *_ = Ga.build(
+        'e_xi e_eta e_phi',
+        X=[
+            a*sinh(xi)*sin(eta)*cos(phi),
+            a*sinh(xi)*sin(eta)*sin(phi),
+            a*cosh(xi)*cos(eta),
+        ],
+        coords=coords,
+        norm=True,
+    )
+    vector = ps3d.mv('A', 'vector', f=True)
+
+    rendered = str(ps3d.grad | vector)
+
+    assert 'D{eta}A__eta' in rendered
+    assert 'D{phi}A__phi' in rendered
+    assert 'D{xi}A__xi' in rendered