Repository: sympy/sympy
Pull request: Introduce general add, mul, and pow function
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Closes #18769
Fixes #18768

#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- core
    - Extensible `add`, `mul` and `power` functions are introduced to allow sympy objects to define what classes should be used for them in place of `Add`, `Mul` and `Pow` (e.g. matrices use `MatAdd`). This is an experimental approach aimed at enabling the behaviour of core routines (expand, collect, etc) to be customised by user-defined types (e.g. `MatAdd` rather than `Add`). This mechanism is still experimental, is not fully implemented across the core and might be changed or removed in a future release of sympy.
<!-- END RELEASE NOTES -->

Problem statement:
Connecting Add, Mul to MatAdd, MatMul
Current matrix module uses `Basic._exec_constructor_postprocessor` method to convert `Add` and `Mul` containing `MatrixExpr` to `MatAdd` and `MatMul`, respectly.  
However, this behavior seems to be somewhat unefficient...

1. It does not pass any option to `MatAdd` and `MatMul`.
    ```python
    >>> from sympy import Add, MatAdd, MatrixSymbol
    >>> A,B = MatrixSymbol('A', 2,2), MatrixSymbol('B', 3,3)
    >>> MatAdd(A,B, check=True)    # Raises error
    Traceback (most recent call last)
    ...
    ShapeError: Matrices A and B are not aligned
    >>> Add(A,B, check=True)    # Does not raise error
    A + B
    ```

2. It takes every steps  (`flatten` and `_from_args` methods of `Add` and `Mul`) to create `Add` or `Mul` object, then takes every steps of `MatAdd` or `MatMul` again.
Current approach behaves in the way that, if new classes with their own postprocessors are added in the future, every postprocessors will be executed one by one.

Suggestion: make the arguments 'preprocessed' rather than 'postprocessed'. This can be done by comparing `_op_priority` of the arguments.

Hints from issue discussion:

Commit message: Apply suggestions from code review

Co-authored-by: Eric Wieser <wieser.eric@gmail.com>

Patch to review:
diff --git a/sympy/core/add.py b/sympy/core/add.py
index a37140b836cc..651f6b627caa 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -5,7 +5,7 @@
 from .parameters import global_parameters
 from .logic import _fuzzy_group, fuzzy_or, fuzzy_not
 from .singleton import S
-from .operations import AssocOp
+from .operations import AssocOp, AssocOpDispatcher
 from .cache import cacheit
 from .numbers import ilcm, igcd
 from .expr import Expr
@@ -66,7 +66,6 @@ def _unevaluated_Add(*args):
         newargs.insert(0, co)
     return Add._from_args(newargs)
 
-
 class Add(Expr, AssocOp):
 
     __slots__ = ()
@@ -1127,6 +1126,7 @@ def __neg__(self):
             return super().__neg__()
         return Add(*[-i for i in self.args])
 
+add = AssocOpDispatcher('add')
 
 from .mul import Mul, _keep_coeff, prod
 from sympy.core.numbers import Rational
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 1e5d6ca9e357..1337d116fd8b 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -163,6 +163,18 @@ def __eq__(self, other):
     # something better and more powerful.  See issue 5510.
     _op_priority = 10.0
 
+    # Handler methods for dispatching
+    # See operations.AssocOpDispatcher
+    @staticmethod
+    def _add_handler(*args, **kwargs):
+        from .add import Add
+        return Add(*args, **kwargs)
+
+    @staticmethod
+    def _mul_handler(*args, **kwargs):
+        from .mul import Mul
+        return Mul(*args, **kwargs)
+
     def __pos__(self):
         return self
 
diff --git a/sympy/core/mul.py b/sympy/core/mul.py
index d4af08fa7390..95969074dae1 100644
--- a/sympy/core/mul.py
+++ b/sympy/core/mul.py
@@ -5,7 +5,7 @@
 from .sympify import sympify
 from .basic import Basic
 from .singleton import S
-from .operations import AssocOp
+from .operations import AssocOp, AssocOpDispatcher
 from .cache import cacheit
 from .logic import fuzzy_not, _fuzzy_group, fuzzy_and
 from .compatibility import reduce
@@ -1901,6 +1901,7 @@ def as_ordered_factors(self, order=None):
     def _sorted_args(self):
         return tuple(self.as_ordered_factors())
 
