diff --git a/example/n_body.cc b/example/n_body.cc
index cb8eab5..6e0c410 100644
--- a/example/n_body.cc
+++ b/example/n_body.cc
@@ -11,9 +11,9 @@
 #include "soa/soa.h"
 
 using namespace std;
+using namespace ikra::soa;
 using ikra::executor::execute;
 using ikra::executor::Iterator;
-using ikra::soa::SoaLayout;
 
 static const int kNumBodies = 25;
 static const double kMaxMass = 1000;
@@ -27,7 +27,8 @@ static const int kMaxRect = 20;
 
 #define RAND (1.0 * rand() / RAND_MAX)
 
-static void render_rect(SDL_Renderer* renderer, double x, double y, double mass) {
+static void render_rect(SDL_Renderer* renderer,
+                        double x, double y, double mass) {
   SDL_Rect rect;
   rect.w = rect.h = mass / kMaxMass * kMaxRect;
   rect.x = (x/2 + 0.5) * kWindowWidth - rect.w/2;
@@ -57,12 +58,10 @@ class Body : public SoaLayout<Body, kNumBodies> {
     if (this == body) return;
 
     double EPS = 0.01;    // Softening parameter (just to avoid infinities).
-    double dx = body->position_[0] - position_[0];
-    double dy = body->position_[1] - position_[1];
-    double dist = sqrt(dx*dx + dy*dy);
+    Array<double, 2> d = body->position_ - position_;
+    double dist = sqrt(d[0]*d[0] + d[1]*d[1]);
     double F = kGravityConstant*mass_*body->mass_ / (dist*dist + EPS*EPS);
-    force_[0] += F*dx / dist;
-    force_[1] += F*dy / dist;
+    force_ += d * F / dist;
   }
 
   void add_force_to_all() {
@@ -75,10 +74,8 @@ class Body : public SoaLayout<Body, kNumBodies> {
   }
 
   void update(double dt) {
-    velocity_[0] += force_[0]*dt / mass_;
-    velocity_[1] += force_[1]*dt / mass_;
-    position_[0] += velocity_[0]*dt;
-    position_[1] += velocity_[1]*dt;
+    velocity_ += force_ * (dt / mass_);
+    position_ += velocity_ * dt;
 
     if (position_[0] < -1 || position_[0] > 1) {
       velocity_[0] = -velocity_[0];
diff --git a/ikra/soa/array_expression.h b/ikra/soa/array_expression.h
new file mode 100644
index 0000000..3717391
--- /dev/null
+++ b/ikra/soa/array_expression.h
@@ -0,0 +1,182 @@
+#ifndef SOA_ARRAY_EXPRESSION_H
+#define SOA_ARRAY_EXPRESSION_H
+
+#include <type_traits>
+
+#include "soa/array_traits.h"
+#include "soa/util.h"
+
+namespace ikra {
+namespace soa {
+
+// Superclass for all kinds of arrays (e.g., SoaArrayField_, ArrayAdd_).
+// `Self` is the type of the derived class (Curiously Recurring Template Pattern).
+template<typename Self>
+class ArrayExpression_ {
+ private:
+  using T = typename ArrayTraits<Self>::T;
+  static const int N = ArrayTraits<Self>::N;
+
+ public:
+  // Force size of this class to be 0.
+  char dummy_[0];
+
+  static const bool kIsArrayExpression = true;
+
+  T operator[](size_t i) const {
+    return static_cast<Self const&>(*this)[i];
+  }
+  T& operator[](size_t i) {
+    return static_cast<Self&>(*this)[i];
+  }
+  size_t size() const { return N; }
+
+  operator Self&() {
+    return static_cast<Self&>(*this);
+  }
+  operator Self const&() {
+    return static_cast<Self const&>(*this);
+  }
+};
+
+// These macros are used to overload assignment operators for array expressions.
+// ElementType and ArraySize should be defined within the context.
+#define IKRA_DEFINE_ARRAY_ASSIGNMENT(op) \
+  template<typename Self> \
+  void operator op##=(ArrayExpression_<Self> const& a) { \
+    for (size_t i = 0; i < ArraySize; ++i) { \
+      (*this)[i] op##= a[i]; \
+    } \
+  }
+#define IKRA_DEFINE_ARRAY_SCALAR_ASSIGNMENT(op) \
+  void operator op##=(ElementType c) { \
+    for (size_t i = 0; i < ArraySize; ++i) { \
+      (*this)[i] op##= c; \
+    } \
+  }
+
+// A wrapper for std::array to make it inherit from ArrayExpression_.
+template<typename ElementType, size_t ArraySize>
+class Array : public ArrayExpression_<Array<ElementType, ArraySize>> {
+ private:
+  std::array<ElementType, ArraySize> array_;
+
+ public:
+  ElementType operator[](size_t i) const { return array_[i]; }
+  ElementType& operator[](size_t i) { return array_[i]; }
+
+  Array() = default;
+  Array(std::array<ElementType, ArraySize> a) : array_(a) {}
+
+  // An array wrapper can be constructed from any ArrayExpression_,
+  // which forces its evaluation.
+  template<typename Self>
+  Array(ArrayExpression_<Self> const& a) {
+    static_assert(std::is_same<ElementType, typename ArrayTraits<Self>::T>::value,
+        "The specified element type is wrong.");
+    static_assert(ArraySize == ArrayTraits<Self>::N,
+        "The specified array size is wrong.");
+
+    for (size_t i = 0; i < ArraySize; ++i) {
+      array_[i] = a[i];
+    }
+  }
+
+  // Define the assignment operator between array wrappers and array expressions.
+  IKRA_DEFINE_ARRAY_ASSIGNMENT();
+  // Defines compound assignment operators between array wrappers and array expressions.
+  IKRA_APPLY_TO_ALL_OPERATORS(IKRA_DEFINE_ARRAY_ASSIGNMENT);
+  // Defines compound assignment operators between array wrappers and scalars.
+  IKRA_APPLY_TO_ALL_OPERATORS(IKRA_DEFINE_ARRAY_SCALAR_ASSIGNMENT);
+};
+
+// Generates classes and operators for the operations bewteen array expressions.
+// `name` will be part of the class name (e.g., Add -> ArrayAdd_),
+// and `op` is the corresponding operator.
+#define IKRA_DEFINE_ARRAY_OPERATION(name, op) \
+  template<typename Left, typename Right> \
+  class Array##name##_ : public ArrayExpression_<Array##name##_<Left, Right>> { \
+   private: \
+    Left  const& left_; \
+    Right const& right_; \
+  \
+    using T = typename ArrayTraits<Left>::T; \
+    static const int N = ArrayTraits<Left>::N; \
+  \
+   public: \
+    Array##name##_(Left const& l, Right const& r) : left_(l), right_(r) { \
+      static_assert(Left::kIsArrayExpression, \
+          "The left operand is not an array expression."); \
+      static_assert(Right::kIsArrayExpression, \
+          "The right operand is not an array expression."); \
+      static_assert(std::is_same<T, typename ArrayTraits<Right>::T>::value, \
+          "The element type of two operands are not the same."); \
+      static_assert(N == ArrayTraits<Right>::N, \
+          "The array sizes of two operands are not the same."); \
+    }; \
+    T operator[](size_t i) const { \
+      return left_[i] op right_[i]; \
+    } \
+  }; \
+  \
+  template<typename Left, typename Right, \
+           bool = Left::kIsArrayExpression, \
+           bool = Right::kIsArrayExpression> \
+  Array##name##_<Left, Right> operator op(Left const& l, Right const& r) { \
+    return Array##name##_<Left, Right>(l, r); \
+  }
+
+// Generates classes and operators for the scalar operation of array expressions.
+// The left template argument is an array while the right one is a scalar.
+// `name` will be part of the class name (e.g., Mul -> ArrayScalarMul_),
+// and `op` is the corresponding operator.
+#define IKRA_DEFINE_ARRAY_SCALAR_OPERATION(name, op) \
+  template<typename Left, typename Right> \
+  class ArrayScalar##name##_ : \
+      public ArrayExpression_<ArrayScalar##name##_<Left, Right>> { \
+   private: \
+    Left const& left_; \
+    Right right_; \
+  \
+   public: \
+    ArrayScalar##name##_(Left const& l, Right r) : left_(l), right_(r) { \
+      static_assert(Left::kIsArrayExpression, \
+          "The left operand is not an array expression."); \
+      static_assert(std::is_same<typename ArrayTraits<Left>::T, Right>::value, \
+          "The element type of two operands are not the same."); \
+    } \
+    Right operator[](size_t i) const { \
+      return left_[i] op right_; \
+    } \
+  }; \
+  \
+  template<typename Left, typename Right, bool = Left::kIsArrayExpression> \
+  typename std::enable_if< \
+    std::is_same<typename ArrayTraits<Left>::T, Right>::value, \
+    ArrayScalar##name##_<Left, Right> \
+  >::type operator op(Left const& l, Right r) { \
+    return ArrayScalar##name##_<Left, Right>(l, r); \
+  } \
+  template<typename Left, typename Right, bool = Right::kIsArrayExpression> \
+  typename std::enable_if< \
+    std::is_same<typename ArrayTraits<Right>::T, Left>::value, \
+    ArrayScalar##name##_<Right, Left> \
+  >::type operator op(Left l, Right const& r) { \
+    return ArrayScalar##name##_<Right, Left>(r, l); \
+  }
+
+IKRA_DEFINE_ARRAY_OPERATION(Add, +);
+IKRA_DEFINE_ARRAY_OPERATION(Sub, -);
+IKRA_DEFINE_ARRAY_OPERATION(Mul, *);
+IKRA_DEFINE_ARRAY_OPERATION(Div, /);
+IKRA_DEFINE_ARRAY_OPERATION(Mod, %);
+IKRA_DEFINE_ARRAY_SCALAR_OPERATION(Add, +);
+IKRA_DEFINE_ARRAY_SCALAR_OPERATION(Sub, -);
+IKRA_DEFINE_ARRAY_SCALAR_OPERATION(Mul, *);
+IKRA_DEFINE_ARRAY_SCALAR_OPERATION(Div, /);
+IKRA_DEFINE_ARRAY_SCALAR_OPERATION(Mod, %);
+
+}  // namespace soa
+}  // namespace ikra
+
+#endif  // SOA_ARRAY_EXPRESSION_H
diff --git a/ikra/soa/array_field.h b/ikra/soa/array_field.h
index f03a15f..ee18293 100644
--- a/ikra/soa/array_field.h
+++ b/ikra/soa/array_field.h
@@ -1,49 +1,63 @@
-#ifndef SOA_ARRAY_H
-#define SOA_ARRAY_H
+#ifndef SOA_ARRAY_FIELD_H
+#define SOA_ARRAY_FIELD_H
 
 #include <type_traits>
 
 #include "soa/constants.h"
 #include "soa/field.h"
