diff --git a/src/algorithm/jde.cpp b/src/algorithm/jde.cpp
index 053c7c2a..5f3eb603 100644
--- a/src/algorithm/jde.cpp
+++ b/src/algorithm/jde.cpp
@@ -187,40 +187,78 @@ void jde::evolve(population &pop) const
 			}
 		}
 
+		std::vector<size_t> elements;
+		if (NP < 21) {
+			// Add elements in order.
+			for (size_t i = 0; i < NP; ++i) {
+				elements.push_back(i);
+			}
+		}
+
 		//Start of the loop through the deme
 		for (size_t i = 0; i < NP; ++i) {
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 2 !!!     */
-				r1 = p_idx();
-			} while (r1==i);
-
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 3 !!!     */
-				r2 = p_idx();
-			} while ((r2==i) || (r2==r1));
-
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 4 !!!     */
-				r3 = p_idx();
-			} while ((r3==i) || (r3==r1) || (r3==r2));
-
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 5 !!!     */
-				r4 = p_idx();
-			} while ((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));
-
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 6 !!!     */
-				r5 = p_idx();
-			} while ((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4));
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 7 !!!     */
-				r6 = p_idx();
-			} while ((r6==i) || (r6==r1) || (r6==r2) || (r6==r3) || (r6==r4) || (r6==r5));
-			do {                       /* Pick a random population member */
-				/* Endless loop for NP < 8 !!!     */
-				r7 = p_idx();
-			} while ((r7==i) || (r7==r1) || (r7==r2) || (r7==r3) || (r7==r4) || (r7==r5) || (r7==r6));
+			if (NP < 21) {
+				bool i_was_sampled = false;
+				size_t sampled_i_pos;
+				// Shuffle 8 elements.
+				// If 'i' is not picked, then we only use 7 of them.
+				for (size_t pos = 0; pos < 8; ++pos) {
+					// uniform_int<> rnd(a,b) gives values in the range [a,b] not [a,b).
+					boost::uniform_int<size_t> rnd_idx(0, NP - pos - 1);
+					size_t rnd = rnd_idx(m_urng);
+					if (elements[rnd] == i) {
+						i_was_sampled = true;
+						// Remember where we put i.
+						sampled_i_pos = NP - pos - 1;
+					}
+					std::swap(elements[rnd], elements[NP - pos - 1]);
+				}
+
+				if (i_was_sampled) {
+					std::swap(elements[sampled_i_pos], elements[NP - 7 - 1]);
+				}
+
+				r1 = elements[NP - 1];
+				r2 = elements[NP - 2];
+				r3 = elements[NP - 3];
+				r4 = elements[NP - 4];
+				r5 = elements[NP - 5];
+				r6 = elements[NP - 6];
+				r7 = elements[NP - 7];
+			} else {
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 2 !!!     */
+					r1 = p_idx();
+				} while (r1==i);
+
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 3 !!!     */
+					r2 = p_idx();
+				} while ((r2==i) || (r2==r1));
+
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 4 !!!     */
+					r3 = p_idx();
+				} while ((r3==i) || (r3==r1) || (r3==r2));
+
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 5 !!!     */
+					r4 = p_idx();
+				} while ((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));
+
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 6 !!!     */
+					r5 = p_idx();
+				} while ((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4));
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 7 !!!     */
+					r6 = p_idx();
+				} while ((r6==i) || (r6==r1) || (r6==r2) || (r6==r3) || (r6==r4) || (r6==r5));
+				do {                       /* Pick a random population member */
+					/* Endless loop for NP < 8 !!!     */
+					r7 = p_idx();
+				} while ((r7==i) || (r7==r1) || (r7==r2) || (r7==r3) || (r7==r4) || (r7==r5) || (r7==r6));
+			}
 
 			// Adapt amplification factor and crossover probability
 			double F=0, CR=0;