diff --git a/RELEASE.md b/RELEASE.md
index 0a897d00..256cbb63 100644
--- a/RELEASE.md
+++ b/RELEASE.md
@@ -29,9 +29,14 @@
 - Added string representator in python for Cluster and Eta 
 - Added roi slice method in python for easy slicing of numpy arrays ``array[roi.slice()]``. 
 - added context manager for ``aare.RawMasterFile``
+- ``NDArray``/``NDView`` expressions now support scalar operands
+  (``2 * a + b / 4``) and can be assigned to an existing ``NDArray``, reusing
+  its buffer.
 
 ### API Changes:
 
+- ``NDArray`` ``+ - * /`` with a scalar now returns a lazy expression instead
+  of an ``NDArray`` and no longer converts the scalar to the element type.
 - ``FastPedestal`` variance is now a private ``double`` intermediate. Removed
   the C++ ``variance()``/``variance_unchecked()`` APIs and Python ``var()``.
   Standard deviation is calculated before conversion to the output type,
@@ -77,6 +82,8 @@
 
 ### Bugfixes:
 - Mismatched operators inhibited vectorization in gcc of NDArray math operators
+- ``NDArray``/``NDView`` math operators and expressions were not vectorized
+  for ``uint8_t`` and 64 bit integers, up to 30x slower than a plain loop.
 - Removed the move constructor and move assignment of
   ``ProducerConsumerQueue``. They left the moved-from queue with a null
   buffer, crashing its destructor if it still held elements, and leaked
diff --git a/benchmarks/ndarray_benchmark.cpp b/benchmarks/ndarray_benchmark.cpp
index 456f4d40..b3a3f610 100644
--- a/benchmarks/ndarray_benchmark.cpp
+++ b/benchmarks/ndarray_benchmark.cpp
@@ -130,4 +130,80 @@ BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State &st) {
     }
 }
 
+BENCHMARK_F(TwoArrays, ScalarMultiplyAddWithOperator)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        NDArray<int, 2> res = a * 2 + b;
+        benchmark::DoNotOptimize(res);
+    }
+}
+BENCHMARK_F(TwoArrays, ScalarMultiplyAddWithIndex)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        NDArray<int, 2> res(a.shape());
+        for (uint32_t i = 0; i < a.size(); i++) {
+            res(i) = a(i) * 2 + b(i);
+        }
+        benchmark::DoNotOptimize(res);
+    }
+}
+
+BENCHMARK_F(TwoArrays, AddToExistingWithOperator)(benchmark::State &st) {
+    NDArray<int, 2> res(a.shape());
+    for (auto _ : st) {
+        // This code gets timed
+        res = a + b;
+        benchmark::DoNotOptimize(res);
+    }
+}
+BENCHMARK_F(TwoArrays, AddToExistingWithIndex)(benchmark::State &st) {
+    NDArray<int, 2> res(a.shape());
+    for (auto _ : st) {
+        // This code gets timed
+        for (uint32_t i = 0; i < a.size(); i++) {
+            res(i) = a(i) + b(i);
+        }
+        benchmark::DoNotOptimize(res);
+    }
+}
+
+// A store to uint8_t or a 64 bit integer can alias the members of NDArray.
+// Compare with a loop over raw pointers to check that this does not block
+// vectorization.
+template <typename T> void InPlaceAddWithOperator(benchmark::State &st) {
+    NDArray<T, 2> a({size, size}, 1);
+    NDArray<T, 2> b({size, size}, 2);
+    for (auto _ : st) {
+        a += b;
+        benchmark::DoNotOptimize(a);
+    }
+}
+template <typename T> void InPlaceAddWithPointer(benchmark::State &st) {
+    NDArray<T, 2> a({size, size}, 1);
+    NDArray<T, 2> b({size, size}, 2);
+    for (auto _ : st) {
+        T *dst = a.data();
+        const T *src = b.data();
+        for (ssize_t i = 0, n = a.size(); i < n; ++i) {
+            dst[i] += src[i];
+        }
+        benchmark::DoNotOptimize(a);
+    }
+}
+template <typename T> void AddToExistingWithOperator(benchmark::State &st) {
+    NDArray<T, 2> a({size, size}, 1);
+    NDArray<T, 2> b({size, size}, 2);
+    NDArray<T, 2> res({size, size}, 0);
+    for (auto _ : st) {
+        res = a + b;
+        benchmark::DoNotOptimize(res);
+    }
+}
+BENCHMARK_TEMPLATE(InPlaceAddWithOperator, uint8_t);
+BENCHMARK_TEMPLATE(InPlaceAddWithPointer, uint8_t);
+BENCHMARK_TEMPLATE(AddToExistingWithOperator, uint8_t);
+BENCHMARK_TEMPLATE(InPlaceAddWithOperator, int64_t);
+BENCHMARK_TEMPLATE(InPlaceAddWithPointer, int64_t);
+BENCHMARK_TEMPLATE(AddToExistingWithOperator, int64_t);
+
 BENCHMARK_MAIN();
