diff --git a/quad.js b/quad.js
index c5fe18f..76e7420 100644
--- a/quad.js
+++ b/quad.js
@@ -93,10 +93,11 @@ module.exports = function cluster (srcPoints, options) {
 
 	return ids
 
-
-
-	// FIXME: it is possible to create one typed array heap and reuse that to avoid memory blow
 	function sort (x, y, diam, ids, level, group) {
+		return sortInner(x, y, diam, new Float64Array(ids), level, group)
+	}
+
+	function sortInner(x, y, diam, ids, level, group) {
 		if (!ids.length) return null
 
 		// save first point as level representative
@@ -132,22 +133,17 @@ module.exports = function cluster (srcPoints, options) {
 		let cx = x + d2, cy = y + d2
 
 		// distribute points by 4 buckets
-		let lolo = [], lohi = [], hilo = [], hihi = []
-
-		for (let i = 1, l = ids.length; i < l; i++) {
-			let idx = ids[i],
-				x = points[idx * 2],
-				y = points[idx * 2 + 1]
-			x < cx ? (y < cy ? lolo.push(idx) : lohi.push(idx)) : (y < cy ? hilo.push(idx) : hihi.push(idx))
-		}
+		const [lo, hi] = partition(ids.subarray(1), points, 0, cx) // partition by x
+		const [lolo, lohi] = partition(lo, points, 1, cy)          // then partition los and his by y
+		const [hilo, hihi] = partition(hi, points, 1, cy)
 
 		group <<= 2
 
 		sublevel.push(
-			sort(x, y, d2, lolo, level, group),
-			sort(x, cy, d2, lohi, level, group + 1),
-			sort(cx, y, d2, hilo, level, group + 2),
-			sort(cx, cy, d2, hihi, level, group + 3)
+			sortInner(x, y, d2, lolo, level, group),
+			sortInner(x, cy, d2, lohi, level, group + 1),
+			sortInner(cx, y, d2, hilo, level, group + 2),
+			sortInner(cx, cy, d2, hihi, level, group + 3)
 		)
 
 		return offset
@@ -284,3 +280,40 @@ function normalize (pts, bounds) {
 
 	return result
 }
+
+// Partitions ids in place into two subarrays,
+// one where the ids reference points which have values below the pivot,
+// and one where those values are >= the pivot.
+function partition(ids, points, offset, pivot) {
+	if (ids.length === 0) {
+		return [ids.subarray(0, 0), ids.subarray(0, 0)]
+	}
+
+	let next = ids[0]
+	let nextLo = 1
+	let nextHi = ids.length - 1
+	while (nextLo <= nextHi) {
+		if (points[next * 2 + offset] < pivot) {
+			let x = ids[nextLo]
+			ids[nextLo++] = next
+			next = x
+		} else {
+			let x = ids[nextHi]
+			ids[nextHi--] = next
+			next = x
+		}
+	}
+
+	if (next < pivot) {
+		ids[0] = next;
+	} else if (nextLo > 0) {
+		ids[0] = ids[nextHi];
+		ids[nextHi] = next;
+		nextHi--;
+		nextLo--;
+	}
+
+	let los = ids.subarray(0, nextLo)
+	let his = ids.subarray(nextLo, ids.length)
+	return [los, his]
+}
diff --git a/test.js b/test.js
index a495c6a..1057d2d 100644
--- a/test.js
+++ b/test.js
@@ -8,7 +8,7 @@ t('quad: offsets case', t => {
   let points = [.15,.8, .2,.15, .6,.6, .6,.45, .8,.1, .9,.6, .91,.61]
   let index = cluster(points, {bounds: [0,0,1,1]})
 
-  t.deepEqual(index.slice(), [0, 1,3,2, 4,5, 6])
+  t.deepEqual(index.slice(), [0, 1,4,6, 3,2,5])
 
   t.end()
 })
@@ -18,7 +18,7 @@ t('quad: output container', t => {
   let arr = []
   let index = cluster(points, {bounds: [0,0,1,1], output: arr})
 
-  t.deepEqual(arr.slice(), [0, 1,3,2, 4,5, 6])
+  t.deepEqual(arr.slice(), [0, 1,4,6, 3,2,5])
   t.equal(index, arr)
 
   t.end()
@@ -42,13 +42,13 @@ t('quad: selection', t => {
 
   let index = cluster(points)
 
-  t.deepEqual(index.range(), [0, 1, 2, 3, 4, 5, 6, 7])
-  t.deepEqual(index.range(1,1,6,6), [1, 2, 3, 4, 5, 6])
-  t.deepEqual(index.range(2,1,5,6), [2, 3, 4, 5])
-  t.deepEqual(index.range(1,2,6,5), [2, 3, 4, 5])
-  t.deepEqual(index.range(1,3,5,6), [3, 4, 5])
+  t.deepEqual(new Set(index.range()), new Set([0, 1, 2, 3, 4, 5, 6, 7]))
+  t.deepEqual(new Set(index.range(1,1,6,6)), new Set([1, 2, 3, 4, 5, 6]))
+  t.deepEqual(new Set(index.range(2,1,5,6)), new Set([2, 3, 4, 5]))
+  t.deepEqual(new Set(index.range(1,2,6,5)), new Set([2, 3, 4, 5]))
+  t.deepEqual(new Set(index.range(1,3,5,6)), new Set([3, 4, 5]))
 
-  t.deepEqual(index.range(5,6,1,3), [3, 4, 5])
+  t.deepEqual(new Set(index.range(5,6,1,3)), new Set([3, 4, 5]))
 
   t.end()
 })
@@ -287,4 +287,3 @@ t('performance', t => {
 
 	t.end()
 })
-