+mul = AssocOpDispatcher('mul')
 
 def prod(a, start=1):
     """Return product of elements of a. Start with int 1 so if only
diff --git a/sympy/core/operations.py b/sympy/core/operations.py
index d36a32e538ae..ca33d6d1eb7c 100644
--- a/sympy/core/operations.py
+++ b/sympy/core/operations.py
@@ -1,3 +1,5 @@
+import itertools
+from operator import attrgetter
 from typing import Tuple
 
 from sympy.utilities.exceptions import SymPyDeprecationWarning
@@ -8,7 +10,9 @@
 from sympy.core.compatibility import ordered
 from sympy.core.logic import fuzzy_and
 from sympy.core.parameters import global_parameters
-from sympy.utilities.iterables import sift
+from sympy.utilities.iterables import sift, tied_max
+from sympy.multipledispatch.dispatcher import (Dispatcher, ambiguity_register_error,
+    str_signature)
 
 
 class AssocOp(Basic):
@@ -21,6 +25,20 @@ class AssocOp(Basic):
 
     This is an abstract base class, concrete derived classes must define
     the attribute `identity`.
+
+    Parameters
+    ==========
+
+    *args :
+        Arguments which are operated
+
+    sympify : bool, optional
+        Default is ``True``. ``False`` may be passed only when all arguments are
+            already sympified, for faster processing.
+
+    evaluate : bool, optional
+        Default is ``True``. If ``False``, do not evaluate the operation.
+
     """
 
     # for performance reason, we don't let is_commutative go to assumptions,
@@ -30,9 +48,13 @@ class AssocOp(Basic):
     _args_type = None
 
     @cacheit
-    def __new__(cls, *args, **options):
+    def __new__(cls, *args, sympify=True, **options):
         from sympy import Order
-        args = list(map(_sympify, args))
+
+        # Allow faster processing by passing ``sympify=False``, if all arguments
+        # are already sympified.
+        if sympify:
+            args = list(map(_sympify, args))
 
         # Disallow non-Expr args in Add/Mul
         typ = cls._args_type
@@ -504,3 +526,188 @@ def args(self):
     @staticmethod
     def _compare_pretty(a, b):
         return (str(a) > str(b)) - (str(a) < str(b))