+#include "soa/array_expression.h"
 
 namespace ikra {
 namespace soa {
 
-// Class for field declarations of type array. This class is intended to be
-// used with T = std::array and forwards all method invocations to the wrapped
-// array object. The array is stored in AoS format.
-template<typename T,
+// Class for field declarations of type array. This class is intended to
+// forwards all method invocations to the wrapped array object.
+// The array is stored in AoS format.
+template<typename ElementType,
+         size_t ArraySize,
          IndexType Capacity,
          uint32_t Offset,
          int AddressMode,
          int StorageMode,
          class Owner>
-class AosArrayField_ : public Field_<T, Capacity, Offset,
-                                     AddressMode, StorageMode, Owner> {
+class AosArrayField_ : public ArrayExpression_<AosArrayField_<ElementType,
+    ArraySize, Capacity, Offset, AddressMode, StorageMode, Owner>> {
+ private:
+  using T = std::array<ElementType, ArraySize>;
+
  public:
   static const int kSize = sizeof(T);
 
   // This operator is just for convenience reasons. The correct way to use it
   // would be "this->operator[](pos)".
-  typename T::reference operator[](typename T::size_type pos) {
+  typename T::reference operator[](typename T::size_type pos) const {
     return this->data_ptr()->operator[](pos);
   }
 
+  // Define the assignment operator between AOS array fields and array expressions.
+  IKRA_DEFINE_ARRAY_ASSIGNMENT();
+  // Defines compound assignment operators between AOS array fields and array expressions.
+  IKRA_APPLY_TO_ALL_OPERATORS(IKRA_DEFINE_ARRAY_ASSIGNMENT);
+  // Defines compound assignment operators between AOS array fields and scalars.
+  IKRA_APPLY_TO_ALL_OPERATORS(IKRA_DEFINE_ARRAY_SCALAR_ASSIGNMENT);
+
 #include "soa/field_shared.inc"
 };
 
 // Class for field declarations of type array. T is the base type of the array.
 // This class is the SoA counter part of AosArrayField_. Array slots are
 // layouted as if they were SoA fields (columns).
-template<typename T,
+template<typename ElementType,
          size_t ArraySize,
          IndexType Capacity,
          uint32_t Offset,
          int AddressMode,
          int StorageMode,
          class Owner>
-class SoaArrayField_ {
+class SoaArrayField_ : public ArrayExpression_<SoaArrayField_<ElementType,
+    ArraySize, Capacity, Offset, AddressMode, StorageMode, Owner>> {
  private:
+  using T = ElementType;
   using Self = SoaArrayField_<T, ArraySize, Capacity, Offset,
                               AddressMode, StorageMode, Owner>;
 
@@ -61,6 +75,13 @@ class SoaArrayField_ {
 
   operator T&() const = delete;
 
+  // Define the assignment operator between SOA array fields and array expressions.
+  IKRA_DEFINE_ARRAY_ASSIGNMENT();
+  // Defines compound assignment operators between SOA array fields and array expressions.
+  IKRA_APPLY_TO_ALL_OPERATORS(IKRA_DEFINE_ARRAY_ASSIGNMENT);
+  // Defines compound assignment operators between SOA array fields and scalars.
+  IKRA_APPLY_TO_ALL_OPERATORS(IKRA_DEFINE_ARRAY_SCALAR_ASSIGNMENT);
+
   // Implement std::array interface.
 
 #if defined(__CUDA_ARCH__) || !defined(__CUDACC__)
@@ -90,6 +111,7 @@ class SoaArrayField_ {
   __ikra_device__ T& back() const {
     return at<ArraySize - 1>();
   }
+
 #else
 
   // A helper class with an overridden operator= method. This class allows
@@ -265,4 +287,4 @@ class SoaArrayField_ {
 }  // namespace soa
 }  // namespace ikra
 
-#endif  // SOA_ARRAY_H
+#endif  // SOA_ARRAY_FIELD_H
diff --git a/ikra/soa/array_traits.h b/ikra/soa/array_traits.h
new file mode 100644
index 0000000..b73ffe1
--- /dev/null
+++ b/ikra/soa/array_traits.h
@@ -0,0 +1,65 @@
+#ifndef SOA_ARRAY_TRAITS_H
+#define SOA_ARRAY_TRAITS_H
+
+#include "soa/constants.h"
+
+namespace ikra {
+namespace soa {
+
+// Code adapted from:
+// https://stackoverflow.com/questions/45801696/no-type-named-type-in-ctrp-derived-class
+template<typename T> struct ArrayTraits;
+
+template<typename ElementType, size_t ArraySize> class Array;
+template<typename ElementType, size_t ArraySize>
+struct ArrayTraits<Array<ElementType, ArraySize>> {
+  using T = ElementType;
+  static const int N = ArraySize;
+};
+
+template<typename ElementType, size_t ArraySize, IndexType Capacity,
+         uint32_t Offset, int AddressMode, int StorageMode, class Owner>
+class AosArrayField_;
+template<typename ElementType, size_t ArraySize, IndexType Capacity,
+         uint32_t Offset, int AddressMode, int StorageMode, class Owner>
+struct ArrayTraits<AosArrayField_<ElementType, ArraySize, Capacity,
+                  Offset,AddressMode, StorageMode, Owner>> {
+  using T = ElementType;
+  static const int N = ArraySize;
+};
+
+template<typename ElementType, size_t ArraySize, IndexType Capacity,
+         uint32_t Offset, int AddressMode, int StorageMode, class Owner>
+class SoaArrayField_;
+template<typename ElementType, size_t ArraySize, IndexType Capacity,
+         uint32_t Offset, int AddressMode, int StorageMode, class Owner>
+struct ArrayTraits<SoaArrayField_<ElementType, ArraySize, Capacity,
+                  Offset,AddressMode, StorageMode, Owner>> {
+  using T = ElementType;
+  static const int N = ArraySize;
+};
+
+#define IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(name) \
+  template<typename Left, typename Right> class Array##name##_; \
+  \
+  template<typename Left, typename Right> \
+  struct ArrayTraits<Array##name##_<Left, Right>> { \
+    using T = typename ArrayTraits<Left>::T; \
+    static const int N = ArrayTraits<Left>::N; \
+  };
+
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(Add);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(Sub);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(Mul);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(Div);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(Mod);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(ScalarAdd);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(ScalarSub);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(ScalarMul);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(ScalarDiv);
+IKRA_DEFINE_TRAITS_OF_EXPRESSION_TEMPLATES(ScalarMod);
+
+}  // namespace soa
+}  // namespace ikra
+
+#endif  // SOA_ARRAY_TRAITS_H
diff --git a/ikra/soa/layout.h b/ikra/soa/layout.h
index 0f934ae..32bac8d 100644
--- a/ikra/soa/layout.h
+++ b/ikra/soa/layout.h
@@ -77,7 +77,7 @@ class SoaLayout : SizeNDummy<AddressMode> {
   // This struct serves as a namespace and contains array field types.
   struct array {
     template<typename T, size_t N, int Offset>
-    using aos = ikra::soa::AosArrayField_<std::array<T, N>, Capacity,
+    using aos = ikra::soa::AosArrayField_<T, N, Capacity,
                                           Offset, AddressMode,
                                           kStorageMode, Self>;
 
@@ -130,10 +130,11 @@ class SoaLayout : SizeNDummy<AddressMode> {
 
   // Create multiple new instances of this class. Data will be allocated inside
   // storage.data.
-  __ikra_device__ void* operator new[](size_t count) {
-    check_sizeof_class();
-    // "count" is the number of new instances.
-    return get_uninitialized(Self::storage().increase_size(count/AddressMode));
+  __ikra_device__ void* operator new[](size_t /*count*/) {
+    // check_sizeof_class();
+    // "count" is the number of bytes.
+    // return get_uninitialized(Self::storage().increase_size(count/AddressMode));
+    assert(false);    // TODO: Implement.
   }
 
   // Return the number of instances of this class.
diff --git a/ikra/soa/storage.h b/ikra/soa/storage.h
index 771a330..c0e016b 100644
--- a/ikra/soa/storage.h
+++ b/ikra/soa/storage.h
@@ -198,7 +198,7 @@ class DynamicStorage_
 
  private:
   template<typename StorageClass>
-  friend class StorageDataOffset;
+  friend struct StorageDataOffset;
 
   // The base pointer (start) of the arena.
   void* arena_base_;
@@ -248,7 +248,7 @@ class StaticStorage_
 
  private:
   template<typename StorageClass>
-  friend class StorageDataOffset;
+  friend struct StorageDataOffset;
 
   // Statically allocated data storage for the arena.
   char arena_base_[ArenaSize];
diff --git a/ikra/soa/util.h b/ikra/soa/util.h
index eaccfe4..ef68024 100644
--- a/ikra/soa/util.h
+++ b/ikra/soa/util.h
@@ -16,6 +16,10 @@ struct NvccWorkaroundIdentityClassHolder {
 // Only needed for clang.
 #define IKRA_fold(x) (__builtin_constant_p(x) ? (x) : (x))
 
+// Applies the macro to all the arithmetic operators.
+#define IKRA_APPLY_TO_ALL_OPERATORS(macro) \
+  macro(+); macro(-); macro(*); macro(/); macro(%);
+
 // This class is used to check if the compiler supports the selected addressing
 // mode.
 template<uintptr_t N>
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 9d5f2b4..b0c653e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -2,13 +2,18 @@ add_executable(array_test soa/array_test.cc)
 target_link_libraries(array_test gtest gtest_main ikra)
 add_test(array_test ${CMAKE_BINARY_DIR}/bin/array_test)
 
+add_executable(array_operations_test soa/array_operations_test.cc)
+target_link_libraries(array_operations_test gtest gtest_main ikra)
+add_test(array_operations_test ${CMAKE_BINARY_DIR}/bin/array_operations_test)
+
 add_executable(basic_class_test soa/basic_class_test.cc)
 target_link_libraries(basic_class_test gtest gtest_main ikra)
 add_test(basic_class_test ${CMAKE_BINARY_DIR}/bin/basic_class_test)
 
-add_executable(pointer_arithmetics_test soa/pointer_arithmetics_test.cc)
-target_link_libraries(pointer_arithmetics_test gtest gtest_main ikra)
-add_test(pointer_arithmetics_test ${CMAKE_BINARY_DIR}/bin/pointer_arithmetics_test)
+# TODO: The new[] operator does not work now.
+# add_executable(pointer_arithmetics_test soa/pointer_arithmetics_test.cc)
+# target_link_libraries(pointer_arithmetics_test gtest gtest_main ikra)
+# add_test(pointer_arithmetics_test ${CMAKE_BINARY_DIR}/bin/pointer_arithmetics_test)
 
 add_executable(executor_test executor_test.cc)
 target_link_libraries(executor_test gtest gtest_main ikra)
diff --git a/test/soa/array_operations_test.cc b/test/soa/array_operations_test.cc
new file mode 100644
index 0000000..9aaa291
--- /dev/null
+++ b/test/soa/array_operations_test.cc
@@ -0,0 +1,177 @@
+#include <numeric>
+
+#include "gtest/gtest.h"
+#include "soa/soa.h"
+
+namespace {
+using namespace ikra::soa;
+
+static const int kClassMaxInst = 1024;
+static const int kTestSize = 9;
+
+template<typename T, size_t N>
+T array_sum(Array<T, N> array) {
+  T sum = 0;
+  for (size_t i = 0; i < N; ++i) {
+    sum += array[i];
+  }
+  return sum;
+}
+
+// Zero addressing mode.
+#define IKRA_TEST_CLASSNAME TestClassZ
+#define IKRA_TEST_ADDRESS_MODE kAddressModeZero
+#include "array_test_layout.inc"
+#undef IKRA_TEST_CLASSNAME
+#undef IKRA_TEST_ADDRESS_MODE
+IKRA_HOST_STORAGE(TestClassZ)
+
+// Valid addressing mode.
+#define IKRA_TEST_CLASSNAME TestClassV
+#define IKRA_TEST_ADDRESS_MODE sizeof(int)
+#include "array_test_layout.inc"
+#undef IKRA_TEST_CLASSNAME
+#undef IKRA_TEST_ADDRESS_MODE
+IKRA_HOST_STORAGE(TestClassV)
+
+template<typename T>
+class ArrayTest : public testing::Test {};
+
+TYPED_TEST_CASE_P(ArrayTest);
+
+TYPED_TEST_P(ArrayTest, ArrayOperations) {
+  TypeParam::initialize_storage();
+  TypeParam** instances = new TypeParam*[kTestSize];
+  Array<Array<int, 3>, kTestSize> arrays;
+
+  // Initialize array fields and arrays.
+  for (int i = 0; i < kTestSize; ++i) {
+    instances[i] = new TypeParam();
+    for (int j = 0; j < 3; ++j) {
+      instances[i]->field1[j] = 2*i;
+      instances[i]->field3[j] = i + j;
+      arrays[i][j] = i - j;
+    }
+  }
+
+  // Test operations between two arrays.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> addition = arrays[i] + arrays[i];
+    EXPECT_EQ(array_sum(addition), 2 * array_sum(arrays[i]));
+    Array<int, 3> subtraction = addition - arrays[i];
+    EXPECT_EQ(array_sum(subtraction), array_sum(arrays[i]));
+    addition += arrays[i];
+    EXPECT_EQ(array_sum(addition), 3 * array_sum(arrays[i]));
+    subtraction -= arrays[i];
+    EXPECT_EQ(array_sum(subtraction), 0);
+  }
+
+  // Test operations between arrays and AOS array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> addition = arrays[i] + instances[i]->field1;
+    EXPECT_EQ(array_sum(addition), 9*i - 3);
+    Array<int, 3> subtraction = instances[i]->field1 - arrays[i];
+    EXPECT_EQ(array_sum(subtraction), 3*i + 3);
+    addition += instances[i]->field1;
+    EXPECT_EQ(array_sum(addition), 15*i - 3);
+    subtraction -= instances[i]->field1;
+    EXPECT_EQ(array_sum(subtraction), -3*i + 3);
+  }
+
+  // Test operations between arrays and SOA array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> addition = arrays[i] + instances[i]->field3;
+    EXPECT_EQ(array_sum(addition), 6*i);
+    Array<int, 3> subtraction = instances[i]->field3 - arrays[i];
+    EXPECT_EQ(array_sum(subtraction), 6);
+    addition += instances[i]->field3;
+    EXPECT_EQ(array_sum(addition), 9*i + 3);
+    subtraction -= instances[i]->field3;
+    EXPECT_EQ(array_sum(subtraction), -3*i + 3);
+  }
+
+  // Test operations between AOS and SOA array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> addition = instances[i]->field3 + instances[i]->field1;
+    EXPECT_EQ(array_sum(addition), 9*i + 3);
+    Array<int, 3> subtraction = instances[i]->field1 - instances[i]->field3;
+    EXPECT_EQ(array_sum(subtraction), 3*i - 3);
+  }
+
+  // Test operations between two AOS array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> addition = instances[i]->field1 + instances[i]->field1;
+    EXPECT_EQ(array_sum(addition), 12*i);
+    Array<int, 3> subtraction = instances[i]->field1 - instances[i]->field1;
+    EXPECT_EQ(array_sum(subtraction), 0);
+    instances[i]->field1 += instances[i]->field1;
+    EXPECT_EQ(instances[i]->field1_sum(), 12*i);
+    instances[i]->field1 -= instances[i]->field1;
+    EXPECT_EQ(instances[i]->field1_sum(), 0);
+  }
+
+  // Test operations between two SOA array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> addition = instances[i]->field3 + instances[i]->field3;
+    EXPECT_EQ(array_sum(addition), 6*i + 6);
+    Array<int, 3> subtraction = instances[i]->field3 - instances[i]->field3;
+    EXPECT_EQ(array_sum(subtraction), 0);
+    instances[i]->field3 += instances[i]->field3;
+    EXPECT_EQ(instances[i]->field3_sum(), 6*i + 6);
+    instances[i]->field3 -= instances[i]->field3;
+    EXPECT_EQ(instances[i]->field3_sum(), 0);
+  }
+}
+
+TYPED_TEST_P(ArrayTest, ArrayScalarOperations) {
+  TypeParam::initialize_storage();
+  TypeParam** instances = new TypeParam*[kTestSize];
+  Array<Array<int, 3>, kTestSize> arrays;
+
+  // Initialize array fields and arrays.
+  for (int i = 0; i < kTestSize; ++i) {
+    instances[i] = new TypeParam();
+    for (int j = 0; j < 3; ++j) {
+      instances[i]->field1[j] = 2*i;
+      instances[i]->field3[j] = i + j;
+      arrays[i][j] = i - j;
+    }
+  }
+
+  // Test multiplications of arrays.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> mul1 = arrays[i] * 3;
+    EXPECT_EQ(array_sum(mul1), 9*i - 9);
+    Array<int, 3> mul2 = 3 * arrays[i];
+    EXPECT_EQ(array_sum(mul2), 9*i - 9);
+    arrays[i] *= 2;
+    EXPECT_EQ(array_sum(arrays[i]), 6*i - 6);
+  }
+
+  // Test multiplications of AOS array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> mul1 = instances[i]->field1 * 3;
+    EXPECT_EQ(array_sum(mul1), 18*i);
+    Array<int, 3> mul2 = 3 * instances[i]->field1;
+    EXPECT_EQ(array_sum(mul2), 18*i);
+    instances[i]->field1 *= 2;
+    EXPECT_EQ(instances[i]->field1_sum(), 12*i);
+  }
+
+  // Test multiplications of AOS array fields.
+  for (int i = 0; i < kTestSize; ++i) {
+    Array<int, 3> mul1 = instances[i]->field3 * 3;
+    EXPECT_EQ(array_sum(mul1), 9*i + 9);
+    Array<int, 3> mul2 = 3 * instances[i]->field3;
+    EXPECT_EQ(array_sum(mul2), 9*i + 9);
+    instances[i]->field3 *= 2;
+    EXPECT_EQ(instances[i]->field3_sum(), 6*i + 6);
+  }
+}
+
+REGISTER_TYPED_TEST_CASE_P(ArrayTest, ArrayOperations, ArrayScalarOperations);
+
+INSTANTIATE_TYPED_TEST_CASE_P(Valid, ArrayTest, TestClassV);
+INSTANTIATE_TYPED_TEST_CASE_P(Zero, ArrayTest, TestClassZ);
+
+}  // namespace
diff --git a/test/soa/array_test.cc b/test/soa/array_test.cc
index 67a7fe6..d611425 100644
--- a/test/soa/array_test.cc
+++ b/test/soa/array_test.cc
@@ -45,8 +45,7 @@ TYPED_TEST_P(ArrayTest, AosArray) {
   for (int i = 0; i < kTestSize; ++i) {
     instances[i]->field1[0] = 1*i + 1;
     instances[i]->field1[1] = 2*i + 2;
-    // TODO: Should use "." for consistency.
-    instances[i]->field1->at(2) = 3*i + 3;
+    instances[i]->field1[2] = 3*i + 3;
   }
 
   // Check result.
@@ -73,7 +72,7 @@ TYPED_TEST_P(ArrayTest, SoaArray) {
   for (int i = 0; i < kTestSize; ++i) {
     instances[i]->field3[0] = 1*i + 1;
     instances[i]->field3[1] = 2*i + 2;
-    instances[i]->field3->at(2) = 3*i + 3;
+    instances[i]->field3[2] = 3*i + 3;
   }
 
   // Check result.
diff --git a/test/soa/benchmarks/build.sh b/test/soa/benchmarks/build.sh
index cd7ab71..24e3ded 100755
--- a/test/soa/benchmarks/build.sh
+++ b/test/soa/benchmarks/build.sh
@@ -64,6 +64,12 @@ do
   objdump -S bin/${out_name} > assembly/${out_name}.S
   echo -n "."
 
+  out_name="${v_compiler}_nbody_ikracpp_with_operations"
+  ${v_compiler} -O3 ${extra_args} nbody/ikracpp_with_operations.cc \
+      -std=c++11 -I../../../ikra -o bin/${out_name}
+  objdump -S bin/${out_name} > assembly/${out_name}.S
+  echo -n "."
+
   # No vectorization with: -fno-tree-vectorize
   out_name="${v_compiler}_nbody_soa"
   ${v_compiler} -O3 ${extra_args} nbody/soa.cc -std=c++11 -o bin/${out_name}
diff --git a/test/soa/benchmarks/nbody/ikracpp_with_operations.cc b/test/soa/benchmarks/nbody/ikracpp_with_operations.cc
new file mode 100644
index 0000000..4c564d4
--- /dev/null
+++ b/test/soa/benchmarks/nbody/ikracpp_with_operations.cc
@@ -0,0 +1,176 @@
+// N-Body Simulation
+// Code adapted from: http://physics.princeton.edu/~fpretori/Nbody/code.htm
+
+#define NDEBUG    // No asserts.
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+
+#include "benchmark.h"
+#include "executor/executor.h"
+#include "soa/soa.h"
+
+using namespace std;
+using namespace ikra::soa;
+using ikra::executor::execute;
+using ikra::executor::Iterator;
+
+static const int kIterations = 5;
+static const int kNumBodies = 8000;
+static const double kMaxMass = 1000;
+static const double kTimeInterval = 0.5;
+
+static const double kGravityConstant = 6.673e-11;   // gravitational constant
+
+#define RAND (1.0 * rand() / RAND_MAX)
+
+class Body : public SoaLayout<Body, kNumBodies> {
+ public:
+  IKRA_INITIALIZE_CLASS
+
+  Body(double mass, double pos_x, double pos_y, double vel_x, double vel_y)
+      : mass_(mass) {
+    position_[0] = pos_x;
+    position_[1] = pos_y;
+    velocity_[0] = vel_x;
+    velocity_[1] = vel_y;
+    this->reset_force();
+  }
+
+  double_ mass_;
+  array_(double, 2) position_;
+  array_(double, 2) velocity_;
+  array_(double, 2) force_;
+
+  void add_force(Body* body) {
+    if (this == body) return;
+
+    double EPS = 0.01;    // Softening parameter (just to avoid infinities).
+    Array<double, 2> d = body->position_ - position_;
+    double dist = sqrt(d[0]*d[0] + d[1]*d[1]);
+    double F = kGravityConstant*mass_*body->mass_ / (dist*dist + EPS*EPS);
+    force_ += d * F / dist;
+  }
+
+  void add_force_to_all() {
+    execute<kNumBodies>(&Body::add_force, this);
+  }
+
+  void reset_force() {
+    force_[0] = 0.0;
+    force_[1] = 0.0;
+  }
+
+  void update(double dt) {
+    velocity_ += force_ * (dt / mass_);
+    position_ += velocity_ * dt;
+
+    if (position_[0] < -1 || position_[0] > 1) {
+      velocity_[0] = -velocity_[0];
+    }
+    if (position_[1] < -1 || position_[1] > 1) {
+      velocity_[1] = -velocity_[1];
+    }
+  }
+
+  void codengen_simple_update(double dt) {
+    for (int i = 0; i < 100; ++i) {
+      velocity_ += force_ * (dt / mass_);
+      position_ += velocity_ * dt;
+    }
+  }
+};
+
+IKRA_HOST_STORAGE(Body)
+
+void instantiation() {
+  srand(42);
+
+  // Create objects.
+  for (int i = 0; i < kNumBodies; ++i) {
+    double mass = (RAND/2 + 0.5) * kMaxMass;
+    double pos_x = RAND*2 - 1;
+    double pos_y = RAND*2 - 1;
+    double vel_x = (RAND - 0.5) / 1000;
+    double vel_y = (RAND - 0.5) / 1000;
+    new Body(mass, pos_x, pos_y, vel_x, vel_y);
+  }
+}
+
+void run_simulation() {
+  for (int i = 0; i < kIterations; ++i) {
+    // Reset forces.
+    execute<kNumBodies>(&Body::reset_force);
+
+    // Update forces.
+    execute<kNumBodies>(&Body::add_force_to_all);
+
+    // Update velocities and positions.
+    execute<kNumBodies>(&Body::update, kTimeInterval);
+  }
+}
+
+void run_simple() {
+  for (int i = 0; i < kIterations*100; ++i) {
+    execute<kNumBodies>(&Body::codengen_simple_update, kTimeInterval);
+  }
+}
+
+int main() {
+    // Initialize object storage.
+  Body::initialize_storage();
+
+  uint64_t time_instantiation = measure<>::execution(instantiation);
+  uint64_t time_simulation = measure<>::execution(run_simulation);
+  uint64_t time_simple = measure<>::execution(run_simple);
+
+  // Calculate checksum
+  int checksum = 11;
+  for (uintptr_t i = 0; i < kNumBodies; i++) {
+    checksum += reinterpret_cast<int>(
+        r_float2int(Body::get(i)->position_[0]));
+    checksum += reinterpret_cast<int>(
+        r_float2int(Body::get(i)->position_[1]));
+    checksum = checksum % 1234567;
+
+    if (i < 10) {
+      printf("VALUE[%lu] = %f, %f\n", i,
+             (double) Body::get(i)->position_[0],
+             (double) Body::get(i)->position_[1]);
+    }
+  }
+
+  printf("instantiation: %lu\nsimulation: %lu\nsimple: %lu\n checksum: %i\n",
+         time_instantiation, time_simulation, time_simple, checksum);
+  return 0;
+}
+
+void codegen_reset_force(Body* body) {
+  body->reset_force();
+}
+
+void codegen_add_force(Body* self, Body* other) {
+  self->add_force(other);
+}
+
+void codegen_update(Body* self, double dt) {
+  self->update(dt);
+}
+
+void codegen_reset_force_manual(Body* body) {
+  body->force_[0] = 0.0;
+  body->force_[1] = 0.0;
+}
+
+void codegen_reset_force_half(Body* body) {
+  body->force_[0] = 0.0;
+}
+
+void codegen_reset_mass(Body* body) {
+  body->mass_ = 0.0;
+}
+
+void codengen_simple_update(Body* body, double dt) {
+  body->codengen_simple_update(dt);
+}