diff --git a/src/tools.c b/src/tools.c
index d25b7401a1e0515c650333b41193d54b5e155d39..7b468426b894044dc8e215f5150d9536fa2d4cdd 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -59,7 +59,6 @@ void factor(int value, int *f1, int *f2) {
 
 void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
               int periodic) {
-
   int i, j, k, count = 0;
   // int mj, mk;
   // double maxratio = 1.0;
@@ -124,7 +123,6 @@ void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
 void pairs_single_density(double *dim, long long int pid,
                           struct part *__restrict__ parts, int N,
                           int periodic) {
-
   int i, k;
   // int mj, mk;
   // double maxratio = 1.0;
@@ -177,19 +175,16 @@ void pairs_single_density(double *dim, long long int pid,
 }
 
 void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
-
   float r2, hi, hj, hig2, hjg2, dx[3];
   struct part *pi, *pj;
 
   /* Implements a double-for loop and checks every interaction */
   for (int i = 0; i < ci->count; ++i) {
-
     pi = &ci->parts[i];
     hi = pi->h;
     hig2 = hi * hi * kernel_gamma2;
 
     for (int j = 0; j < cj->count; ++j) {
-
       pj = &cj->parts[j];
 
       /* Pairwise distance */
@@ -201,7 +196,6 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
 
       /* Hit or miss? */
       if (r2 < hig2) {
-
         /* Interact */
         runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj);
       }
@@ -210,13 +204,11 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
 
   /* Reverse double-for loop and checks every interaction */
   for (int j = 0; j < cj->count; ++j) {
-
     pj = &cj->parts[j];
     hj = pj->h;
     hjg2 = hj * hj * kernel_gamma2;
 
     for (int i = 0; i < ci->count; ++i) {
-
       pi = &ci->parts[i];
 
       /* Pairwise distance */
@@ -228,7 +220,6 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
 
       /* Hit or miss? */
       if (r2 < hjg2) {
-
         /* Interact */
         runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi);
       }
@@ -242,13 +233,11 @@ void self_all_density(struct runner *r, struct cell *ci) {
 
   /* Implements a double-for loop and checks every interaction */
   for (int i = 0; i < ci->count; ++i) {
-
     pi = &ci->parts[i];
     hi = pi->h;
     hig2 = hi * hi * kernel_gamma2;
 
     for (int j = i + 1; j < ci->count; ++j) {
-
       pj = &ci->parts[j];
       hj = pj->h;
       hjg2 = hj * hj * kernel_gamma2;
@@ -264,14 +253,12 @@ void self_all_density(struct runner *r, struct cell *ci) {
 
       /* Hit or miss? */
       if (r2 < hig2) {
-
         /* Interact */
         runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj);
       }
 
       /* Hit or miss? */
       if (r2 < hjg2) {
-
         dxi[0] = -dxi[0];
         dxi[1] = -dxi[1];
         dxi[2] = -dxi[2];
@@ -285,7 +272,6 @@ void self_all_density(struct runner *r, struct cell *ci) {
 
 void pairs_single_grav(double *dim, long long int pid,
                        struct gpart *__restrict__ parts, int N, int periodic) {
-
   int i, k;
   // int mj, mk;
   // double maxratio = 1.0;
@@ -344,7 +330,6 @@ void pairs_single_grav(double *dim, long long int pid,
  */
 
 void density_dump(int N) {
-
   int k;
   float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
   struct part /**pi[4],  *pj[4],*/ Pi[4], Pj[4];
@@ -382,7 +367,6 @@ void density_dump(int N) {
 void engine_single_density(double *dim, long long int pid,
                            struct part *__restrict__ parts, int N,
                            int periodic) {
-
   int i, k;
   double r2, dx[3];
   float fdx[3], ih;
@@ -429,7 +413,6 @@ void engine_single_density(double *dim, long long int pid,
 
 void engine_single_force(double *dim, long long int pid,
                          struct part *__restrict__ parts, int N, int periodic) {
-
   int i, k;
   double r2, dx[3];
   float fdx[3];
@@ -471,3 +454,29 @@ void engine_single_force(double *dim, long long int pid,
           p.a_hydro[1], p.a_hydro[2]);
   fflush(stdout);
 }
+
+/**
+ * Returns a random number (uniformly distributed) in [a,b[
+ */
+double random_uniform(double a, double b) {
+  return (rand() / (double)RAND_MAX) * (b - a) + a;
+}
+
+/**
+ * @brief Randomly shuffle an array of particles.
+ */
+void shuffle_particles(struct part *parts, const int count) {
+  if (count > 1) {
+    for (int i = 0; i < count - 1; i++) {
+      int j = i + random_uniform(0., (double)(count - 1 - i));
+
+      struct part particle = parts[j];
+
+      parts[j] = parts[i];
+
+      parts[i] = particle;
+    }
+
+  } else
+    error("Array not big enough to shuffle!");
+}
diff --git a/src/tools.h b/src/tools.h
index ccffc77ceb8a967fd40c3737651ba75d529eee0f..8e0212652922fb4afd1ac89a64d4a01c71ecd4a9 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -38,4 +38,7 @@ void self_all_density(struct runner *r, struct cell *ci);
 void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
               int periodic);
 
+double random_uniform(double a, double b);
+void shuffle_particles(struct part *parts, const int count);
+
 #endif /* SWIFT_TOOL_H */
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 7915511eed50a229a94eda6bb338607099303421..9b18dd19231fc4925d8e95ad197ff36b0ea5a105 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -31,14 +31,6 @@ enum velocity_types {
   velocity_rotating
 };
 
-/**
- * @brief Returns a random number (uniformly distributed) in [a,b[
- */
-double random_uniform(double a, double b) {
-  return (rand() / (double)RAND_MAX) * (b - a) + a;
-}
-
-
 /**
  * @brief Constructs a cell and all of its particle in a valid state prior to
  * a DOPAIR or DOSELF calcuation.
@@ -46,10 +38,12 @@ double random_uniform(double a, double b) {
  * @param n The cube root of the number of particles.
  * @param offset The position of the cell offset from (0,0,0).
  * @param size The cell size.
- * @param h The smoothing length of the particles in units of the inter-particle separation.
+ * @param h The smoothing length of the particles in units of the inter-particle
+ *separation.
  * @param density The density of the fluid.
  * @param partId The running counter of IDs.
- * @param pert The perturbation to apply to the particles in the cell in units of the inter-particle separation.
+ * @param pert The perturbation to apply to the particles in the cell in units
+ *of the inter-particle separation.
  * @param vel The type of velocity field (0, random, divergent, rotating)
  */
 struct cell *make_cell(size_t n, double *offset, double size, double h,
@@ -127,6 +121,8 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
   cell->ti_end_min = 1;
   cell->ti_end_max = 1;
 
+  shuffle_particles(cell->parts, cell->count);
+
   cell->sorted = 0;
   cell->sort = NULL;
   cell->sortsize = 0;
@@ -145,7 +141,6 @@ void clean_up(struct cell *ci) {
  * @brief Initializes all particles field to be ready for a density calculation
  */
 void zero_particle_fields(struct cell *c) {
-
   for (size_t pid = 0; pid < c->count; pid++) {
     c->parts[pid].rho = 0.f;
     c->parts[pid].rho_dh = 0.f;
@@ -157,7 +152,6 @@ void zero_particle_fields(struct cell *c) {
  * @brief Ends the loop by adding the appropriate coefficients
  */
 void end_calculation(struct cell *c) {
-
   for (size_t pid = 0; pid < c->count; pid++) {
     hydro_end_density(&c->parts[pid], 1);
   }
@@ -168,7 +162,6 @@ void end_calculation(struct cell *c) {
  */
 void dump_particle_fields(char *fileName, struct cell *main_cell,
                           struct cell **cells) {
-
   FILE *file = fopen(fileName, "w");
 
   /* Write header */
@@ -205,7 +198,6 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
   for (int i = 0; i < 3; ++i) {
     for (int j = 0; j < 3; ++j) {
       for (int k = 0; k < 3; ++k) {
-
         struct cell *cj = cells[i * 9 + j * 3 + k];
         if (cj == main_cell) continue;
 
@@ -242,7 +234,6 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
 
 /* And go... */
 int main(int argc, char *argv[]) {
-
   size_t runs = 0, particles = 0;
   double h = 1.2348, size = 1., rho = 1.;
   double perturbation = 0.;
@@ -336,7 +327,6 @@ int main(int argc, char *argv[]) {
   for (int i = 0; i < 3; ++i) {
     for (int j = 0; j < 3; ++j) {
       for (int k = 0; k < 3; ++k) {
-
         double offset[3] = {i * size, j * size, k * size};
         cells[i * 9 + j * 3 + k] = make_cell(particles, offset, size, h, rho,
                                              &partId, perturbation, vel);
@@ -349,7 +339,6 @@ int main(int argc, char *argv[]) {
 
   ticks time = 0;
   for (size_t i = 0; i < runs; ++i) {
-
     /* Zero the fields */
     for (int j = 0; j < 27; ++j) zero_particle_fields(cells[j]);
 
diff --git a/tests/testPair.c b/tests/testPair.c
index 6e46b577ca63a8d3c2edce888a7485af0949813d..07e576332f81b0471f35f09dfcda21ad535adc8b 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -24,13 +24,6 @@
 #include <unistd.h>
 #include "swift.h"
 
-/**
- * Returns a random number (uniformly distributed) in [a,b[
- */
-double random_uniform(double a, double b) {
-  return (rand() / (double)RAND_MAX) * (b - a) + a;
-}
-
 /* n is both particles per axis and box size:
  * particles are generated on a mesh with unit spacing
  */
@@ -93,6 +86,8 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
   cell->ti_end_min = 1;
   cell->ti_end_max = 1;
 
+  shuffle_particles(cell->parts, cell->count);
+
   cell->sorted = 0;
   cell->sort = NULL;
   cell->sortsize = 0;
@@ -111,7 +106,6 @@ void clean_up(struct cell *ci) {
  * @brief Initializes all particles field to be ready for a density calculation
  */
 void zero_particle_fields(struct cell *c) {
-
   for (size_t pid = 0; pid < c->count; pid++) {
     c->parts[pid].rho = 0.f;
     c->parts[pid].rho_dh = 0.f;
@@ -123,7 +117,6 @@ void zero_particle_fields(struct cell *c) {
  * @brief Dump all the particles to a file
  */
 void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
-
   FILE *file = fopen(fileName, "w");
 
   /* Write header */
@@ -254,7 +247,6 @@ int main(int argc, char *argv[]) {
 
   time = 0;
   for (size_t i = 0; i < runs; ++i) {
-
     /* Zero the fields */
     zero_particle_fields(ci);
     zero_particle_fields(cj);