\ No newline at end of file
diff --git a/include/aare/ArrayExpr.hpp b/include/aare/ArrayExpr.hpp
index 8eab0d2e..8520208d 100644
--- a/include/aare/ArrayExpr.hpp
+++ b/include/aare/ArrayExpr.hpp
@@ -5,6 +5,8 @@
 #include <cassert>
 #include <cstddef>
 #include <cstdint>
+#include <functional>
+#include <type_traits>
 
 namespace aare {
 
@@ -20,80 +22,95 @@ template <typename E, ssize_t Ndim> class ArrayExpr {
     }
 };
 
-template <typename A, typename B, ssize_t Ndim>
-class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
-    const A &arr1_;
-    const B &arr2_;
-
-  public:
-    ArrayAdd(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
-        assert(arr1.size() == arr2.size());
-    }
-    auto operator[](size_t i) const { return arr1_[i] + arr2_[i]; }
-    size_t size() const { return arr1_.size(); }
-    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
-};
-
-template <typename A, typename B, ssize_t Ndim>
-class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
-    const A &arr1_;
-    const B &arr2_;
-
-  public:
-    ArraySub(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
-        assert(arr1.size() == arr2.size());
-    }
-    auto operator[](size_t i) const { return arr1_[i] - arr2_[i]; }
-    size_t size() const { return arr1_.size(); }
-    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
-};
-
-template <typename A, typename B, ssize_t Ndim>
-class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>, Ndim> {
-    const A &arr1_;
-    const B &arr2_;
-
-  public:
-    ArrayMul(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
-        assert(arr1.size() == arr2.size());
-    }
-    auto operator[](size_t i) const { return arr1_[i] * arr2_[i]; }
-    size_t size() const { return arr1_.size(); }
-    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
-};
-
-template <typename A, typename B, ssize_t Ndim>
-class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
-    const A &arr1_;
-    const B &arr2_;
-
-  public:
-    ArrayDiv(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
-        assert(arr1.size() == arr2.size());
-    }
-    auto operator[](size_t i) const { return arr1_[i] / arr2_[i]; }
-    size_t size() const { return arr1_.size(); }
-    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
-};
-
-template <typename A, typename B, ssize_t Ndim>
-auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
-    return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+/**
+ * @brief What an expression stores for one of its operands. An NDArray is
+ * replaced by a view, everything else is copied. The expression then holds its
+ * data pointers by value, so the compiler can keep them in registers even if
+ * the evaluation writes to a type that may alias them (e.g. uint8_t).
+ */
+template <typename E, ssize_t Ndim>
+auto operand(const ArrayExpr<E, Ndim> &expr) {
+    if constexpr (E::is_leaf)
+        return static_cast<const E &>(expr).view();
+    else
+        return static_cast<const E &>(expr);
 }
 