+
+class AssocOpDispatcher:
+    """
+    Handler dispatcher for associative operators
+
+    .. notes::
+       This approach is experimental, and can be replaced or deleted in the future.
+       See https://github.com/sympy/sympy/pull/19463.
+
+    Explanation
+    ===========
+
+    If arguments of different kind are passed, their ``_op_priority`` are compared to
+    choose the type with the highest priority. If this is ambiguous, priority relations
+    which are registered by ``register_priority`` method are referred to.
+    Once a type is selected, handler method of the type, which is a static method, is used
+    to perform the operation.
+
+    Priority registration is unordered. You cannot make ``A*B`` and ``B*A`` use different
+    handler method. All logic dealing with the order of arguments must be implemented in
+    the handler method.
+
+    Examples
+    ========
+
+    >>> from sympy import Add, Expr, Symbol
+    >>> from sympy.core.add import add
+
+    >>> class NewExpr(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         return NewAdd(*args, **kwargs)
+    >>> class NewAdd(NewExpr, Add):
+    ...     pass
+
+    >>> @add.register_priority(Expr, NewExpr)
+    ... @add.register_priority(NewExpr, NewExpr)
+    ... def _(_1, _2):
+    ...     return NewExpr
+
+    >>> a, b = Symbol('a'), NewExpr()
+    >>> add(a, b) == NewAdd(a, b)
+    True
+
+    Registered priority selector can return any type.
+
+    >>> class Foo(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         raise NotImplementedError
+    >>> class Bar(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         raise NotImplementedError
+    >>> class FooBar(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         return NewAdd(*args, **kwargs)
+
+    >>> @add.register_priority(Foo, Bar)
+    ... def _(_1, _2):
+    ...     return FooBar
+
+    >>> add(Foo(), Bar()) == NewAdd(Foo(), Bar())
+    True
+
+    """
+    def __init__(self, name, doc=None):
+        self.name = name
+        self.doc = doc
+        self.handlermethod = "_%s_handler" % name
+        self._handlergetter = attrgetter(self.handlermethod)
+        self._dispatcher = Dispatcher(name)
+
+    def __repr__(self):
+        return "<dispatched %s>" % self.name
+
+    def register_priority(self, cls1, cls2, on_ambiguity=ambiguity_register_error):
+        """
+        Register the priority between two types, in both straight and reversed order.
+        Registered function must return the class which has the handler that is to be used.
+        """
+        def _(func):
+            self._dispatcher.add((cls1, cls2), func, on_ambiguity=on_ambiguity)
+            self._dispatcher.add((cls2, cls1), func, on_ambiguity=on_ambiguity)
+            return func
+        return _
+
+    @cacheit
+    def __call__(self, *args, sympify=True, **kwargs):
+        """
+        Parameters
+        ==========
+
+        *args :
+            Arguments which are operated
+
+        sympify : bool, optional
+            Default is ``True``. ``False`` may be passed only when all arguments are
+            already sympified, for faster processing.
+        """
+        if sympify:
+            args = tuple(map(_sympify, args))
+        types = frozenset(map(type, args))
+
+        # If sympify=True was passed, no need to sympify again so pass sympify=False.
+        # If sympify=False was passed, all args are already sympified so pass sympify=False.
+        return self.dispatch(types)(*args, sympify=False, **kwargs)
+
+    @cacheit
+    def dispatch(self, types):
+        """
+        Select the type with highest priority, and return its handler.
+        """
+
+        # Quick exit for the case where all arguments are of same type
+        # or all handlermethods are same
+        handlers = set(map(self._handlergetter, types))
+        if len(handlers) == 1:
+            h, = handlers
+            return h
+
+        # Sieve the types with low op_priority or no handler
+        # This makes it faster than using dispatcher on every types.
+        sieved_data = tied_max(
+            types, key=lambda typ: getattr(typ, '_op_priority', 0)
+        )
+
+        # Only one type remains
+        if len(sieved_data) == 1:
+            typ, = sieved_data
+            return self._handlergetter(typ)
+
+        # Recursively select with registered binary priority
+        for i, typ in enumerate(sieved_data):
+            if i == 0:
+                result_type = typ
+            else:
+                dispatched_selector = self._dispatcher.dispatch(result_type, typ)
+                selected_type = dispatched_selector(result_type, typ)
+
+                if not isinstance(selected_type, type):
+                    raise RuntimeError(
+                        "Dispatcher for {!r} and {!r} must return a type, but got {!r}".format(
+                        result_type, typ, selected_type
+                    ))
+
+                handler = self._handlergetter(selected_type)
+                if handler is None:
+                    raise RuntimeError(
+                        "Dispatcher for {!r} and {!r} returned {!r} which does not have {!r} method.".format(
+                        result_type, typ, selected_type, self.handlermethod
+                    ))
+
+                result_type = selected_type
+
+        # return handler method of selected argument
+        return handler
+
+    @property
+    def __doc__(self):
+        docs = [
+            "Multiply dispatched associative operator: %s" % self.name,
+            "Note that support for this is experimental, see the docs for :class:`AssocOpDispatcher` for details"
+        ]
+
+        if self.doc:
+            docs.append(self.doc)
+
+        docs.append("Registered priority relations are as follows:")
+
+        other = []
+        for func, sigs in itertools.groupby(self._dispatcher.ordering[::-1], lambda sig: self._dispatcher.funcs[sig]):
+            if func.__doc__:
+                s = 'Inputs: %s\n' % ', '.join('<%s>' % str_signature(sig) for sig in sigs)
+                s += '-' * len(s) + '\n'
+                s += func.__doc__.strip()
+                docs.append(s)
+            else:
+                other += list(sigs)
+
+        if other:
+            docs.append('Other signatures:\n    ' + '\n    '.join(other))
+
+        return '\n\n'.join(docs)
diff --git a/sympy/core/power.py b/sympy/core/power.py
index 0a3b8b564c42..6a610be5ad60 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -12,6 +12,7 @@
 from .parameters import global_parameters
 from sympy.utilities.iterables import sift
 from sympy.utilities.exceptions import SymPyDeprecationWarning
+from sympy.multipledispatch import Dispatcher
 
 from mpmath.libmp import sqrtrem as mpmath_sqrtrem
 
