diff --git a/src/hpcombi.cpp b/src/hpcombi.cpp
index 4aa2399d..30708cc0 100644
--- a/src/hpcombi.cpp
+++ b/src/hpcombi.cpp
@@ -569,7 +569,7 @@ The functionality described on this page is only available if
 
       thing.def(
           "__mul__",
-          // The next line is not a type, but is consistent with the
+          // The next line is not a typo, but is consistent with the
           // other transformations in libsemigroups_pybind11, since
           // function composition in HPCombi is backwards.
           [](PTransf16 const& x, PTransf16 const& y) { return y * x; },
@@ -1227,6 +1227,16 @@ The functionality described on this page is only available if
 
       thing.def("__copy__", [](Transf16 const& v) { return Transf16(v); });
 
+      thing.def(
+          "__mul__",
+          // The next line is not a typo, but is consistent with the
+          // other transformations in libsemigroups_pybind11, since
+          // function composition in HPCombi is backwards.
+          // Also this method is required because the return type of the one for
+          // PTransf16 is always PTransf16.
+          [](Transf16 const& x, Transf16 const& y) { return y * x; },
+          py::is_operator());
+
       ////////////////////////////////////////////////////////////////////////
       // Constructors
       ////////////////////////////////////////////////////////////////////////
@@ -1415,6 +1425,16 @@ The functionality described on this page is only available if
       thing.def("__repr__",
                 [](Perm16 const& self) { return repr(self, "Perm16"); });
 
+      thing.def(
+          "__mul__",
+          // The next line is not a typo, but is consistent with the
+          // other transformations in libsemigroups_pybind11, since
+          // function composition in HPCombi is backwards.
+          // Also this method is required because the return type of the one for
+          // PTransf16 is always PTransf16.
+          [](Perm16 const& x, Perm16 const& y) { return y * x; },
+          py::is_operator());
+
       ////////////////////////////////////////////////////////////////////////
       // Static methods
       ////////////////////////////////////////////////////////////////////////
@@ -2157,6 +2177,16 @@ The functionality described on this page is only available if
                                    "'PPerm16' and 'int");
       });
 
+      thing.def(
+          "__mul__",
+          // The next line is not a typo, but is consistent with the
+          // other transformations in libsemigroups_pybind11, since
+          // function composition in HPCombi is backwards.
+          // Also this method is required because the return type of the one for
+          // PTransf16 is always PTransf16.
+          [](PPerm16 const& x, PPerm16 const& y) { return y * x; },
+          py::is_operator());
+
       ////////////////////////////////////////////////////////////////////////
       // Static methods
       ////////////////////////////////////////////////////////////////////////
diff --git a/tests/test_hpcombi.py b/tests/test_hpcombi.py
index 0c2c488b..a26c46ca 100644
--- a/tests/test_hpcombi.py
+++ b/tests/test_hpcombi.py
@@ -306,6 +306,10 @@ def test_hpcombi_ptransf_nb_fix_points():
         assert PTransf16([1, 3, 2, 255, 10]).nb_fix_points() == 12
         assert PTransf16.one().nb_fix_points() == 16
 
+    def test_hpcombi_ptransf_mul_return_type():
+        x = PTransf16([1, 3, 2, 255, 10])
+        assert isinstance(x * x, PTransf16)
+
     ########################################################################
     # Transf16
     ########################################################################
@@ -361,6 +365,10 @@ def test_hpcombi_transf16_one():
         assert Transf16.one() == Transf16(list(range(16)))
         assert isinstance(Transf16.one(), Transf16)
 
+    def test_hpcombi_transf_mul_return_type():
+        x = Transf16([1, 0, 2])
+        assert isinstance(x * x, Transf16)
+
     ########################################################################
     # Perm16
     ########################################################################
@@ -497,6 +505,10 @@ def test_hpcombi_perm16_left_weak_leq_length():
     def test_hpcombi_perm16_unrankSJT():  # pylint: disable=invalid-name
         assert Perm16.unrankSJT(2) == Perm16([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 15, 13, 14])
 
+    def test_hpcombi_perm16_mul_return_type():
+        x = Perm16([1, 0, 2])
+        assert isinstance(x * x, Perm16)
+
     ########################################################################
     # PPerm16
     ########################################################################
@@ -561,3 +573,7 @@ def test_hpcombi_pperm16_inverse_ref():
         assert x * x.inverse_ref() == x.left_one()
         assert x.inverse_ref() * x == x.right_one()
         assert x**-1 == x.inverse_ref()
+
+    def test_hpcombi_pperm16_mul_return_type():
+        x = PPerm16([1, 0, 2])
+        assert isinstance(x * x, PPerm16)