diff --git a/docs/manual/source/conf.py b/docs/manual/source/conf.py
index c525ee2..c10faeb 100644
--- a/docs/manual/source/conf.py
+++ b/docs/manual/source/conf.py
@@ -27,7 +27,7 @@ copyright = '2019-2025, Zoltán Vörös and contributors'
 author = 'Zoltán Vörös'
 
 # The full version, including alpha/beta/rc tags
-release = '6.7.3'
+release = '6.9.0'
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/docs/manual/source/ulab-ndarray.rst b/docs/manual/source/ulab-ndarray.rst
index 0403633..46288a2 100644
--- a/docs/manual/source/ulab-ndarray.rst
+++ b/docs/manual/source/ulab-ndarray.rst
@@ -1814,12 +1814,12 @@ array.
 Binary operators
 ================
 
-``ulab`` implements the ``+``, ``-``, ``*``, ``/``, ``**``, ``<``,
-``>``, ``<=``, ``>=``, ``==``, ``!=``, ``+=``, ``-=``, ``*=``, ``/=``,
-``**=`` binary operators, as well as the ``AND``, ``OR``, ``XOR``
-bit-wise operators that work element-wise. Note that the bit-wise
-operators will raise an exception, if either of the operands is of
-``float`` or ``complex`` type.
+``ulab`` implements the ``+``, ``-``, ``*``, ``/``, ``**``, ``%``,
+``<``, ``>``, ``<=``, ``>=``, ``==``, ``!=``, ``+=``, ``-=``, ``*=``,
+``/=``, ``**=``, ``%=`` binary operators, as well as the ``AND``,
+``OR``, ``XOR`` bit-wise operators that work element-wise. Note that the
+bit-wise operators will raise an exception, if either of the operands is
+of ``float`` or ``complex`` type.
 
 Broadcasting is available, meaning that the two operands do not even
 have to have the same shape. If the lengths along the respective axes
diff --git a/docs/ulab-convert.ipynb b/docs/ulab-convert.ipynb
index 1c7a843..53f28bc 100644
--- a/docs/ulab-convert.ipynb
+++ b/docs/ulab-convert.ipynb
@@ -14,7 +14,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2022-02-09T06:27:15.118699Z",
@@ -61,7 +61,7 @@
     "author = 'Zoltán Vörös'\n",
     "\n",
     "# The full version, including alpha/beta/rc tags\n",
-    "release = '6.7.3'\n",
+    "release = '6.9.0'\n",
     "\n",
     "\n",
     "# -- General configuration ---------------------------------------------------\n",
@@ -217,7 +217,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2022-02-09T06:27:21.647179Z",
@@ -258,7 +258,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2022-02-09T06:27:42.024028Z",
diff --git a/tests/2d/numpy/modulo.py b/tests/2d/numpy/modulo.py
new file mode 100644
index 0000000..b5171ba
--- /dev/null
+++ b/tests/2d/numpy/modulo.py
@@ -0,0 +1,26 @@
+try:
+    from ulab import numpy as np
+except:
+    import numpy as np
+
+
+dtypes = (np.uint8, np.int8, np.uint16, np.int16, np.float)
+
+for dtype1 in dtypes:
+    x1 = np.array(range(6), dtype=dtype1).reshape((2, 3))
+    for dtype2 in dtypes:
+        x2 = np.array(range(1, 4), dtype=dtype2)
+        print(x1 % x2)
+
+print()
+print('=' * 30)
+print('inplace modulo')
+print('=' * 30)
+print()
+
+for dtype1 in dtypes:
+    x1 = np.array(range(6), dtype=dtype1).reshape((2, 3))
+    for dtype2 in dtypes:
+        x2 = np.array(range(1, 4), dtype=dtype2)
+        x1 %= x2
+        print(x1)
diff --git a/tests/2d/numpy/modulo.py.exp b/tests/2d/numpy/modulo.py.exp
new file mode 100644
index 0000000..73aa463
--- /dev/null
+++ b/tests/2d/numpy/modulo.py.exp
@@ -0,0 +1,105 @@
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=uint8)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=uint16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int8)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 0, 1],
+       [0, 2, 0]], dtype=uint8)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=uint16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+
+==============================
+inplace modulo
+==============================
+
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=uint8)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0, 0, 1],
+       [0, 2, 0]], dtype=uint8)
+array([[0, 0, 1],
+       [0, 0, 0]], dtype=int16)
+array([[0.0, 0.0, 1.0],
+       [0.0, 0.0, 0.0]], dtype=float64)
+array([[0.0, 0.0, 1.0],
+       [0.0, 0.0, 0.0]], dtype=float64)
+array([[0.0, 0.0, 1.0],
+       [0.0, 0.0, 0.0]], dtype=float64)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0, 1, 2],
+       [0, 0, 2]], dtype=int16)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)
+array([[0.0, 1.0, 2.0],
+       [0.0, 0.0, 2.0]], dtype=float64)