@@ -1732,6 +1733,8 @@ def _eval_difference_delta(self, n, step):
             new_e = e.subs(n, n + step)
             return (b**(new_e - e) - 1) * self
 
+power = Dispatcher('power')
+power.add((object, object), Pow)
 
 from .add import Add
 from .numbers import Integer
diff --git a/sympy/core/tests/test_operations.py b/sympy/core/tests/test_operations.py
index 7555398f4395..a2424b3afd61 100644
--- a/sympy/core/tests/test_operations.py
+++ b/sympy/core/tests/test_operations.py
@@ -1,8 +1,9 @@
-from sympy import Integer, S, symbols, Mul
+from sympy import Integer, S, Symbol, symbols, Expr
 from sympy.core.operations import AssocOp, LatticeOp
 from sympy.testing.pytest import raises
 from sympy.core.sympify import SympifyError
-from sympy.core.add import Add
+from sympy.core.add import Add, add
+from sympy.core.mul import Mul, mul
 
 # create the simplest possible Lattice class
 
@@ -64,3 +65,69 @@ class MyAssoc(AssocOp):
     assert v.args == (u, d)
     # like Add, any unevaluated outer call will flatten inner args
     assert MyAssoc(a, v).args == (a, b, c, d)
+
+def test_add_dispatcher():
+
+    class NewBase1(Expr):
+        @staticmethod
+        def _add_handler(*args, **kwargs):
+            return NewAdd1(*args, **kwargs)
+    class NewAdd1(NewBase1, Add):
+        pass
+    class NewBase2(Expr):
+        _op_priority = 11
+        @staticmethod
+        def _add_handler(*args, **kwargs):
+            return NewAdd2(*args, **kwargs)
+    class NewAdd2(NewBase2, Add):
+        pass
+    a, b, c = Symbol('a'), NewBase1(), NewBase2()
+
+    @add.register_priority(Expr, NewBase1)
+    @add.register_priority(NewBase1, NewBase1)
+    def _(_1, _2):
+        return NewBase1
+
+    # Add called as fallback
+    assert add(1, 2) == Add(1, 2)
+    assert add(a, a) == Add(a, a)
+
+    # selection by registered priority
+    assert add(a,b,a) == NewAdd1(2*a, b)
+
+    # selection by _op_priority
+    assert add(a,c,a) == NewAdd2(2*a, c)
+    assert add(a,b,c) == NewAdd2(a,b,c)
+
+def test_mul_dispatcher():
+
+    class NewBase1(Expr):
+        @staticmethod
+        def _mul_handler(*args, **kwargs):
+            return NewMul1(*args, **kwargs)
+    class NewMul1(NewBase1, Mul):
+        pass
+    class NewBase2(Expr):
+        _op_priority = 11
+        @staticmethod
+        def _mul_handler(*args, **kwargs):
+            return NewMul2(*args, **kwargs)
+    class NewMul2(NewBase2, Mul):
+        pass
+    a, b, c = Symbol('a'), NewBase1(), NewBase2()
+
+    @mul.register_priority(Expr, NewBase1)
+    @mul.register_priority(NewBase1, NewBase1)
+    def _(_1, _2):
+        return NewBase1
+
+    # Mul called as fallback
+    assert mul(1, 2) == Mul(1, 2)
+    assert mul(a, a) == Mul(a, a)
+
+    # selection by registered priority
+    assert mul(a,b,a) == NewMul1(a**2, b)
+
+    # selection by _op_priority
+    assert mul(a,c,a) == NewMul2(a**2, c)
+    assert mul(a,b,c) == NewMul2(a,b,c)
diff --git a/sympy/core/tests/test_power.py b/sympy/core/tests/test_power.py
index 6aa6ef2e57cc..413fbd453217 100644
--- a/sympy/core/tests/test_power.py
+++ b/sympy/core/tests/test_power.py
@@ -12,7 +12,7 @@
 from sympy.series.order import O
 from sympy.sets import FiniteSet
 from sympy.core.expr import unchanged
-
+from sympy.core.power import power
 from sympy.testing.pytest import warns_deprecated_sympy
 
 
@@ -555,3 +555,27 @@ def test_issue_18762():
     e, p = symbols('e p')
     g0 = sqrt(1 + e**2 - 2*e*cos(p))
     assert len(g0.series(e, 1, 3).args) == 4
