diff --git a/CHANGES.md b/CHANGES.md
index 9e4679a05..9fc60f454 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -5,7 +5,7 @@
 * SIRF/STIR
   - The implementation of the creation of `sirf.STIR.ImageData` from `sirf.STIR.AcquisitionData` has been revised to ensure compatibility of `ImageData` dimensions and voxel sizes with `AcquisitionData`.
 * Python interface
-  - Restored functionality for algebraic operations mixing `STIR.ImageData` and numpy arrays. (Note that sirf objects need to be on the "left" of the operation.)
+  - Restored functionality for algebraic operations mixing SIRF data containers and numpy arrays and corrected the description of the result type in User Guide.
 
 ## v3.9.0
 
diff --git a/doc/UserGuide.md b/doc/UserGuide.md
index 4272a1adf..9fe4c6302 100644
--- a/doc/UserGuide.md
+++ b/doc/UserGuide.md
@@ -110,7 +110,7 @@ Some classes are _derived_ from other classes, which means that they have (_inhe
 
 ### SIRF data algebra <a name="SIRF_data_algebra"></a>
 
-SIRF Python interface supports algebraic operations (`+`, `-`, `*` and `/`): e.g. elements of the data array stored in the object `a*b` are the products of the respective elements in `a` and `b`. Either or both `a` and `b` can be SIRF data objects of the same kind (either both `ImageData` or both `AcquisitionData`) or `numpy` arrays or scalars. One should be aware though that if `a` is a SIRF object then, just as one would expect, the product `a*b` will be a SIRF object of the same kind, but if `a` is a `numpy` object (array or scalar) then Python will try to convert `b` to a `numpy` object before computing `a*b`, and only if this fails it will compute `b*a` instead. To avoid confusion, the users are advised to check the type of `a*b` or, better still, always place a SIRF object on the left side of the product.
+SIRF Python interface supports algebraic operations (`+`, `-`, `*` and `/`): e.g. elements of the data array stored in the object `a*b` are the products of the respective elements in `a` and `b`. Either or both `a` and `b` can be SIRF data objects or `numpy` arrays or scalars. If both are SIRF data objects, they must be of the same type (e.g. either both `ImageData` or both `AcquisitionData`), and the result will be a SIRF data object of the same type, the same kind of result being produced if one is a SIRF data object and the other is either `numpy` array or scalar. If none is a SIRF data object, then the type of result will be determined by Python.
 
 ### Error handling <a name="error_handling"></a>
 
diff --git a/src/common/SIRF.py b/src/common/SIRF.py
index c126e9854..5a2940d80 100644
--- a/src/common/SIRF.py
+++ b/src/common/SIRF.py
@@ -89,6 +89,7 @@ class ArrayContainer(ABC):
     accessible via the address of the first data item.
     """
     warnings.simplefilter('once', DeprecationWarning)
+    __array_priority__ = 10
 
     def asarray(self, xp=numpy, copy=None, **kwargs):
         """Returns view (or fallback copy) of self"""
@@ -131,6 +132,10 @@ def __add__(self, other):
         '''
         return self.add(other)
 
+    def __radd__(self, other):
+        return self + other
+        #return other.add(self)
+
     def __sub__(self, other):
         '''
         Overloads - for data containers.
@@ -141,6 +146,9 @@ def __sub__(self, other):
         '''
         return self.subtract(other)
 
+    def __rsub__(self, other):
+        return -other + self
+
     def __mul__(self, other):
         '''
         Overloads * for data containers multiplication by a scalar or another
@@ -171,6 +179,9 @@ def __truediv__(self, other):
         '''
         return self.divide(other)
 
+    def __rtruediv__(self, other):
+        return other*self.power(-1)
+
     def __iadd__(self, other):
         self.add(other, out=self)
         return self
diff --git a/src/common/Utilities.py b/src/common/Utilities.py
index 7b44ff9a5..f20dad62c 100644
--- a/src/common/Utilities.py
+++ b/src/common/Utilities.py
@@ -790,24 +790,50 @@ def data_container_algebra_tests(test, x, eps=1e-4):
     t = numpy.linalg.norm(az)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    z = ay*x
+    az = z.as_array()
+    s = numpy.linalg.norm(az - ax * ay)
+    t = numpy.linalg.norm(az)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     y = x + 1
     ay = y.as_array()
     s = numpy.linalg.norm(ay - (ax + 1))
     t = numpy.linalg.norm(ay)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    y = 1 + x
+    test.check_if_equal(type(y), type(x))
+    ay = y.as_array()
+    s = numpy.linalg.norm(ay - (ax + 1))
+    t = numpy.linalg.norm(ay)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     y = x + ax
     ay = y.as_array()
     s = numpy.linalg.norm(ay - (ax + ax))
     t = numpy.linalg.norm(ay)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    y = ax + x
+    test.check_if_equal(type(y), type(x))
+    ay = y.as_array()
+    s = numpy.linalg.norm(ay - (ax + ax))
+    t = numpy.linalg.norm(ay)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     t = numpy.linalg.norm(ay)
     y = x - ax
     ay = y.as_array()
     s = numpy.linalg.norm(ay)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    y = ax - x
+    test.check_if_equal(type(y), type(x))
+    ay = y.as_array()
+    s = numpy.linalg.norm(ay)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     y *= 0
     x.add(1, out=y)
     ay = y.as_array()
@@ -827,12 +853,26 @@ def data_container_algebra_tests(test, x, eps=1e-4):
     t = numpy.linalg.norm(az)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    z = ay/x
+    test.check_if_equal(type(z), type(x))
+    az = z.as_array()
+    s = numpy.linalg.norm(az - ay/ax)
+    t = numpy.linalg.norm(az)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     z = x/2
     az = z.as_array()
     s = numpy.linalg.norm(az - ax/2)
     t = numpy.linalg.norm(az)
     test.check_if_zero_within_tolerance(s, eps * t)
 
+    z = 2/x
+    test.check_if_equal(type(z), type(x))
+    az = z.as_array()
+    s = numpy.linalg.norm(az - 2/ax)
+    t = numpy.linalg.norm(az)
+    test.check_if_zero_within_tolerance(s, eps * t)
+
     z *= 0
     x.divide(y, out=z)
     az = z.as_array()