-template <typename A, typename B, ssize_t Ndim>
-auto operator-(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
-    return ArraySub<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+template <typename Op, typename A, typename B, ssize_t Ndim>
+class ArrayBinaryOp : public ArrayExpr<ArrayBinaryOp<Op, A, B, Ndim>, Ndim> {
+    A arr1_;
+    B arr2_;
+
+  public:
+    ArrayBinaryOp(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.shape() == arr2.shape());
+    }
+    auto operator[](size_t i) const { return Op{}(arr1_[i], arr2_[i]); }
+    size_t size() const { return arr1_.size(); }
+    std::array<ssize_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+/** @brief Scalar operand, takes size and shape from the expression it is
+ * combined with. */
+template <typename T, ssize_t Ndim>
+class ArrayScalar : public ArrayExpr<ArrayScalar<T, Ndim>, Ndim> {
+    T value_;
+    size_t size_;
+    std::array<ssize_t, Ndim> shape_;
+
+  public:
+    template <typename E>
+    ArrayScalar(T value, const ArrayExpr<E, Ndim> &like)
+        : value_(value), size_(like.size()), shape_(like.shape()) {}
+    T operator[](size_t) const { return value_; }
+    size_t size() const { return size_; }
+    std::array<ssize_t, Ndim> shape() const { return shape_; }
+};
+
+// Builds the expression for expr op expr, expr op scalar and scalar op expr
+template <typename Op, typename A, typename B, ssize_t Ndim>
+auto binary_op(const ArrayExpr<A, Ndim> &lhs, const ArrayExpr<B, Ndim> &rhs) {
+    auto a = operand(lhs);
+    auto b = operand(rhs);
+    return ArrayBinaryOp<Op, decltype(a), decltype(b), Ndim>(a, b);
 }
 
-template <typename A, typename B, ssize_t Ndim>
-auto operator*(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
-    return ArrayMul<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+template <typename Op, typename A, typename T, ssize_t Ndim,
+          typename = std::enable_if_t<std::is_arithmetic_v<T>>>
+auto binary_op(const ArrayExpr<A, Ndim> &lhs, T rhs) {
+    return binary_op<Op>(lhs, ArrayScalar<T, Ndim>(rhs, lhs));
 }
 
-template <typename A, typename B, ssize_t Ndim>
-auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
-    return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+template <typename Op, typename T, typename B, ssize_t Ndim,
+          typename = std::enable_if_t<std::is_arithmetic_v<T>>>
+auto binary_op(T lhs, const ArrayExpr<B, Ndim> &rhs) {
+    return binary_op<Op>(ArrayScalar<T, Ndim>(lhs, rhs), rhs);
+}
+
+// The operators only take part in overload resolution if binary_op accepts
+// the operands, meaning at least one of them is an ArrayExpr
+template <typename L, typename R>
+auto operator+(const L &lhs, const R &rhs)
+    -> decltype(binary_op<std::plus<>>(lhs, rhs)) {
+    return binary_op<std::plus<>>(lhs, rhs);
+}
+
+template <typename L, typename R>
+auto operator-(const L &lhs, const R &rhs)
+    -> decltype(binary_op<std::minus<>>(lhs, rhs)) {
+    return binary_op<std::minus<>>(lhs, rhs);
+}
+
+template <typename L, typename R>
+auto operator*(const L &lhs, const R &rhs)
+    -> decltype(binary_op<std::multiplies<>>(lhs, rhs)) {
+    return binary_op<std::multiplies<>>(lhs, rhs);
+}
+
+template <typename L, typename R>
+auto operator/(const L &lhs, const R &rhs)
+    -> decltype(binary_op<std::divides<>>(lhs, rhs)) {
+    return binary_op<std::divides<>>(lhs, rhs);
 }
 
 } // namespace aare
\ No newline at end of file
diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp
index fb6c3bb5..68f88ad6 100644
--- a/include/aare/NDArray.hpp
+++ b/include/aare/NDArray.hpp
@@ -27,6 +27,8 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     T *data_;
 
   public:
+    static constexpr bool is_leaf = true;
+
     ///////////////////////////////////////////////////////////////////////////////
     // Constructors
     //
@@ -135,9 +137,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      */
     template <typename E>
     NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
-        for (size_t i = 0; i < size_; ++i) {
-            data_[i] = expr[i];
-        }
+        *this = expr;
     }
 
     /**
@@ -324,6 +324,26 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         return *this;
     }
 
+    /**
+     * @brief Evaluate an ArrayExpr into the NDArray. Reuses the existing
+     * buffer if the shape matches.
+     */
+    template <typename E> NDArray &operator=(const ArrayExpr<E, Ndim> &expr) {
+        if (expr.shape() != shape_) {
+            // Evaluate before releasing the buffer, expr might refer to it
+            NDArray tmp(expr.shape());
+            tmp = expr;
+            return *this = std::move(tmp);
+        }
+        // Local copies, see elementwise()
+        const auto e = operand(expr);
+        T *data = data_;
+        for (size_t i = 0, n = size_; i < n; ++i) {
+            data[i] = e[i];
+        }
+        return *this;
+    }
+
     ///////////////////////////////////////////////////////////////////////////////
     // Math operators
     //