+
+def test_power_dispatcher():
+
+    class NewBase(Expr):
+        pass
+    class NewPow(NewBase, Pow):
+        pass
+    a, b = Symbol('a'), NewBase()
+
+    @power.register(Expr, NewBase)
+    @power.register(NewBase, Expr)
+    @power.register(NewBase, NewBase)
+    def _(a, b):
+        return NewPow(a, b)
+
+    # Pow called as fallback
+    assert power(2, 3) == 8*S.One
+    assert power(a, 2) == Pow(a, 2)
+    assert power(a, a) == Pow(a, a)
+
+    # NewPow called by dispatch
+    assert power(a, b) == NewPow(a, b)
+    assert power(b, a) == NewPow(b, a)
+    assert power(b, b) == NewPow(b, b)
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 83596decde7d..066b174a0a37 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -5,6 +5,8 @@
 import collections
 
 from sympy.core import S, Symbol, Tuple, Integer, Basic, Expr, Mul, Add
+from sympy.core.add import add
+from sympy.core.mul import mul
 from sympy.core.decorators import call_highest_priority
 from sympy.core.compatibility import SYMPY_INTS, default_sort_key
 from sympy.core.sympify import SympifyError, _sympify
@@ -593,6 +595,14 @@ def applyfunc(self, func):
         from .applyfunc import ElementwiseApplyFunction
         return ElementwiseApplyFunction(func, self)
 
+    @staticmethod
+    def _add_handler(*args, **kwargs):
+        return MatAdd(*args, **kwargs)
+
+    @staticmethod
+    def _mul_handler(*args, **kwargs):
+        return MatMul(*args, **kwargs)
+
 @dispatch(MatrixExpr, Expr)
 def _eval_is_eq(lhs, rhs): # noqa:F811
     return False
@@ -648,6 +658,13 @@ def _postprocessor(expr):
     "Add": [get_postprocessor(Add)],
 }
 
+@add.register_priority(Expr, MatrixExpr)
+@add.register_priority(MatrixExpr, MatrixExpr)
+@mul.register_priority(Expr, MatrixExpr)
+@mul.register_priority(MatrixExpr, MatrixExpr)
+def _(_1, _2):
+    "Use MatrixExpr's handler method."
+    return MatrixExpr
 
 def _matrix_derivative(expr, x):
     from sympy import Derivative
diff --git a/sympy/matrices/expressions/tests/test_matadd.py b/sympy/matrices/expressions/tests/test_matadd.py
index b4d0ec046106..d13890557cca 100644
--- a/sympy/matrices/expressions/tests/test_matadd.py
+++ b/sympy/matrices/expressions/tests/test_matadd.py
@@ -2,24 +2,27 @@
 from sympy.matrices.expressions.special import GenericZeroMatrix, ZeroMatrix
 from sympy.matrices import eye, ImmutableMatrix
 from sympy.core import Add, Basic, S
+from sympy.core.add import add
 from sympy.testing.pytest import XFAIL, raises
 
 X = MatrixSymbol('X', 2, 2)
 Y = MatrixSymbol('Y', 2, 2)
 
 def test_evaluate():
-    assert MatAdd(X, X, evaluate=True) == MatAdd(X, X).doit()
+    assert MatAdd(X, X, evaluate=True) == add(X, X, evaluate=True) == MatAdd(X, X).doit()
 
 def test_sort_key():
-    assert MatAdd(Y, X).doit().args == (X, Y)
+    assert MatAdd(Y, X).doit().args == add(Y, X).doit().args == (X, Y)
 
 
 def test_matadd_sympify():
     assert isinstance(MatAdd(eye(1), eye(1)).args[0], Basic)
+    assert isinstance(add(eye(1), eye(1)).args[0], Basic)
 
 
 def test_matadd_of_matrices():
     assert MatAdd(eye(2), 4*eye(2), eye(2)).doit() == ImmutableMatrix(6*eye(2))
+    assert add(eye(2), 4*eye(2), eye(2)).doit() == ImmutableMatrix(6*eye(2))
 
 
 def test_doit_args():
