diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 118234bf1..166863e70 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -307,7 +307,8 @@ class Rotation:
             p_m = self.quaternion[...,1:]
             q_o = other.quaternion[...,0:1]
             p_o = other.quaternion[...,1:]
-            q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,)))
+            qmo = q_m*q_o
+            q = (qmo - np.einsum('...i,...i',p_m,p_o).reshape(qmo.shape))
             p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
             return self.copy(Rotation(np.block([q,p]))._standardize())
         else:
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 6e135212f..6b9cd2b1a 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -975,6 +975,13 @@ class TestRotation:
             assert np.allclose(rot_broadcast.quaternion[...,i,:], rot.quaternion)
 
 
+    @pytest.mark.parametrize('shape',[(3,2),(4,6)])
+    def test_broadcastcomposition(self,shape):
+        a = Rotation.from_random(shape[0])
+        b = Rotation.from_random(shape[1])
+        assert (a[:,np.newaxis]*b[np.newaxis,:]).allclose(a.broadcast_to(shape)*b.broadcast_to(shape))
+
+
     @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion,        np.array([-1,0,0,0])),
                                                  (Rotation.from_quaternion,        np.array([1,1,1,0])),
                                                  (Rotation.from_Euler_angles,      np.array([1,4,0])),