@@ -337,10 +357,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
             throw(std::runtime_error(
                 "Shape of NDArray must match for operator +="));
 
-        for (size_t i = 0; i < size_; ++i) {
-            data_[i] += other.data_[i];
-        }
-        return *this;
+        return elementwise(other.data_, std::plus<T>());
     }
 
     /**
@@ -351,10 +368,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
             throw(std::runtime_error(
                 "Shape of NDArray must match for operator -="));
 
-        for (size_t i = 0; i < size_; ++i) {
-            data_[i] -= other.data_[i];
-        }
-        return *this;
+        return elementwise(other.data_, std::minus<T>());
     }
 
     /**
@@ -365,10 +379,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
             throw(std::runtime_error(
                 "Shape of NDArray must match for operator *="));
 
-        for (size_t i = 0; i < size_; ++i) {
-            data_[i] *= other.data_[i];
-        }
-        return *this;
+        return elementwise(other.data_, std::multiplies<T>());
     }
 
     /**
@@ -379,13 +390,9 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      */
     template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
         // check shape
-        if (shape_ == other.shape()) {
-            for (size_t i = 0; i < size_; ++i) {
-                data_[i] /= other(i);
-            }
-            return *this;
-        }
-        throw(std::runtime_error("Shape of NDArray must match"));
+        if (shape_ != other.shape())
+            throw(std::runtime_error("Shape of NDArray must match"));
+        return elementwise(other.data(), std::divides<>());
     }
 
     /**
@@ -400,36 +407,28 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      * @brief Add a scalar value to all elements in the NDArray.
      */
     NDArray &operator+=(const T &value) {
-        for (size_t i = 0; i < size_; ++i)
-            data_[i] += value;
-        return *this;
+        return elementwise(value, std::plus<T>());
     }
 
     /**
      * @brief Subtract a scalar value to all elements in the NDArray.
      */
     NDArray &operator-=(const T &value) {
-        for (size_t i = 0; i < size_; ++i)
-            data_[i] -= value;
-        return *this;
+        return elementwise(value, std::minus<T>());
     }
 
     /**
      * @brief Multiply all elements in the NDArray with a scalar value
      */
     NDArray &operator*=(const T &value) {
-        for (size_t i = 0; i < size_; ++i)
-            data_[i] *= value;
-        return *this;
+        return elementwise(value, std::multiplies<T>());
     }
 
     /**
      * @brief Divide all elements in the NDArray with a scalar value
      */
     NDArray &operator/=(const T &value) {
-        for (size_t i = 0; i < size_; ++i)
-            data_[i] /= value;
-        return *this;
+        return elementwise(value, std::divides<T>());
     }
 
     /**
@@ -437,53 +436,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      * Used for example to mask out gain bits for Jungfrau detectors.
      */
     NDArray &operator&=(const T &mask) {
-        for (auto it = begin(); it != end(); ++it)
-            *it &= mask;
-        return *this;
-    }
-
-    /**
-     * @brief Operator +  with a scalar value. Returns a new NDArray.
-     *
-     * TODO! Expression template version of this?
-     */
-    NDArray operator+(const T &value) {
-        NDArray result = *this;
-        result += value;
-        return result;
-    }
-
-    /**
-     * @brief Operator -  with a scalar value. Returns a new NDArray.
-     *
-     * TODO! Expression template version of this?
-     */
-    NDArray operator-(const T &value) {
-        NDArray result = *this;
-        result -= value;
-        return result;
-    }
-
-    /**
-     * @brief Operator *  with a scalar value. Returns a new NDArray.
-     *
-     * TODO! Expression template version of this?
-     */
-    NDArray operator*(const T &value) {
-        NDArray result = *this;
-        result *= value;
-        return result;
-    }
-
-    /**
-     * @brief Operator /  with a scalar value. Returns a new NDArray.
-     *
-     * TODO! Expression template version of this?
-     */
-    NDArray operator/(const T &value) {
-        NDArray result = *this;
-        result /= value;
-        return result;
+        return elementwise(mask, std::bit_and<T>());
     }
 
     /**
@@ -509,19 +462,16 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      * @brief Compute the square root of all elements in the NDArray.
      */
     void sqrt() {
-        for (size_t i = 0; i < size_; ++i) {
-            data_[i] = std::sqrt(data_[i]);
+        T *data = data_;
+        for (size_t i = 0, n = size_; i < n; ++i) {
+            data[i] = std::sqrt(data[i]);
         }
     }
 
     /*
      * @brief Prefix increment operator. Increments all elements by 1.
      */
-    NDArray &operator++() {
-        for (size_t i = 0; i < size_; ++i)
-            data_[i] += T{1};
-        return *this;
-    }
+    NDArray &operator++() { return elementwise(T{1}, std::plus<T>()); }
 
     /** @brief Create a mutable view of the NDArray. */
     NDView<T, Ndim> view() { return NDView<T, Ndim>{data_, shape_}; }
@@ -532,6 +482,31 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     }
 
   private:
+    /**
+     * @brief Apply data_[i] = op(data_[i], value) to all elements.
+     *
+     * The members are copied to locals since a store to data_[i] may alias
+     * them for some T (uint8_t, int64_t, ...). The compiler then has to reload
+     * them every iteration, which blocks vectorization.
+     */
+    template <typename Op> NDArray &elementwise(T value, Op op) {
+        T *data = data_;
+        for (size_t i = 0, n = size_; i < n; ++i) {
+            data[i] = op(data[i], value);
+        }
+        return *this;
+    }
+
+    /** @brief Apply data_[i] = op(data_[i], other[i]) to all elements. */
+    template <typename V, typename Op>
+    NDArray &elementwise(const V *other, Op op) {
+        T *data = data_;
+        for (size_t i = 0, n = size_; i < n; ++i) {
+            data[i] = op(data[i], other[i]);
+        }
+        return *this;
+    }
+
     /**
      * @brief Reset the NDArray to an empty state. Dropping the ownership of
      * the data. Used internally for move operations to avoid double free or
diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp
index c507b41f..715127c9 100644
--- a/include/aare/NDView.hpp
+++ b/include/aare/NDView.hpp
@@ -173,8 +173,7 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
     }
 
     NDView &operator=(const T val) {
-        for (auto it = begin(); it != end(); ++it)
-            *it = val;
+        std::fill(begin(), end(), val);
         return *this;
     }
 
@@ -232,15 +231,19 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
 
     template <class BinaryOperation>
     NDView &elemenwise(T val, BinaryOperation op) {
-        for (uint64_t i = 0; i != size_; ++i) {
-            buffer_[i] = op(buffer_[i], val);
+        // Local copies since a store to buffer_[i] may alias the members
+        T *buffer = buffer_;
+        for (uint64_t i = 0, n = size_; i != n; ++i) {
+            buffer[i] = op(buffer[i], val);
         }
         return *this;
     }
     template <class BinaryOperation>
     NDView &elemenwise(const NDView &other, BinaryOperation op) {
-        for (uint64_t i = 0; i != size_; ++i) {
-            buffer_[i] = op(buffer_[i], other.buffer_[i]);
+        T *buffer = buffer_;
+        const T *other_buffer = other.buffer_;
+        for (uint64_t i = 0, n = size_; i != n; ++i) {
+            buffer[i] = op(buffer[i], other_buffer[i]);
         }
         return *this;
     }
diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp
index 3a3cafc0..3657c7b0 100644
--- a/src/NDArray.test.cpp
+++ b/src/NDArray.test.cpp
@@ -412,7 +412,7 @@ TEST_CASE("Elementwise operations on images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
         double v = 1.0;
-        auto C = A - v;
+        NDArray<double> C = A - v;
         REQUIRE(C.data() != A.data());
 
         // Value of C matches
@@ -429,7 +429,7 @@ TEST_CASE("Elementwise operations on images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
         double v = 1.0;
-        auto C = A + v;
+        NDArray<double> C = A + v;
         REQUIRE(C.data() != A.data());
 
         // Value of C matches
@@ -446,7 +446,7 @@ TEST_CASE("Elementwise operations on images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
         double v = 3.7;
-        auto C = A / v;
+        NDArray<double> C = A / v;
         REQUIRE(C.data() != A.data());
 
         // Value of C matches
@@ -463,12 +463,12 @@ TEST_CASE("Elementwise operations on images") {
         NDArray<double> A(shape, a_val);
         NDArray<double> B(shape, b_val);
         double v = 3.7;
-        auto C = A / v;
+        NDArray<double> C = A * v;
         REQUIRE(C.data() != A.data());
 
         // Value of C matches
         for (uint32_t i = 0; i < C.size(); ++i) {
-            REQUIRE(C(i) == a_val / v);
+            REQUIRE(C(i) == a_val * v);
         }
 
         // Value of A is not changed
@@ -478,6 +478,51 @@ TEST_CASE("Elementwise operations on images") {
     }
 }
 
+TEST_CASE("Expressions with scalars") {
+    NDArray<int> A({5, 5}, 3);
+    NDArray<int> B({5, 5}, 8);
+
+    NDArray<int> C = 2 * A + B / 4 - 1;
+    for (ssize_t i = 0; i < C.size(); ++i) {
+        REQUIRE(C(i) == 7);
+    }
+
+    // The scalar is not converted to the element type of the array
+    NDArray<double> D = A * 0.5;
+    for (ssize_t i = 0; i < D.size(); ++i) {
+        REQUIRE(D(i) == 1.5);
+    }
+}
+
+TEST_CASE("Assign an expression to an existing NDArray") {
+    NDArray<int> A({5, 5}, 3);
+    NDArray<int> B({5, 5}, 8);
+
+    SECTION("same shape reuses the buffer") {
+        NDArray<int> C({5, 5}, 0);
+        auto *ptr = C.data();
+        C = A + B;
+        REQUIRE(C.data() == ptr);
+        for (ssize_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == 11);
+        }
+    }
+    SECTION("the target can be part of the expression") {
+        A = A + B;
+        for (ssize_t i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == 11);
+        }
+    }
+    SECTION("different shape reallocates") {
+        NDArray<int> C({2, 3}, 0);
+        C = A + B;
+        REQUIRE(C.shape() == A.shape());
+        for (ssize_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == 11);
+        }
+    }
+}
+
 TEST_CASE("Assign an std::array to a 1D NDArray") {
     NDArray<int, 1> a{{5}, 0};
     std::array<int, 5> b{1, 2, 3, 4, 5};