@@ -28,7 +31,8 @@ def test_doit_args():
     assert MatAdd(A, MatPow(B, 2)).doit() == A + B**2
     assert MatAdd(A, MatMul(A, B)).doit() == A + A*B
     assert (MatAdd(A, X, MatMul(A, B), Y, MatAdd(2*A, B)).doit() ==
-            MatAdd(3*A + A*B + B, X, Y))
+    add(A, X, MatMul(A, B), Y, add(2*A, B)).doit() ==
+    MatAdd(3*A + A*B + B, X, Y))
 
 
 def test_generic_identity():
@@ -40,5 +44,5 @@ def test_zero_matrix_add():
     assert Add(ZeroMatrix(2, 2), ZeroMatrix(2, 2)) == ZeroMatrix(2, 2)
 
 @XFAIL
-def test_matrix_add_with_scalar():
+def test_matrix_Add_with_scalar():
     raises(TypeError, lambda: Add(0, ZeroMatrix(2, 2)))
diff --git a/sympy/matrices/expressions/tests/test_matmul.py b/sympy/matrices/expressions/tests/test_matmul.py
index 340dd3456db0..529be8559f13 100644
--- a/sympy/matrices/expressions/tests/test_matmul.py
+++ b/sympy/matrices/expressions/tests/test_matmul.py
@@ -1,4 +1,5 @@
 from sympy.core import I, symbols, Basic, Mul, S
+from sympy.core.mul import mul
 from sympy.functions import adjoint, transpose
 from sympy.matrices import (Identity, Inverse, Matrix, MatrixSymbol, ZeroMatrix,
         eye, ImmutableMatrix)
@@ -153,6 +154,11 @@ def test_construction_with_Mul():
     assert Mul(C, D) == MatMul(C, D)
     assert Mul(D, C) == MatMul(D, C)
 
+def test_construction_with_mul():
+    assert mul(C, D) == MatMul(C, D)
+    assert mul(D, C) == MatMul(D, C)
+    assert mul(C, D) != MatMul(D, C)
+
 def test_generic_identity():
     assert MatMul.identity == GenericIdentity()
     assert MatMul.identity != S.One
diff --git a/sympy/multipledispatch/dispatcher.py b/sympy/multipledispatch/dispatcher.py
index 9c5f38651183..6cb30c1d4209 100644
--- a/sympy/multipledispatch/dispatcher.py
+++ b/sympy/multipledispatch/dispatcher.py
@@ -28,6 +28,34 @@ def ambiguity_warn(dispatcher, ambiguities):
     warn(warning_text(dispatcher.name, ambiguities), AmbiguityWarning)
 
 
+def ambiguity_register_error(dispatcher, ambiguities):
+    """
+    Automatically register function which raises error when ambiguous types
+    are registered.
+
+    Parameters
+    ----------
+    dispatcher : Dispatcher
+        The dispatcher on which the ambiguity was detected
+    ambiguities : set
+        Set of type signature pairs that are ambiguous within this dispatcher
+
+    See Also:
+        Dispatcher.add
+        ambiguity_warn
+    """
+    for amb in ambiguities:
+        @dispatcher.register(*super_signature(amb))
+        def _(*args, **kwargs):
+            """ Unimplemented handler """
+            types = tuple(type(a) for a in args)
+            raise NotImplementedError(
+                "Ambiguous signature for %s: <%s>" % (
+                    dispatcher.name, str_signature(types)
+                )
+            )
+
+
 _unresolved_dispatchers = set() # type: Set[Dispatcher]
 _resolve = [True]
 
diff --git a/sympy/multipledispatch/tests/test_dispatcher.py b/sympy/multipledispatch/tests/test_dispatcher.py
index b980a8b19afc..8edaae2eb6eb 100644
--- a/sympy/multipledispatch/tests/test_dispatcher.py
+++ b/sympy/multipledispatch/tests/test_dispatcher.py
@@ -1,6 +1,7 @@
 from sympy.multipledispatch.dispatcher import (Dispatcher, MDNotImplementedError,
                                          MethodDispatcher, halt_ordering,
-                                         restart_ordering)
+                                         restart_ordering,
+                                         ambiguity_register_error)
 from sympy.testing.pytest import raises, warns
 
 
@@ -262,3 +263,13 @@ def _(a):
         raise MDNotImplementedError()
 
     assert raises(NotImplementedError, lambda: f(1.0))
