Multi output primitives (#330)

* Multi-output primitives

---------

Co-authored-by: Angelos Katharopoulos <a_katharopoulos@apple.com>

python/src/ops.cpp

@@ -303,6 +303,33 @@ void init_ops(py::module_& m) {
         Returns:
             array: The quotient ``a / b``.
       )pbdoc");
+  m.def(
+      "divmod",
+      [](const ScalarOrArray& a_, const ScalarOrArray& b_, StreamOrDevice s) {
+        auto [a, b] = to_arrays(a_, b_);
+        return divmod(a, b, s);
+      },
+      "a"_a,
+      "b"_a,
+      py::pos_only(),
+      py::kw_only(),
+      "stream"_a = none,
+      R"pbdoc(
+        divmod(a: Union[scalar, array], b: Union[scalar, array], stream: Union[None, Stream, Device] = None) -> array
+
+        Element-wise quotient and remainder.
+
+        The fuction ``divmod(a, b)`` is equivalent to but faster than
+        ``(a // b, a % b)``. The function uses numpy-style broadcasting
+        semantics. Either or both input arrays can also be scalars.
+
+        Args:
+            a (array): Input array or scalar.
+            b (array): Input array or scalar.
+
+        Returns:
+            tuple(array, array): The quotient ``a // b`` and remainder ``a % b``.
+      )pbdoc");
   m.def(
       "floor_divide",
       [](const ScalarOrArray& a_, const ScalarOrArray& b_, StreamOrDevice s) {