Do not simplify in overloaded operators, type mappers

.github/workflows/ci.yml

@@ -10,8 +10,8 @@ on:
         - cron:  '17 3 * * 0'
 
 concurrency:
-  group: ${{ github.head_ref || github.ref_name }}
-  cancel-in-progress: true
+    group: ${{ github.head_ref || github.ref_name }}
+    cancel-in-progress: true
 
 jobs:
     typos:
@@ -77,7 +77,7 @@ jobs:
         runs-on: ubuntu-latest
         strategy:
             matrix:
-                python-version: ["3.8", "3.9", "3.x"]
+                python-version: ["3.10", "3.12", "3.x"]
         steps:
         -   uses: actions/checkout@v4
         -

doc/algorithms.rst

@@ -2,12 +2,3 @@ Algorithms
 ==========
 
 .. automodule:: pymbolic.algorithm
-
-.. autofunction:: integer_power
-.. autofunction:: extended_euclidean
-.. autofunction:: gcd
-.. autofunction:: lcm
-.. autofunction:: fft
-.. autofunction:: ifft
-.. autofunction:: sym_fft
-

doc/conf.py

@@ -31,3 +31,8 @@
     "ExpressionT": "ExpressionT",
     "ArithmeticExpressionT": "ArithmeticExpressionT",
 }
+
+import sys
+
+
+sys._BUILDING_SPHINX_DOCS = True

doc/mappers.rst

@@ -44,5 +44,10 @@ Analysis tools
 
 .. automodule:: pymbolic.mapper.analysis
 
+Simplification
+^^^^^^^^^^^^^^
+
+.. automodule:: pymbolic.mapper.flattener
+
 
 .. vim: sw=4

pymbolic/__init__.py

@@ -38,8 +38,6 @@
 from .mapper import flattener
 from . import primitives
 