+
+def test_ambiguity_register_error():
+    f = Dispatcher('add')
+
+    # suppress warning for registering ambiguous signal
+    f.add((int, object), lambda x,y: x+y, ambiguity_register_error)
+    f.add((object, int), lambda x,y: x+y, ambiguity_register_error)
+
+    # raises error if ambiguous signal is passed
+    assert raises(NotImplementedError, lambda: f(1,2))
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index dd36dea304b5..da4258439c7d 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -819,6 +819,37 @@ def sift(seq, keyfunc, binary=False):
     return T, F
 
 
+def tied_max(iterable, *, key):
+    """
+    Return the list of items from *iterable* which returns
+    highest value when *key* is applied.
+
+    Examples
+    ========
+
+    >>> from sympy.utilities.iterables import tied_max
+    >>> iterable = ['abc', 'de', 'fgh', 'i']
+    >>> tied_max(iterable, key=len)
+    ['abc', 'fgh']
+
+    """
+    it = iter(iterable)
+    try:
+        first = next(it)
+    except StopIteration:
+        return []
+    max_key = key(first)
+    max_vals = [first]
+    for val in it:
+        k = key(val)
+        if k == max_key:
+            max_vals.append(val)
+        elif k > max_key:
+            max_key = k
+            max_vals = [val]
+    return max_vals
+
+
 def take(iter, n):
     """Return ``n`` items from ``iter`` iterator. """
     return [ value for _, value in zip(range(n), iter) ]


Code review comments:
diff_hunk: @@ -21,6 +25,20 @@ class AssocOp(Basic):
 
     This is an abstract base class, concrete derived classes must define
     the attribute `identity`.
+
+    Parameters
+    ==========
+
+    *args :
+        Arguments which are operated
+
+    sympify : bool, optional
path: sympy/core/operations.py
text: @user1:
Perhaps this should be made private (`_sympify` and not documented). This is good for core library code, but I don't think it would be a good idea for most users to use this.

diff_hunk: @@ -504,3 +526,188 @@ def args(self):
     @staticmethod
     def _compare_pretty(a, b):
         return (str(a) > str(b)) - (str(a) < str(b))
+
+class AssocOpDispatcher:
+    """
+    Handler dispatcher for associative operators
+
+    .. notes::
+       This approach is experimental, and can be replaced or deleted in the future.
+       See https://github.com/sympy/sympy/pull/19463.
+
+    Explanation
+    ===========
+
+    If arguments of different kind are passed, their ``_op_priority`` are compared to
+    choose the type with the highest priority. If this is ambiguous, priority relations
+    which are registered by ``register_priority`` method are referred to.
+    Once a type is selected, handler method of the type, which is a static method, is used
+    to perform the operation.
+
+    Priority registration is unordered. You cannot make ``A*B`` and ``B*A`` use different
+    handler method. All logic dealing with the order of arguments must be implemented in
+    the handler method.
+
+    Examples
+    ========
+
+    >>> from sympy import Add, Expr, Symbol
+    >>> from sympy.core.add import add
+
+    >>> class NewExpr(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         return NewAdd(*args, **kwargs)
+    >>> class NewAdd(NewExpr, Add):
+    ...     pass
+
+    >>> @add.register_priority(Expr, NewExpr)
+    ... @add.register_priority(NewExpr, NewExpr)
+    ... def _(_1, _2):
+    ...     return NewExpr
+
+    >>> a, b = Symbol('a'), NewExpr()
+    >>> add(a, b) == NewAdd(a, b)
+    True
+
+    Registered priority selector can return any type.
+
+    >>> class Foo(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         raise NotImplementedError
+    >>> class Bar(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         raise NotImplementedError
+    >>> class FooBar(Expr):
+    ...     @staticmethod
+    ...     def _add_handler(*args, **kwargs):
+    ...         return NewAdd(*args, **kwargs)
+
+    >>> @add.register_priority(Foo, Bar)
+    ... def _(_1, _2):
path: sympy/core/operations.py
text: @user1:
This relates to my other question, but I'm not a fan of using `_` for decorated function names. 

@author:
Maybe I can modify this to `register_priority(type1, type2, result_type)` while keeping `(type1, type2)` be dispatched.

@author:
Done.

@user2:
Using `_` as a decorated function name seems to be a blessed idiom for the cases where the function name is irrelevant because it's being used in e.g. single dispatch. I think flake8 special cases the name for this purpose.