Implement divmod on Expr

Closes sympy#7889.

Note there is a separate issue with a discrepancy between the builtin divmod
and divmod on SymPy number types that is not fixed here (sympy#15561).

sympy/core/expr.py

@@ -217,6 +217,19 @@ def __rfloordiv__(self, other):
         from sympy.functions.elementary.integers import floor
         return floor(other / self)
 
+
+    @_sympifyit('other', NotImplemented)
+    @call_highest_priority('__rdivmod__')
+    def __divmod__(self, other):
+        from sympy.functions.elementary.integers import floor
+        return floor(self / other), Mod(self, other)
+
+    @_sympifyit('other', NotImplemented)
+    @call_highest_priority('__divmod__')
+    def __rdivmod__(self, other):
+        from sympy.functions.elementary.integers import floor
+        return floor(other / self), Mod(other, self)
+
     def __int__(self):
         # Although we only need to round to the units position, we'll
         # get one more digit so the extra testing below can be avoided

sympy/core/tests/test_arit.py

@@ -1985,3 +1985,8 @@ def test_Add_is_zero():
 
 def test_issue_14392():
     assert (sin(zoo)**2).as_real_imag() == (nan, nan)
+
+def test_divmod():
+    assert divmod(x, y) == (x//y, x % y)
+    assert divmod(x, 3) == (x//3, x % 3)
+    assert divmod(3, x) == (3//x, 3 % x)