-from .polynomial import Polynomial
-
 from .primitives import (Variable as var,  # noqa: N813
     Variable,
     Expression,
@@ -73,7 +71,6 @@
     "Expression",
     "ExpressionT",
     "NumberT",
-    "Polynomial",
     "ScalarT",
     "Variable",
     "compile",

pymbolic/algorithm.py

@@ -1,3 +1,15 @@
+"""
+.. autofunction:: integer_power
+.. autofunction:: extended_euclidean
+.. autofunction:: gcd
+.. autofunction:: lcm
+.. autofunction:: fft
+.. autofunction:: ifft
+.. autofunction:: sym_fft
+.. autofunction:: reduced_row_echelon_form
+.. autofunction:: solve_affine_equations_for
+"""
+
 from __future__ import annotations
 
 
@@ -23,7 +35,16 @@
 THE SOFTWARE.
 """
 
-from pytools import memoize
+import operator
+import sys
+from typing import TYPE_CHECKING, overload
+from warnings import warn
+
+from pytools import MovedFunctionDeprecationWrapper, memoize
+
+
+if TYPE_CHECKING or getattr(sys, "_BUILDING_SPHINX_DOCS", None):
+    import numpy as np
 
 
 # {{{ integer powers
@@ -281,10 +302,47 @@ def csr_matrix_multiply(S, x):  # noqa
     return result
 
 
-# {{{ gaussian elimination
+# {{{ reduced_row_echelon_form
 
-def gaussian_elimination(mat, rhs):
+@overload
+def reduced_row_echelon_form(
+            mat: np.ndarray,
+            *, integral: bool | None = None,
+        ) -> np.ndarray:
+    ...
+
+
+@overload
+def reduced_row_echelon_form(
+            mat: np.ndarray,
+            rhs: np.ndarray,
+            *, integral: bool | None = None,
+        ) -> tuple[np.ndarray, np.ndarray]:
+    ...
+
+
+def reduced_row_echelon_form(
+            mat: np.ndarray,
+            rhs: np.ndarray | None = None,
+            integral: bool | None = None,
+        ) -> tuple[np.ndarray, np.ndarray] | np.ndarray:
     m, n = mat.shape
+
+    mat = mat.copy()
+    if rhs is not None:
+        rhs = rhs.copy()
+
+    if integral is None:
+        warn(
+             "Not specifying 'integral' is deprecated, please add it as an argument. "
+             "This will stop being supported in 2025.",
+             DeprecationWarning, stacklevel=2)
+
+    if integral:
+        div_func = operator.floordiv
+    else:
+        div_func = operator.truediv
+
     i = 0
     j = 0
 
@@ -303,8 +361,9 @@ def gaussian_elimination(mat, rhs):
             # swap rows i and nonz
             mat[i], mat[nonz_row] = \
                     (mat[nonz_row].copy(), mat[i].copy())
-            rhs[i], rhs[nonz_row] = \
-                    (rhs[nonz_row].copy(), rhs[i].copy())
+            if rhs is not None:
+                rhs[i], rhs[nonz_row] = \
+                        (rhs[nonz_row].copy(), rhs[i].copy())
 
             for u in range(0, m):
                 if u == i:
@@ -314,32 +373,51 @@ def gaussian_elimination(mat, rhs):
                     continue
 
                 ell = lcm(mat[u, j], mat[i, j])
-                u_fac = ell//mat[u, j]
-                i_fac = ell//mat[i, j]
+                u_fac = div_func(ell, mat[u, j])
+                i_fac = div_func(ell, mat[i, j])
 
                 mat[u] = u_fac*mat[u] - i_fac*mat[i]
-                rhs[u] = u_fac*rhs[u] - i_fac*rhs[i]
+                if rhs is not None:
+                    rhs[u] = u_fac*rhs[u] - i_fac*rhs[i]
 
                 assert mat[u, j] == 0
 
             i += 1
 
         j += 1
 
-    for i in range(m):
-        g = gcd_many(*(
-            [a for a in mat[i] if a]
-            +
-            [a for a in rhs[i] if a]))
+    if integral:
+        for i in range(m):
+            g = gcd_many(*(
+                [a for a in mat[i] if a]
+                +
+                [a for a in rhs[i] if a] if rhs is not None else []))
+
+            mat[i] //= g
+            if rhs is not None:
+                rhs[i] //= g
+
+    import numpy as np
+
+    from pymbolic.mapper.flattener import flatten
+    vec_flatten = np.vectorize(flatten, otypes=[object])
 
-        mat[i] //= g
-        rhs[i] //= g
+    for i in range(m):
+        mat[i] = vec_flatten(mat[i])
+        if rhs is not None:
+            rhs[i] = vec_flatten(rhs[i])
 
-    return mat, rhs
+    if rhs is None:
+        return mat
+    else:
+        return mat, rhs
 
 # }}}
 
 
+gaussian_elimination = MovedFunctionDeprecationWrapper(reduced_row_echelon_form, "2025")
+
+
 # {{{ symbolic (linear) equation solving
 
 def solve_affine_equations_for(unknowns, equations):
@@ -393,7 +471,7 @@ def solve_affine_equations_for(unknowns, equations):
 
     # }}}
 
-    mat, rhs_mat = gaussian_elimination(mat, rhs_mat)
+    mat, rhs_mat = reduced_row_echelon_form(mat, rhs_mat, integral=True)
 
     # FIXME /!\ Does not check for overdetermined system.
 
@@ -411,7 +489,8 @@ def solve_affine_equations_for(unknowns, equations):
         div = mat[nonz_row, j]
 
         unknown_val = int(rhs_mat[nonz_row, -1]) // div
-        for parameter, coeff in zip(parameters_list, rhs_mat[nonz_row]):
+        for parameter, coeff in zip(
+                    parameters_list, rhs_mat[nonz_row, :-1], strict=True):
             unknown_val += (int(coeff) // div) * parameter
 
         result[unknown] = unknown_val

pymbolic/compiler.py

@@ -26,7 +26,7 @@
 import math
 
 import pymbolic
-from pymbolic.mapper.stringifier import PREC_NONE, PREC_POWER, PREC_SUM, StringifyMapper
+from pymbolic.mapper.stringifier import PREC_NONE, StringifyMapper
 
 
 class CompileMapper(StringifyMapper):
@@ -45,34 +45,6 @@ def map_constant(self, expr, enclosing_prec):
 
         return repr(expr)
 
-    def map_polynomial(self, expr, enclosing_prec):
-        # Use Horner's scheme to evaluate the polynomial
-
-        sbase = self(expr.base, PREC_POWER)
-
-        def stringify_exp(exp):
-            if exp == 0:
-                return ""
-            elif exp == 1:
-                return f"*{sbase}"
-            else:
-                return f"*{sbase}**{exp}"
-
-        result = ""
-        rev_data = expr.data[::-1]
-        for i, (exp, coeff) in enumerate(rev_data):
-            if i+1 < len(rev_data):
-                next_exp = rev_data[i+1][0]
-            else:
-                next_exp = 0
-            result = "({}+{}){}".format(result, self(coeff, PREC_SUM),
-                    stringify_exp(exp-next_exp))
-
-        if enclosing_prec > PREC_SUM and len(expr.data) > 1:
-            return f"({result})"
-        else:
-            return result
-
     def map_numpy_array(self, expr, enclosing_prec):
         def stringify_leading_dimension(ary):
             if len(ary.shape) == 1: