diff --git a/benchmark/mergeSortedXY.ts b/benchmark/mergeSortedXY.ts
new file mode 100644
index 000000000..c69a532c4
--- /dev/null
+++ b/benchmark/mergeSortedXY.ts
@@ -0,0 +1,195 @@
+/* eslint-disable no-console */
+import Benchmark from 'benchmark';
+import type { DataXY } from 'cheminfo-types';
+import { XSadd } from 'ml-xsadd';
+
+import { xGetSortOrder } from '../src/x/xGetSortOrder.ts';
+
+// The two strategies of mergeSortedXY, copied here so one process compares them on
+// the same data. The heap costs O(log k) per point, ordering the concatenated
+// points costs the same whatever the number of spectra.
+
+const TOTAL_POINTS = 2_000_000;
+
+/**
+ * Counts the points of every spectrum.
+ * @param data - spectra.
+ * @returns the number of points of all of them.
+ */
+function getTotalLength(data: DataXY[]): number {
+  let total = 0;
+  for (const spectrum of data) {
+    total += spectrum.x.length;
+  }
+  return total;
+}
+
+/**
+ * Merges by repeatedly taking the smallest x of a binary heap of the spectra.
+ * @param data - spectra, each one with ascending x values.
+ * @returns every point of every spectrum, ascending by x.
+ */
+function mergeByHeap(data: DataXY[]) {
+  const total = getTotalLength(data);
+  const x = new Float64Array(total);
+  const y = new Float64Array(total);
+  const spectrumIndex = new Uint32Array(total);
+
+  const ids = new Uint32Array(data.length);
+  const keys = new Float64Array(data.length);
+  const positions = new Uint32Array(data.length);
+  let size = 0;
+  for (let i = 0; i < data.length; i++) {
+    const spectrum = data[i];
+    if (spectrum.x.length === 0) continue;
+    const key = spectrum.x[0];
+    let child = size++;
+    while (child > 0) {
+      const parent = (child - 1) >> 1;
+      if (keys[parent] <= key) break;
+      keys[child] = keys[parent];
+      ids[child] = ids[parent];
+      child = parent;
+    }
+    keys[child] = key;
+    ids[child] = i;
+  }
+
+  let at = 0;
+  while (size > 0) {
+    const id = ids[0];
+    const spectrum = data[id];
+    const position = positions[id];
+    const currentX = keys[0];
+    x[at] = currentX;
+    y[at] = spectrum.y[position];
+    spectrumIndex[at] = id;
+    at++;
+
+    const next = position + 1;
+    positions[id] = next;
+    let key: number;
+    let newId: number;
+    if (next < spectrum.x.length) {
+      key = spectrum.x[next];
+      newId = id;
+    } else {
+      size--;
+      if (size === 0) break;
+      key = keys[size];
+      newId = ids[size];
+    }
+
+    let parent = 0;
+    for (;;) {
+      const left = 2 * parent + 1;
+      if (left >= size) break;
+      const right = left + 1;
+      const child = right < size && keys[right] < keys[left] ? right : left;
+      if (keys[child] >= key) break;
+      keys[parent] = keys[child];
+      ids[parent] = ids[child];
+      parent = child;
+    }
+    keys[parent] = key;
+    ids[parent] = newId;
+  }
+
+  return { x, y, spectrumIndex };
+}
+
+/**
+ * Merges by concatenating the spectra and ordering the result by x.
+ * @param data - spectra, each one with ascending x values.
+ * @returns every point of every spectrum, ascending by x.
+ */
+function mergeBySorting(data: DataXY[]) {
+  const total = getTotalLength(data);
+  const concatenatedX = new Float64Array(total);
+  const concatenatedY = new Float64Array(total);
+  const concatenatedIndex = new Uint32Array(total);
+  let at = 0;
+  for (let i = 0; i < data.length; i++) {
+    const spectrum = data[i];
+    concatenatedX.set(spectrum.x, at);
+    concatenatedY.set(spectrum.y, at);
+    concatenatedIndex.fill(i, at, at + spectrum.x.length);
+    at += spectrum.x.length;
+  }
+
+  const order = xGetSortOrder(concatenatedX);
+  const x = new Float64Array(total);
+  const y = new Float64Array(total);
+  const spectrumIndex = new Uint32Array(total);
+  for (let i = 0; i < total; i++) {
+    const from = order[i];
+    x[i] = concatenatedX[from];
+    y[i] = concatenatedY[from];
+    spectrumIndex[i] = concatenatedIndex[from];
+  }
+
+  return { x, y, spectrumIndex };
+}
+
+/**
+ * Builds centroided spectra sharing a peak list, each one drifting slightly.
+ * @param numberOfSpectra - how many spectra.
+ * @param peaksPerSpectrum - how many peaks in each.
+ * @returns the spectra, each with ascending x values.
+ */
+function makeSpectra(
+  numberOfSpectra: number,
+  peaksPerSpectrum: number,
+): DataXY[] {
+  const { random } = new XSadd(42);
+  const spacing = 1900 / peaksPerSpectrum;
+  const data: DataXY[] = [];
+  for (let spectrum = 0; spectrum < numberOfSpectra; spectrum++) {
+    const x = new Float64Array(peaksPerSpectrum);
+    const y = new Float64Array(peaksPerSpectrum);
+    for (let i = 0; i < peaksPerSpectrum; i++) {
+      x[i] = 100 + i * spacing + (random() - 0.5) * spacing * 0.4;
+      y[i] = 10 ** (random() * 5);
+    }
+    data.push({ x, y });
+  }
+  return data;
+}
+
+console.log(`${TOTAL_POINTS} points, split over a growing number of spectra`);
+for (const numberOfSpectra of [5, 15, 30, 100, 1000, 5000]) {
+  const peaksPerSpectrum = Math.round(TOTAL_POINTS / numberOfSpectra);
+  const data = makeSpectra(numberOfSpectra, peaksPerSpectrum);
+  const total = numberOfSpectra * peaksPerSpectrum;
+
+  const fromHeap = mergeByHeap(data);
+  const fromSorting = mergeBySorting(data);
+  let identical = true;
+  for (let i = 0; i < total; i++) {
+    if (fromHeap.x[i] !== fromSorting.x[i]) {
+      identical = false;
+      break;
+    }
+  }
+
+  const nanoseconds: Record<string, number> = {};
+  new Benchmark.Suite()
+    .add('heap', () => mergeByHeap(data), { minSamples: 30 })
+    .add('sorting', () => mergeBySorting(data), { minSamples: 30 })
+    .on('cycle', (event: Benchmark.Event) => {
+      const { target } = event;
+      const { stats, name } = target;
+      if (!stats || !name) return;
+      nanoseconds[name] = (stats.mean * 1e9) / total;
+    })
+    .run();
+
+  const { heap, sorting } = nanoseconds;
+  console.log(
+    `${String(numberOfSpectra).padStart(5)} spectra x ${String(peaksPerSpectrum).padStart(6)} peaks` +
+      `  same x: ${identical}` +
+      `  heap ${heap.toFixed(1).padStart(5)} ns/point` +
+      `  sorting ${sorting.toFixed(1).padStart(5)} ns/point` +
+      `  ${heap > sorting ? `sorting ${(heap / sorting).toFixed(2)}x` : `heap ${(sorting / heap).toFixed(2)}x`}`,
+  );
+}
diff --git a/benchmark/xGetSortOrder.ts b/benchmark/xGetSortOrder.ts
new file mode 100644
index 000000000..519cf2c2a
--- /dev/null
+++ b/benchmark/xGetSortOrder.ts
@@ -0,0 +1,302 @@
+/* eslint-disable no-console */
+import Benchmark from 'benchmark';
+import type { DataXY } from 'cheminfo-types';
+import { XSadd } from 'ml-xsadd';
+
+import type { Point } from '../src/types/index.ts';
+import { xGetApproximateSortOrder } from '../src/x/xGetApproximateSortOrder.ts';
+import { xGetSortOrder } from '../src/x/xGetSortOrder.ts';
+
+// The four shapes of sorting the library does, each one with the implementation currently
+// in src copied here, so one process compares them on the same data:
+//   1. permute two parallel arrays  (xySortX, xyUniqueX)
+//   2. read the nth largest value   (xyGetNMaxY, xyFilterTopYValues)
+//   3. materialize the sorted values (getSortedFloat64, xNoiseSanPlot, getSlots)
+//   4. order an array of points    (xyObjectSortX, xyObjectBestPoints)
+
+// The library keeps the comparator below RADIX_ORDER_MIN_LENGTH: the crossover measured on
+// V8 and JavaScriptCore, for all three input shapes, is around 150 values.
+const SIZES = [1_000, 10_000, 100_000, 1_000_000];
+
+/**
+ * A spectrum whose x values are in random order, as xySortX gets them when the points come
+ * from several sources.
+ * @param size - number of points.
+ * @returns unordered x and the y that goes with them.
+ */
+function makeShuffledSpectrum(size: number): DataXY<Float64Array> {
+  const { random } = new XSadd(42);
+  const x = new Float64Array(size);
+  const y = new Float64Array(size);
+  for (let i = 0; i < size; i++) {
+    x[i] = i * 0.01 + 100;
+    y[i] = 10 ** (random() * 5);
+  }
+  for (let i = size - 1; i > 0; i--) {
+    const j = Math.floor(random() * (i + 1));
+    let swap = x[i];
+    x[i] = x[j];
+    x[j] = swap;
+    swap = y[i];
+    y[i] = y[j];
+    y[j] = swap;
+  }
+  return { x, y };
+}
+
+// 1. Permute two parallel arrays -------------------------------------------------------------
+
+/**
+ * xySortX as it is today: one object per point, sorted by a comparator, then scattered back.
+ * @param data - points, x in any order.
+ * @returns the same points, ascending by x.
+ */
+function sortXYByObjects(data: DataXY): DataXY<Float64Array> {
+  const { x, y } = data;
+  const xyObject = Array.from(x, (val, index) => ({
+    x: val,
+    y: y[index],
+  }));
+  xyObject.sort((a, b) => a.x - b.x);
+
+  const response = {
+    x: new Float64Array(x.length),
+    y: new Float64Array(y.length),
+  };
+  for (let i = 0; i < x.length; i++) {
+    response.x[i] = xyObject[i].x;
+    response.y[i] = xyObject[i].y;
+  }
+  return response;
+}
+
+/**
+ * The same permutation read off xGetSortOrder, with no object allocated.
+ * @param data - points, x in any order.
+ * @returns the same points, ascending by x.
+ */
+function sortXYByOrder(data: DataXY): DataXY<Float64Array> {
+  const { x, y } = data;
+  const length = x.length;
+  const order = xGetSortOrder(x);
+  const sortedX = new Float64Array(length);
+  const sortedY = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    const at = order[i];
+    sortedX[i] = x[at];
+    sortedY[i] = y[at];
+  }
+  return { x: sortedX, y: sortedY };
+}
+
+// 2. Read the nth largest value --------------------------------------------------------------
+
+/**
+ * The threshold xyGetNMaxY uses today: a sorted copy of y, reversed, read at n - 1.
+ * @param y - intensities.
+ * @param n - number of points to keep.
+ * @returns the nth largest intensity.
+ */
+function thresholdBySortedCopy(y: Float64Array, n: number): number {
+  const floatY = Float64Array.from(y);
+  floatY.sort();
+  floatY.reverse();
+  return floatY[n - 1];
+}
+
+/**
+ * The same threshold read through the order, leaving y where it is.
+ * @param y - intensities.
+ * @param n - number of points to keep.
+ * @returns the nth largest intensity.
+ */
+function thresholdByOrder(y: Float64Array, n: number): number {
+  const order = xGetSortOrder(y, { descending: true });
+  return y[order[n - 1]];
+}
+
+/**
+ * The same threshold, ranking only by the high word of each intensity.
+ * @param y - intensities.
+ * @param n - number of points to keep.
+ * @returns the nth largest intensity, to 2^-20 relative.
+ */
+function thresholdByApproximateOrder(y: Float64Array, n: number): number {
+  const order = xGetApproximateSortOrder(y, { descending: true });
+  return y[order[n - 1]];
+}
+
+// 3. Materialize the sorted values -----------------------------------------------------------
+
+/**
+ * getSortedFloat64 as it is today.
+ * @param array - values.
+ * @returns the values, ascending.
+ */
+function sortedValuesByTypedSort(array: Float64Array): Float64Array {
+  const sorted = Float64Array.from(array);
+  sorted.sort();
+  return sorted;
+}
+
+/**
+ * The same values gathered through the order.
+ * @param array - values.
+ * @returns the values, ascending.
+ */
+function sortedValuesByOrder(array: Float64Array): Float64Array {
+  const length = array.length;
+  const order = xGetSortOrder(array);
+  const sorted = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    sorted[i] = array[order[i]];
+  }
+  return sorted;
+}
+
+// 4. Order an array of points ----------------------------------------------------------------
+
+/**
+ * xyObjectSortX as it is today, on a copy so the benchmark never re-sorts sorted points.
+ * @param points - points, x in any order.
+ * @returns the points, ascending by x.
+ */
+function sortPointsByComparator(points: Point[]): Point[] {
+  const sorted = points.slice();
+  sorted.sort((a, b) => a.x - b.x);
+  return sorted;
+}
+
+/**
+ * The same points gathered through the order of their x values.
+ * @param points - points, x in any order.
+ * @returns the points, ascending by x.
+ */
+function sortPointsByOrder(points: Point[]): Point[] {
+  const length = points.length;
+  const x = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    x[i] = points[i].x;
+  }
+  const order = xGetSortOrder(x);
+  const sorted = new Array<Point>(length);
+  for (let i = 0; i < length; i++) {
+    sorted[i] = points[order[i]];
+  }
+  // xyObjectSortX orders in place, so the gathered points are written back over the input.
+  const inPlace = points.slice();
+  for (let i = 0; i < length; i++) {
+    inPlace[i] = sorted[i];
+  }
+  return inPlace;
+}
+
+/**
+ * Runs one A/B suite and prints ns per element for each case.
+ * @param title - what is being compared.
+ * @param size - number of elements one call goes through.
+ * @param cases - the implementations to compare, the current one first.
+ */
+function compare(
+  title: string,
+  size: number,
+  cases: Array<[string, () => unknown]>,
+): void {
+  const nanoseconds: Record<string, number> = {};
+  const rmes: Record<string, number> = {};
+  const suite = new Benchmark.Suite();
+  for (const [name, run] of cases) {
+    suite.add(name, run, { minSamples: 30 });
+  }
+  suite
+    .on('cycle', (event: Benchmark.Event) => {
+      const { stats, name } = event.target;
+      if (!stats || !name) return;
+      nanoseconds[name] = (stats.mean * 1e9) / size;
+      rmes[name] = stats.rme;
+    })
+    .run();
+
+  const reference = nanoseconds[cases[0][0]];
+  const line = cases
+    .map(([name]) => {
+      const ns = nanoseconds[name];
+      const speedup =
+        name === cases[0][0] ? '' : ` ${(reference / ns).toFixed(2)}x`;
+      return `${name} ${ns.toFixed(2)} ns/pt (±${rmes[name].toFixed(1)}%)${speedup}`;
+    })
+    .join(' | ');
+  console.log(`  ${title.padEnd(22)} ${line}`);
+}
+
+for (const size of SIZES) {
+  const data = makeShuffledSpectrum(size);
+  const { x, y } = data;
+  const n = Math.max(1, Math.round(size / 100));
+
+  const byObjects = sortXYByObjects(data);
+  const byOrder = sortXYByOrder(data);
+  let samePermutation = true;
+  for (let i = 0; i < size; i++) {
+    if (byObjects.x[i] !== byOrder.x[i] || byObjects.y[i] !== byOrder.y[i]) {
+      samePermutation = false;
+      break;
+    }
+  }
+
+  const exactThreshold = thresholdBySortedCopy(y, n);
+  const orderThreshold = thresholdByOrder(y, n);
+  const approximateThreshold = thresholdByApproximateOrder(y, n);
+
+  const typedSorted = sortedValuesByTypedSort(x);
+  const orderSorted = sortedValuesByOrder(x);
+  let sameValues = true;
+  for (let i = 0; i < size; i++) {
+    if (typedSorted[i] !== orderSorted[i]) {
+      sameValues = false;
+      break;
+    }
+  }
+
+  const points: Point[] = new Array(size);
+  for (let i = 0; i < size; i++) {
+    points[i] = { x: x[i], y: y[i] };
+  }
+  const byComparator = sortPointsByComparator(points);
+  const pointsByOrder = sortPointsByOrder(points);
+  let samePoints = true;
+  for (let i = 0; i < size; i++) {
+    if (byComparator[i] !== pointsByOrder[i]) {
+      samePoints = false;
+      break;
+    }
+  }
+
+  console.log(`\nn = ${size}`);
+  compare('permute x and y', size, [
+    ['objects', () => sortXYByObjects(data)],
+    ['order', () => sortXYByOrder(data)],
+  ]);
+  console.log(`    same points: ${samePermutation}`);
+
+  compare(`${n}th largest y`, size, [
+    ['sorted copy', () => thresholdBySortedCopy(y, n)],
+    ['order', () => thresholdByOrder(y, n)],
+    ['approximate order', () => thresholdByApproximateOrder(y, n)],
+  ]);
+  console.log(
+    `    threshold: sorted copy ${exactThreshold} | order ${orderThreshold} | approximate ${approximateThreshold}`,
+  );
+
+  compare('sorted values', size, [
+    ['typed sort', () => sortedValuesByTypedSort(x)],
+    ['order + gather', () => sortedValuesByOrder(x)],
+  ]);
+  console.log(`    same values: ${sameValues}`);
+
+  compare('order points', size, [
+    ['comparator', () => sortPointsByComparator(points)],
+    ['order', () => sortPointsByOrder(points)],
+  ]);
+  console.log(`    same order: ${samePoints}`);
+}
diff --git a/package-lock.json b/package-lock.json
index 86fec393a..d7fd0314b 100644
--- a/package-lock.json
+++ b/package-lock.json
@@ -10,7 +10,7 @@
       "license": "MIT",
       "dependencies": {
         "binary-search": "^1.3.6",
-        "cheminfo-types": "^1.15.0",
+        "cheminfo-types": "^1.16.0",
         "fft.js": "^4.0.4",
         "is-any-array": "^3.0.0",
         "ml-matrix": "^6.14.0",
@@ -3848,9 +3848,9 @@
       }
     },
     "node_modules/cheminfo-types": {
-      "version": "1.15.0",
-      "resolved": "https://registry.npmjs.org/cheminfo-types/-/cheminfo-types-1.15.0.tgz",
-      "integrity": "sha512-shv45WN2u0yN9EHH1bisNrv+fy4Cw+eLM5lOoriP67mePrwbHZ1kJqg90C8GEU7K1A8gJsicEoVZHcuBbuul/w==",
+      "version": "1.16.0",
+      "resolved": "https://registry.npmjs.org/cheminfo-types/-/cheminfo-types-1.16.0.tgz",
+      "integrity": "sha512-j0MiXFnLRVEZg3vUr7SXm6F3HuDNVE5arfQyDMiy1w2fAcuOs3QRgVsw7aOIACvr0aqbvn85cucwt+yOu64IiQ==",
       "license": "MIT"
     },
     "node_modules/ci-info": {
diff --git a/package.json b/package.json
index c292553a6..4b3398596 100644
--- a/package.json
+++ b/package.json
@@ -18,6 +18,7 @@
     "./xy2": "./lib/xy2/index.js",
     "./xyArray": "./lib/xyArray/index.js",
     "./xyObject": "./lib/xyObject/index.js",
+    "./xys": "./lib/xys/index.js",
     "./zones": "./lib/zones/index.js"
   },
   "files": [
@@ -40,7 +41,7 @@
   },
   "dependencies": {
     "binary-search": "^1.3.6",
-    "cheminfo-types": "^1.15.0",
+    "cheminfo-types": "^1.16.0",
     "fft.js": "^4.0.4",
     "is-any-array": "^3.0.0",
     "ml-matrix": "^6.14.0",
diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap
index c44dcf1ec..8bdabba27 100644
--- a/src/__tests__/__snapshots__/index.test.ts.snap
+++ b/src/__tests__/__snapshots__/index.test.ts.snap
@@ -30,10 +30,13 @@ exports[`existence of exported functions 1`] = `
   "xDistributionStats",
   "xDivide",
   "xDotProduct",
+  "xDoubleTypedArrayLength",
   "xEnsureFloat64",
   "xEqualIntegrationVectorSimilarity",
   "xFindClosestIndex",
+  "xGetApproximateSortOrder",
   "xGetFromToIndex",
+  "xGetSortOrder",
   "xGetTargetIndex",
   "xHilbertTransform",
   "xHistogram",
@@ -138,8 +141,12 @@ exports[`existence of exported functions 1`] = `
   "xreimSortX",
   "xreimZeroFilling",
   "xyArrayAlign",
+  "xyArrayAlignByIntensity",
   "xyArrayAlignToFirst",
+  "xyArrayFilterMinYValue",
+  "xyArrayMaxY",
   "xyArrayMerge",
+  "xyArrayMergeByIntensity",
   "xyArrayWeightedMerge",
   "xyObjectBestPoints",
   "xyObjectCheck",
@@ -154,6 +161,7 @@ exports[`existence of exported functions 1`] = `
   "xyObjectSortX",
   "xyObjectSumY",
   "xyObjectToXY",
+  "xysFilter",
   "zonesNormalize",
   "zonesWithPoints",
   "zonesToXEquallySpaced",
@@ -195,6 +203,7 @@ exports[`existence of exported functions 1`] = `
   "getCombinations",
   "getCombinationsIterator",
   "getRescaler",
+  "isLittleEndian",
   "isPowerOfTwo",
   "nextPowerOfTwo",
   "recursiveResolve",
diff --git a/src/index.ts b/src/index.ts
index 113946a47..d95d37ede 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -12,6 +12,7 @@ export * from './xreim/index.ts';
 
 export * from './xyArray/index.ts';
 export * from './xyObject/index.ts';
+export * from './xys/index.ts';
 
 export * from './zones/index.ts';
 
diff --git a/src/utils/__tests__/isLittleEndian.test.ts b/src/utils/__tests__/isLittleEndian.test.ts
new file mode 100644
index 000000000..34c59cf4e
--- /dev/null
+++ b/src/utils/__tests__/isLittleEndian.test.ts
@@ -0,0 +1,14 @@
+import { expect, test } from 'vitest';
+
+import { isLittleEndian } from '../isLittleEndian.ts';
+
+test('matches the byte order seen through a DataView', () => {
+  const view = new DataView(new ArrayBuffer(2));
+  new Uint16Array(view.buffer)[0] = 1;
+
+  expect(isLittleEndian()).toBe(view.getUint16(0, true) === 1);
+});
+
+test('is stable across calls', () => {
+  expect(isLittleEndian()).toBe(isLittleEndian());
+});
diff --git a/src/utils/index.ts b/src/utils/index.ts
index d7b069db1..e1c100ebb 100644
--- a/src/utils/index.ts
+++ b/src/utils/index.ts
@@ -6,6 +6,7 @@ export { clearFFTCache, setFFTCacheMaxSize } from './fftCache.ts';
 export * from './getCombinations.ts';
 export * from './getCombinationsIterator.ts';
 export * from './getRescaler.ts';
+export * from './isLittleEndian.ts';
 export * from './isPowerOfTwo.ts';
 export * from './nextPowerOfTwo.ts';
 export * from './recursiveResolve.ts';
diff --git a/src/utils/isLittleEndian.ts b/src/utils/isLittleEndian.ts
new file mode 100644
index 000000000..7807f9aa3
--- /dev/null
+++ b/src/utils/isLittleEndian.ts
@@ -0,0 +1,11 @@
+const LITTLE_ENDIAN = new Uint8Array(new Uint16Array([1]).buffer)[0] === 1;
+
+/**
+ * Check whether the host stores multi-byte values least significant byte first.
+ * Needed whenever a typed array is reinterpreted through another view, because
+ * the position of the words of a value then depends on the platform.
+ * @returns true on a little-endian host.
+ */
+export function isLittleEndian(): boolean {
+  return LITTLE_ENDIAN;
+}
diff --git a/src/x/__tests__/xDoubleTypedArrayLength.test.ts b/src/x/__tests__/xDoubleTypedArrayLength.test.ts
new file mode 100644
index 000000000..6e4981505
--- /dev/null
+++ b/src/x/__tests__/xDoubleTypedArrayLength.test.ts
@@ -0,0 +1,33 @@
+import { expect, test } from 'vitest';
+
+import { xDoubleTypedArrayLength } from '../xDoubleTypedArrayLength.ts';
+
+test('keeps the values and zeroes the added positions', () => {
+  expect(
+    xDoubleTypedArrayLength(Float64Array.from([1.5, -2, 3])),
+  ).toStrictEqual(Float64Array.from([1.5, -2, 3, 0, 0, 0]));
+});
+
+test('returns an array of the same kind', () => {
+  expect(xDoubleTypedArrayLength(Uint32Array.from([7, 8]))).toStrictEqual(
+    Uint32Array.from([7, 8, 0, 0]),
+  );
+  expect(xDoubleTypedArrayLength(Int16Array.from([-4]))).toStrictEqual(
+    Int16Array.from([-4, 0]),
+  );
+  expect(xDoubleTypedArrayLength(Uint8Array.from([255]))).toStrictEqual(
+    Uint8Array.from([255, 0]),
+  );
+});
+
+test('an empty array becomes one element long, so doubling progresses', () => {
+  expect(xDoubleTypedArrayLength(new Float64Array(0))).toStrictEqual(
+    new Float64Array(1),
+  );
+});
+
+test('does not return the array it was given', () => {
+  const array = Float64Array.from([1, 2]);
+
+  expect(xDoubleTypedArrayLength(array)).not.toBe(array);
+});
diff --git a/src/x/__tests__/xGetApproximateSortOrder.test.ts b/src/x/__tests__/xGetApproximateSortOrder.test.ts
new file mode 100644
index 000000000..1075c0d00
--- /dev/null
+++ b/src/x/__tests__/xGetApproximateSortOrder.test.ts
@@ -0,0 +1,119 @@
+import { XSadd } from 'ml-xsadd';
+import { expect, test } from 'vitest';
+
+import { xGetApproximateSortOrder } from '../xGetApproximateSortOrder.ts';
+
+/**
+ * Applies an order to the values.
+ * @param values - the values.
+ * @param order - the indices, in order.
+ * @returns the reordered values.
+ */
+function apply(values: number[], order: Uint32Array): number[] {
+  return Array.from(order, (index) => values[index]);
+}
+
+test('ascending order of well separated values', () => {
+  const values = [3, 1, 4, 1, 5, 9, 2, 6];
+
+  expect(apply(values, xGetApproximateSortOrder(values))).toStrictEqual([
+    1, 1, 2, 3, 4, 5, 6, 9,
+  ]);
+});
+
+test('descending order of well separated values', () => {
+  const values = [3, 1, 4, 1, 5];
+
+  expect(
+    apply(values, xGetApproximateSortOrder(values, { descending: true })),
+  ).toStrictEqual([5, 4, 3, 1, 1]);
+});
+
+test('negative values order below positive ones', () => {
+  const values = [1, -1, 0, -3.5, 2.5, -0.5];
+
+  expect(apply(values, xGetApproximateSortOrder(values))).toStrictEqual([
+    -3.5, -1, -0.5, 0, 1, 2.5,
+  ]);
+});
+
+test('equal values keep their input order', () => {
+  const values = [2, 1, 2, 1, 2, 1];
+
+  expect(Array.from(xGetApproximateSortOrder(values))).toStrictEqual([
+    1, 3, 5, 0, 2, 4,
+  ]);
+});
+
+test('empty and single element arrays', () => {
+  expect(xGetApproximateSortOrder([])).toStrictEqual(new Uint32Array(0));
+  expect(Array.from(xGetApproximateSortOrder([42]))).toStrictEqual([0]);
+});
+
+test('values closer than the key resolution keep their input order', () => {
+  // The high word ignores the low 32 bits of the mantissa, so these three share a key.
+  const values = [1, 1 + 2 ** -40, 1 + 2 ** -30];
+
+  expect(
+    Array.from(xGetApproximateSortOrder(values, { descending: true })),
+  ).toStrictEqual([0, 1, 2]);
+});
+
+test('orders to 20 bits of relative precision on the 16 bit digit path', () => {
+  const length = 70000;
+  const values = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    values[i] = Math.abs(Math.sin(i)) * 1e6 + 1;
+  }
+  const ordered = Array.from(
+    xGetApproximateSortOrder(values),
+    (i) => values[i],
+  );
+
+  expect(ordered).toHaveLength(length);
+
+  let worst = 0;
+  for (let i = 1; i < length; i++) {
+    if (ordered[i] < ordered[i - 1]) {
+      worst = Math.max(worst, (ordered[i - 1] - ordered[i]) / ordered[i - 1]);
+    }
+  }
+
+  expect(worst).toBeLessThan(2 ** -20);
+});
+
+test('stays within the key resolution of a native descending sort', () => {
+  const length = 5000;
+  const values = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    values[i] = Math.abs(Math.cos(i)) * 1e4 + 1;
+  }
+  const ordered = Float64Array.from(
+    xGetApproximateSortOrder(values, { descending: true }),
+    (i) => values[i],
+  );
+  const native = values.toSorted().toReversed();
+  let worst = 0;
+  for (let i = 0; i < length; i++) {
+    worst = Math.max(worst, Math.abs(ordered[i] - native[i]) / native[i]);
+  }
+
+  expect(worst).toBeLessThan(2 ** -20);
+});
+
+test('matches a comparator sort on seeded values spanning many magnitudes', () => {
+  const random = new XSadd(7);
+  const values = new Float64Array(10000);
+  for (let i = 0; i < values.length; i++) {
+    values[i] = (random.random() - 0.5) * 10 ** ((random.random() - 0.5) * 40);
+  }
+
+  const expected = Array.from(
+    { length: values.length },
+    (_, index) => index,
+  ).toSorted((first, second) => values[second] - values[first]);
+
+  expect(
+    Array.from(xGetApproximateSortOrder(values, { descending: true })),
+  ).toStrictEqual(expected);
+});
diff --git a/src/x/__tests__/xGetSortOrder.test.ts b/src/x/__tests__/xGetSortOrder.test.ts
new file mode 100644
index 000000000..2bcddf112
--- /dev/null
+++ b/src/x/__tests__/xGetSortOrder.test.ts
@@ -0,0 +1,139 @@
+import { expect, test } from 'vitest';
+
+import { xGetSortOrder } from '../xGetSortOrder.ts';
+
+/**
+ * Applies an order to the values.
+ * @param values - the values.
+ * @param order - the indices, in order.
+ * @returns the reordered values.
+ */
+function apply(values: number[], order: Uint32Array): number[] {
+  return Array.from(order, (index) => values[index]);
+}
+
+test('ascending order of a plain array', () => {
+  const values = [3, 1, 4, 1, 5, 9, 2, 6];
+
+  expect(Array.from(xGetSortOrder(values))).toStrictEqual([
+    1, 3, 6, 0, 2, 4, 7, 5,
+  ]);
+  expect(apply(values, xGetSortOrder(values))).toStrictEqual([
+    1, 1, 2, 3, 4, 5, 6, 9,
+  ]);
+});
+
+test('descending order', () => {
+  const values = [3, 1, 4, 1, 5];
+
+  expect(
+    apply(values, xGetSortOrder(values, { descending: true })),
+  ).toStrictEqual([5, 4, 3, 1, 1]);
+});
+
+test('negative values order below positive ones', () => {
+  const values = [1, -1, 0, -3.5, 2.5, -0.5];
+
+  expect(apply(values, xGetSortOrder(values))).toStrictEqual([
+    -3.5, -1, -0.5, 0, 1, 2.5,
+  ]);
+});
+
+test('equal values keep their input order', () => {
+  const values = [2, 1, 2, 1, 2, 1];
+
+  expect(Array.from(xGetSortOrder(values))).toStrictEqual([1, 3, 5, 0, 2, 4]);
+});
+
+test('equal values keep their input order when descending', () => {
+  const values = [2, 1, 2, 1, 2, 1];
+
+  expect(Array.from(xGetSortOrder(values, { descending: true }))).toStrictEqual(
+    [0, 2, 4, 1, 3, 5],
+  );
+});
+
+test('empty and single element arrays', () => {
+  expect(xGetSortOrder([])).toStrictEqual(new Uint32Array(0));
+  expect(Array.from(xGetSortOrder([42]))).toStrictEqual([0]);
+});
+
+test('accepts a Float64Array without modifying it', () => {
+  const values = Float64Array.from([3, 1, 2]);
+
+  expect(Array.from(xGetSortOrder(values))).toStrictEqual([1, 2, 0]);
+  expect(values).toStrictEqual(Float64Array.from([3, 1, 2]));
+});
+
+test('orders infinities, both zeros and subnormals like the native sort', () => {
+  const values = [
+    3,
+    -Infinity,
+    0,
+    -0,
+    Infinity,
+    -3,
+    Number.MIN_VALUE,
+    -Number.MIN_VALUE,
+    5e-324,
+  ];
+  const native = Float64Array.from(values).toSorted();
+
+  expect(Float64Array.from(apply(values, xGetSortOrder(values)))).toStrictEqual(
+    native,
+  );
+});
+
+test('NaN sorts last, like the native sort', () => {
+  const values = [3, Number.NaN, 1, 2];
+
+  expect(apply(values, xGetSortOrder(values))).toStrictEqual([
+    1,
+    2,
+    3,
+    Number.NaN,
+  ]);
+});
+
+test('matches the native sort on the 8 bit digit path', () => {
+  const values = new Float64Array(5000);
+  for (let i = 0; i < values.length; i++) {
+    values[i] = Math.sin(i) * 1e6;
+  }
+  const sorted = values.toSorted();
+  const order = xGetSortOrder(values);
+
+  expect(order).toHaveLength(5000);
+  expect(Float64Array.from(order, (index) => values[index])).toStrictEqual(
+    sorted,
+  );
+});
+
+test('matches the native sort on the 16 bit digit path', () => {
+  const length = 70000;
+  const values = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    values[i] = Math.sin(i) * 1e6 - 5e5;
+  }
+  const sorted = values.toSorted();
+  const order = xGetSortOrder(values);
+
+  expect(order).toHaveLength(length);
+  expect(Float64Array.from(order, (index) => values[index])).toStrictEqual(
+    sorted,
+  );
+});
+
+test('descending is the exact reverse of ascending when values are distinct', () => {
+  const values = new Float64Array(3000);
+  for (let i = 0; i < values.length; i++) {
+    values[i] = Math.cos(i) * 1234.5;
+  }
+  const ascending = Float64Array.from(xGetSortOrder(values), (i) => values[i]);
+  const descending = Float64Array.from(
+    xGetSortOrder(values, { descending: true }),
+    (i) => values[i],
+  );
+
+  expect(descending).toStrictEqual(ascending.toReversed());
+});
diff --git a/src/x/index.ts b/src/x/index.ts
index d185ae2f8..0e3cdcd34 100644
--- a/src/x/index.ts
+++ b/src/x/index.ts
@@ -20,10 +20,13 @@ export * from './xCumulative.ts';
 export * from './xDistributionStats.ts';
 export * from './xDivide.ts';
 export * from './xDotProduct.ts';
+export * from './xDoubleTypedArrayLength.ts';
 export * from './xEnsureFloat64.ts';
 export * from './xEqualIntegrationVectorSimilarity.ts';
 export * from './xFindClosestIndex.ts';
+export * from './xGetApproximateSortOrder.ts';
 export * from './xGetFromToIndex.ts';
+export * from './xGetSortOrder.ts';
 export * from './xGetTargetIndex.ts';
 export * from './xHilbertTransform.ts';
 export * from './xHistogram.ts';
diff --git a/src/x/utils/radixSortOrder.ts b/src/x/utils/radixSortOrder.ts
new file mode 100644
index 000000000..8ca8a0d73
--- /dev/null
+++ b/src/x/utils/radixSortOrder.ts
@@ -0,0 +1,192 @@
+import type { NumberArray } from 'cheminfo-types';
+
+import { isLittleEndian } from '../../utils/isLittleEndian.ts';
+
+// Internal helpers: not re-exported from `x/index.ts`, so they stay out of the public API.
+
+// The counting pass costs O(radix) whatever the length, so the digit stays narrow until the
+// array is long enough to amortize a wide one.
+const SMALL_DIGIT_BITS = 8;
+const LARGE_DIGIT_BITS = 16;
+const LARGE_ARRAY_LENGTH = 1 << LARGE_DIGIT_BITS;
+
+/**
+ * Length from which ordering through the radix beats a comparator, for a caller that has to
+ * build the keys and read the values back through the order. The measured crossover is around
+ * 150 values on both V8 and JavaScriptCore, whether the values come from a typed array, a
+ * plain array or the x of an array of points.
+ */
+export const RADIX_ORDER_MIN_LENGTH = 1 << SMALL_DIGIT_BITS;
+
+/**
+ * Orders the values by their full 64 bits, as an LSD radix sort.
+ * @param array - the values to order.
+ * @param descending - order from the largest to the smallest.
+ * @returns the indices of the values, in order.
+ */
+export function radixSortOrderExact(
+  array: NumberArray,
+  descending: boolean,
+): Uint32Array {
+  const values = ensureFloat64(array);
+  const length = values.length;
+  let order = new Uint32Array(length);
+  if (length < 2) return order;
+
+  const digitBits = getDigitBits(length);
+  const radix = 1 << digitBits;
+  const digitMask = radix - 1;
+  const words = new Uint32Array(values.buffer, values.byteOffset, length * 2);
+  const flip = descending ? 0xffffffff : 0;
+  const highWord = isLittleEndian() ? 1 : 0;
+  const lowWord = 1 - highWord;
+
+  let keysLow = new Uint32Array(length);
+  let keysHigh = new Uint32Array(length);
+  for (let i = 0; i < length; i++) {
+    const high = words[i * 2 + highWord];
+    const low = words[i * 2 + lowWord];
+    // Monotonic unsigned key of the double: a negative flips every bit, a positive only the
+    // sign bit. Complementing that key reverses the order.
+    if ((high & 0x80000000) === 0) {
+      keysHigh[i] = high ^ 0x80000000 ^ flip;
+      keysLow[i] = low ^ flip;
+    } else {
+      keysHigh[i] = ~high ^ flip;
+      keysLow[i] = ~low ^ flip;
+    }
+    order[i] = i;
+  }
+
+  let nextOrder = new Uint32Array(length);
+  let nextLow = new Uint32Array(length);
+  let nextHigh = new Uint32Array(length);
+  const counts = new Uint32Array(radix);
+
+  for (let shift = 0; shift < 64; shift += digitBits) {
+    const onHigh = shift >= 32;
+    const wordShift = onHigh ? shift - 32 : shift;
+    const keys = onHigh ? keysHigh : keysLow;
+
+    counts.fill(0);
+    for (let i = 0; i < length; i++) {
+      counts[(keys[i] >>> wordShift) & digitMask]++;
+    }
+    // A digit shared by every value would make this pass the identity.
+    if (counts[(keys[0] >>> wordShift) & digitMask] === length) continue;
+    cumulate(counts);
+
+    for (let i = 0; i < length; i++) {
+      const low = keysLow[i];
+      const high = keysHigh[i];
+      const at = counts[((onHigh ? high : low) >>> wordShift) & digitMask]++;
+      nextOrder[at] = order[i];
+      nextLow[at] = low;
+      nextHigh[at] = high;
+    }
+    // A destructuring swap makes JavaScriptCore give up the typed array specialization of
+    // these variables, which costs 8x on the whole sort.
+    let swap = order;
+    order = nextOrder;
+    nextOrder = swap;
+    swap = keysLow;
+    keysLow = nextLow;
+    nextLow = swap;
+    swap = keysHigh;
+    keysHigh = nextHigh;
+    nextHigh = swap;
+  }
+
+  return order;
+}
+
+/**
+ * Orders the values by the high word of each double, as an LSD radix sort. Values differing by
+ * less than their 2^-20 relative resolution share a key and keep their input order.
+ * @param array - the values to order.
+ * @param descending - order from the largest to the smallest.
+ * @returns the indices of the values, in order.
+ */
+export function radixSortOrderApproximate(
+  array: NumberArray,
+  descending: boolean,
+): Uint32Array {
+  const values = ensureFloat64(array);
+  const length = values.length;
+  let order = new Uint32Array(length);
+  if (length < 2) return order;
+
+  const digitBits = getDigitBits(length);
+  const radix = 1 << digitBits;
+  const digitMask = radix - 1;
+  const words = new Uint32Array(values.buffer, values.byteOffset, length * 2);
+  const flip = descending ? 0xffffffff : 0;
+  const highWord = isLittleEndian() ? 1 : 0;
+
+  let keys = new Uint32Array(length);
+  for (let i = 0; i < length; i++) {
+    const high = words[i * 2 + highWord];
+    keys[i] = ((high & 0x80000000) === 0 ? high ^ 0x80000000 : ~high) ^ flip;
+    order[i] = i;
+  }
+
+  let nextOrder = new Uint32Array(length);
+  let nextKeys = new Uint32Array(length);
+  const counts = new Uint32Array(radix);
+
+  for (let shift = 0; shift < 32; shift += digitBits) {
+    counts.fill(0);
+    for (let i = 0; i < length; i++) {
+      counts[(keys[i] >>> shift) & digitMask]++;
+    }
+    if (counts[(keys[0] >>> shift) & digitMask] === length) continue;
+    cumulate(counts);
+
+    for (let i = 0; i < length; i++) {
+      const key = keys[i];
+      const at = counts[(key >>> shift) & digitMask]++;
+      nextOrder[at] = order[i];
+      nextKeys[at] = key;
+    }
+    // See radixSortOrderExact: a destructuring swap is 8x slower on JavaScriptCore.
+    let swap = order;
+    order = nextOrder;
+    nextOrder = swap;
+    swap = keys;
+    keys = nextKeys;
+    nextKeys = swap;
+  }
+
+  return order;
+}
+
+/**
+ * Turns the counts of each digit into the position its first entry goes to.
+ * @param counts - counts per digit, replaced in place by the start offsets.
+ */
+function cumulate(counts: Uint32Array): void {
+  let start = 0;
+  for (let i = 0; i < counts.length; i++) {
+    const count = counts[i];
+    counts[i] = start;
+    start += count;
+  }
+}
+
+/**
+ * Picks the digit width, keeping the radix no wider than the array it scans.
+ * @param length - number of values to order.
+ * @returns the number of bits of one digit.
+ */
+function getDigitBits(length: number): number {
+  return length >= LARGE_ARRAY_LENGTH ? LARGE_DIGIT_BITS : SMALL_DIGIT_BITS;
+}
+
+/**
+ * Returns the values as a Float64Array, without copying one that already is.
+ * @param array - values.
+ * @returns the values, readable as doubles.
+ */
+function ensureFloat64(array: NumberArray): Float64Array {
+  return array instanceof Float64Array ? array : Float64Array.from(array);
+}
diff --git a/src/x/xDoubleTypedArrayLength.ts b/src/x/xDoubleTypedArrayLength.ts
new file mode 100644
index 000000000..3576ff9c6
--- /dev/null
+++ b/src/x/xDoubleTypedArrayLength.ts
@@ -0,0 +1,22 @@
+import type { NumberArray } from 'cheminfo-types';
+
+type TypedNumberArray = Exclude<NumberArray, number[]>;
+
+/**
+ * Copies the values into an array of the same kind of twice the length, the added positions
+ * left at zero. A typed array cannot be resized, so a buffer that fills up has to allocate a
+ * longer one and copy into it; doubling is what keeps the total cost of filling a buffer of
+ * unknown size linear. An empty array becomes one element long, so that doubling it repeatedly
+ * makes progress.
+ * @param array - values to copy.
+ * @returns an array of the same kind, twice as long, holding the same values.
+ */
+export function xDoubleTypedArrayLength<ArrayType extends TypedNumberArray>(
+  array: ArrayType,
+): ArrayType {
+  const Constructor = array.constructor as new (length: number) => ArrayType;
+  const grown = new Constructor(Math.max(1, array.length * 2));
+  grown.set(array);
+
+  return grown;
+}
diff --git a/src/x/xGetApproximateSortOrder.ts b/src/x/xGetApproximateSortOrder.ts
new file mode 100644
index 000000000..49aa665fa
--- /dev/null
+++ b/src/x/xGetApproximateSortOrder.ts
@@ -0,0 +1,35 @@
+import type { NumberArray } from 'cheminfo-types';
+
+import { radixSortOrderApproximate } from './utils/radixSortOrder.ts';
+
+export interface XGetApproximateSortOrderOptions {
+  /**
+   * Order the values from the largest to the smallest.
+   * @default false
+   */
+  descending?: boolean;
+}
+
+/**
+ * Returns the indices that approximately order the values, without moving them.
+ * Only the high word of each double is used as the sort key, which halves the work of
+ * `xGetSortOrder` but orders to 20 bits of relative precision instead of exactly: two values
+ * closer than 2^-20 relative share a key and keep their input order. Reading the result back as
+ * values therefore does not give a monotonic sequence — on random data about one adjacent pair
+ * in fourteen is transposed, by at most 2^-20 relative.
+ * Never feed this order to anything that assumes sorted values, a binary search in particular;
+ * use `xGetSortOrder` there. It is meant for ranking by magnitude, such as taking the n most
+ * intense points of a spectrum, where that resolution is far below the noise.
+ * The sort is stable, so equal values keep their input order.
+ * @param array - values to order.
+ * @param options - options.
+ * @returns the indices of the values, in approximate order.
+ */
+export function xGetApproximateSortOrder(
+  array: NumberArray,
+  options: XGetApproximateSortOrderOptions = {},
+): Uint32Array {
+  const { descending = false } = options;
+
+  return radixSortOrderApproximate(array, descending);
+}
diff --git a/src/x/xGetSortOrder.ts b/src/x/xGetSortOrder.ts
new file mode 100644
index 000000000..a2a2dcc59
--- /dev/null
+++ b/src/x/xGetSortOrder.ts
@@ -0,0 +1,30 @@
+import type { NumberArray } from 'cheminfo-types';
+
+import { radixSortOrderExact } from './utils/radixSortOrder.ts';
+
+export interface XGetSortOrderOptions {
+  /**
+   * Order the values from the largest to the smallest.
+   * @default false
+   */
+  descending?: boolean;
+}
+
+/**
+ * Returns the indices that order the values, without moving them.
+ * The sort is a radix sort over the full 64 bits of each double: it costs a fixed number of
+ * linear passes instead of the n log n comparisons of `Array.prototype.sort`, and orders
+ * exactly, NaN, infinities, subnormals and the two zeros included.
+ * The sort is stable, so equal values keep their input order.
+ * @param array - values to order.
+ * @param options - options.
+ * @returns the indices of the values, in order.
+ */
+export function xGetSortOrder(
+  array: NumberArray,
+  options: XGetSortOrderOptions = {},
+): Uint32Array {
+  const { descending = false } = options;
+
+  return radixSortOrderExact(array, descending);
+}
diff --git a/src/xy/__tests__/xySortX.test.ts b/src/xy/__tests__/xySortX.test.ts
index 6bb690460..618c25313 100644
--- a/src/xy/__tests__/xySortX.test.ts
+++ b/src/xy/__tests__/xySortX.test.ts
@@ -1,3 +1,4 @@
+import { XSadd } from 'ml-xsadd';
 import { expect, test } from 'vitest';
 
 import { xySortX } from '../xySortX.ts';
@@ -60,3 +61,21 @@ test('typed XY arrays', () => {
     y: Float64Array.from([8, 9, 1, 7, 2, 0, 5, 10, 6]),
   });
 });
+
+// The short arrays above are ordered by a comparator, the long one by xGetSortOrder: both
+// have to give the same stable order, ties included.
+test('same stable order whatever the length', () => {
+  for (const length of [100, 1000]) {
+    const { random } = new XSadd(42);
+    const x = Array.from({ length }, () => Math.round(random() * 20));
+    const y = Array.from({ length }, (value, index) => index);
+    const expected = Array.from({ length }, (value, index) => index).toSorted(
+      (a, b) => x[a] - x[b],
+    );
+
+    const result = xySortX({ x, y });
+
+    expect(result.x).toStrictEqual(Float64Array.from(expected, (at) => x[at]));
+    expect(result.y).toStrictEqual(Float64Array.from(expected, (at) => y[at]));
+  }
+});
diff --git a/src/xy/xySortX.ts b/src/xy/xySortX.ts
index baf0c9e29..04eb6948c 100644
--- a/src/xy/xySortX.ts
+++ b/src/xy/xySortX.ts
@@ -1,10 +1,12 @@
 import type { DataXY } from 'cheminfo-types';
 
-import { xIsMonotonic } from '../x/index.ts';
+import { xGetSortOrder, xIsMonotonic } from '../x/index.ts';
+import { RADIX_ORDER_MIN_LENGTH } from '../x/utils/radixSortOrder.ts';
 
 /**
  * This function performs a quick sort of the x array while transforming the y array to preserve the coordinates.
  * @param data - object that contains property x (Array) and y (Array).
+ * @returns the points, ascending by x.
  */
 export function xySortX(data: DataXY): DataXY<Float64Array> {
   const { x, y } = data;
@@ -22,16 +24,27 @@ export function xySortX(data: DataXY): DataXY<Float64Array> {
     };
   }
 
+  const response = {
+    x: new Float64Array(x.length),
+    y: new Float64Array(y.length),
+  };
+
+  if (x.length >= RADIX_ORDER_MIN_LENGTH) {
+    const order = xGetSortOrder(x);
+    for (let i = 0; i < order.length; i++) {
+      const at = order[i];
+      response.x[i] = x[at];
+      response.y[i] = y[at];
+    }
+    return response;
+  }
+
   const xyObject = Array.from(x, (val, index) => ({
     x: val,
     y: y[index],
   }));
   xyObject.sort((a, b) => a.x - b.x);
 
-  const response = {
-    x: new Float64Array(x.length),
-    y: new Float64Array(y.length),
-  };
   for (let i = 0; i < x.length; i++) {
     response.x[i] = xyObject[i].x;
     response.y[i] = xyObject[i].y;
diff --git a/src/xyArray/__tests__/xyArrayAlign.test.ts b/src/xyArray/__tests__/xyArrayAlign.test.ts
index 379a91e38..241988fce 100644
--- a/src/xyArray/__tests__/xyArrayAlign.test.ts
+++ b/src/xyArray/__tests__/xyArrayAlign.test.ts
@@ -29,7 +29,11 @@ test('The y values must be present everywhere', () => {
   const result = xyArrayAlign(data, { delta: 0.15, requiredY: true });
 
   expect(result).toStrictEqual({
-    x: [3.025],
-    ys: [[1], [1], [2]],
+    x: Float64Array.from([3.025]),
+    ys: [
+      Float64Array.from([1]),
+      Float64Array.from([1]),
+      Float64Array.from([2]),
+    ],
   });
 });
diff --git a/src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts b/src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts
new file mode 100644
index 000000000..ea4ce7a29
--- /dev/null
+++ b/src/xyArray/__tests__/xyArrayAlignByIntensity.test.ts
@@ -0,0 +1,146 @@
+import { expect, test } from 'vitest';
+
+import { xyArrayAlign } from '../xyArrayAlign.ts';
+import { xyArrayAlignByIntensity } from '../xyArrayAlignByIntensity.ts';
+
+test('the most intense peak takes the peaks of the other spectra', () => {
+  const data = [
+    { x: [10, 10.4, 10.8], y: [1, 10, 1] },
+    { x: [10.2, 10.6], y: [2, 2] },
+  ];
+  const result = xyArrayAlignByIntensity(data, { delta: 0.25 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([10, 10.4, 10.8]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [1, 10, 1],
+    [0, 4, 0],
+  ]);
+});
+
+test('same length spectra, simple integers', () => {
+  const data = [
+    { x: [1, 2, 3], y: [1, 1, 1] },
+    { x: [0.1, 1.1, 2.1, 3.1, 4.1], y: [1, 1, 1, 1, 1] },
+    { x: [2.9, 3.1, 3.9, 4.9], y: [1, 1, 1, 1] },
+  ];
+  const result = xyArrayAlignByIntensity(data, { delta: 0.15 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([
+    0.1, 1.05, 2.05, 2.95, 3.1, 3.9, 4.1, 4.9,
+  ]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [0, 1, 1, 1, 0, 0, 0, 0],
+    [1, 1, 1, 0, 1, 0, 1, 0],
+    [0, 0, 0, 1, 1, 1, 0, 1],
+  ]);
+});
+
+test('a slot never grows wider than 2 * delta', () => {
+  // Each x is within delta of the previous one, so a chain of close neighbours
+  // ends up in a single slot spanning the whole range.
+  const x = Array.from({ length: 21 }, (_, index) => 100 + index * 0.004);
+  const y = x.map(() => 1);
+  const chained = xyArrayAlign([{ x, y }], { delta: 0.005 });
+
+  expect(chained.x).toHaveLength(1);
+  expect(chained.x[0]).toBeCloseTo(100.04, 10);
+
+  const result = xyArrayAlignByIntensity([{ x, y }], { delta: 0.005 });
+
+  // Each slot holds the peak that opened it and the one 0.004 further, never more.
+  expect(Array.from(result.x)).toBeDeepCloseTo([
+    100.002, 100.01, 100.018, 100.026, 100.034, 100.042, 100.05, 100.058,
+    100.066, 100.074, 100.08,
+  ]);
+  expect(Array.from(result.ys[0])).toStrictEqual([
+    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,
+  ]);
+});
+
+test('delta may be a function of x', () => {
+  const data = [
+    { x: [100, 100.5, 1000], y: [1, 1, 1] },
+    { x: [100.6, 1005], y: [1, 1] },
+  ];
+  const result = xyArrayAlignByIntensity(data, { delta: (x) => x * 0.006 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([100.36666666666666, 1002.5]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [2, 1],
+    [1, 1],
+  ]);
+});
+
+test('no data', () => {
+  expect(xyArrayAlignByIntensity([])).toStrictEqual({
+    x: new Float64Array(0),
+    ys: [],
+  });
+  expect(xyArrayAlignByIntensity([{ x: [], y: [] }])).toStrictEqual({
+    x: new Float64Array(0),
+    ys: [new Float64Array(0)],
+  });
+});
+
+test('x values must be ascending', () => {
+  expect(() =>
+    xyArrayAlignByIntensity([
+      { x: [1, 2], y: [1, 1] },
+      { x: [3, 2], y: [1, 1] },
+    ]),
+  ).toThrow('x values of the spectrum 1 must be ascending');
+});
+
+// The four spectra share the same x, but the big peak sits at 1 for the first two
+// and at 1.5 for the last two.
+const twoCompetingPeaks = [
+  { x: [1, 1.5, 10], y: [100, 1, 1] },
+  { x: [1, 1.5, 10], y: [10, 1, 1] },
+  { x: [1, 1.5, 10], y: [1, 10, 1] },
+  { x: [1, 1.5, 10], y: [1, 100, 1] },
+];
+
+test('competing peaks closer than delta end up in one slot', () => {
+  // The 100 at x = 1 opens the first slot, and delta = 1 reaches 1.5, so every
+  // point of both peaks joins it: 224 of intensity at (112 * 1 + 112 * 1.5) / 224.
+  const result = xyArrayAlignByIntensity(twoCompetingPeaks);
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([1.25, 10]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [101, 1],
+    [11, 1],
+    [11, 1],
+    [101, 1],
+  ]);
+});
+
+test('competing peaks further apart than delta keep their own slot', () => {
+  // 1.5 is now out of reach of the slot opened at 1, so each peak keeps its own
+  // and the spectra come back exactly as they went in.
+  const result = xyArrayAlignByIntensity(twoCompetingPeaks, { delta: 0.3 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([1, 1.5, 10]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [100, 1, 1],
+    [10, 1, 1],
+    [1, 10, 1],
+    [1, 100, 1],
+  ]);
+});
+
+test('the weighted centre is pulled towards the more intense spectrum', () => {
+  // Every peak of the second spectrum sits 0.5 above its counterpart in the first and
+  // within delta, so each pair shares a slot. The centre is the intensity weighted
+  // average, so it lands near the strong spectrum rather than halfway: a peak of 10 at
+  // 1 and one of 1 at 1.5 give (10 * 1 + 1 * 1.5) / 11. The peaks at 10 coincide.
+  const result = xyArrayAlignByIntensity([
+    { x: [1, 5, 10], y: [10, 10, 1] },
+    { x: [1.5, 5.5, 10], y: [1, 1, 1] },
+  ]);
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([11.5 / 11, 55.5 / 11, 10]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [10, 10, 1],
+    [1, 1, 1],
+  ]);
+});
diff --git a/src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts b/src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts
new file mode 100644
index 000000000..920dc6d56
--- /dev/null
+++ b/src/xyArray/__tests__/xyArrayFilterMinYValue.test.ts
@@ -0,0 +1,66 @@
+import { expect, test } from 'vitest';
+
+import { xyArrayFilterMinYValue } from '../xyArrayFilterMinYValue.ts';
+
+const data = [
+  { x: [1, 2, 3], y: [100, 5, 1] },
+  { x: [4, 5], y: [50, 2] },
+];
+
+test('the threshold comes from the largest y of all the spectra', () => {
+  // 5% of 100, so the 2 of the second spectrum goes even though it is 4% of its own max
+  const result = xyArrayFilterMinYValue(data, 0.05);
+
+  expect(result.map((spectrum) => Array.from(spectrum.x))).toStrictEqual([
+    [1, 2],
+    [4],
+  ]);
+  expect(result.map((spectrum) => Array.from(spectrum.y))).toStrictEqual([
+    [100, 5],
+    [50],
+  ]);
+});
+
+test('a spectrum may lose all of its points', () => {
+  const result = xyArrayFilterMinYValue(data, 0.6);
+
+  expect(result.map((spectrum) => Array.from(spectrum.x))).toStrictEqual([
+    [1],
+    [],
+  ]);
+});
+
+test('a spectrum may keep all of its points', () => {
+  const spectra = [
+    { x: [1, 2], y: [100, 90] },
+    { x: [3, 4], y: [50, 1] },
+  ];
+  const result = xyArrayFilterMinYValue(spectra, 0.05);
+
+  expect(result[0].x).toStrictEqual([1, 2]);
+  expect(result[1].x).toStrictEqual([3]);
+});
+
+test('without a value nothing is filtered', () => {
+  expect(xyArrayFilterMinYValue(data)).toBe(data);
+});
+
+test('no data', () => {
+  expect(xyArrayFilterMinYValue([], 0.5)).toStrictEqual([]);
+  expect(xyArrayFilterMinYValue([{ x: [], y: [] }], 0.5)).toStrictEqual([
+    { x: [], y: [] },
+  ]);
+});
+
+test('the filtered spectra are plain arrays', () => {
+  const filtered = xyArrayFilterMinYValue(
+    [{ x: Float64Array.from([1, 2]), y: Float64Array.from([100, 1]) }],
+    0.5,
+  );
+
+  // The overload gives number[], so the arrays can be pushed to.
+  filtered[0].x.push(3);
+
+  expect(filtered[0].x).toStrictEqual([1, 3]);
+  expect(filtered[0].y).toStrictEqual([100]);
+});
diff --git a/src/xyArray/__tests__/xyArrayMaxY.test.ts b/src/xyArray/__tests__/xyArrayMaxY.test.ts
new file mode 100644
index 000000000..973d765c3
--- /dev/null
+++ b/src/xyArray/__tests__/xyArrayMaxY.test.ts
@@ -0,0 +1,54 @@
+import { expect, test } from 'vitest';
+
+import { xyArrayMaxY } from '../xyArrayMaxY.ts';
+
+test('the largest y of all the spectra', () => {
+  expect(
+    xyArrayMaxY([
+      { x: [1, 2], y: [10, 50] },
+      { x: [3], y: [30] },
+    ]),
+  ).toBe(50);
+});
+
+test('empty spectra are skipped', () => {
+  expect(
+    xyArrayMaxY([
+      { x: [], y: [] },
+      { x: [1], y: [7] },
+      { x: [], y: [] },
+    ]),
+  ).toBe(7);
+});
+
+test('every y may be negative', () => {
+  expect(
+    xyArrayMaxY([
+      { x: [1], y: [-5] },
+      { x: [2], y: [-2] },
+    ]),
+  ).toBe(-2);
+});
+
+test('typed arrays', () => {
+  expect(
+    xyArrayMaxY([
+      { x: Float64Array.from([1]), y: Float64Array.from([3]) },
+      { x: Float64Array.from([2]), y: Float64Array.from([9]) },
+    ]),
+  ).toBe(9);
+});
+
+test('no point at all', () => {
+  expect(() => xyArrayMaxY([])).toThrow('can not process empty arrays');
+  expect(() => xyArrayMaxY([{ x: [], y: [] }])).toThrow(
+    'can not process empty arrays',
+  );
+});
+
+test('the y must contain numbers', () => {
+  expect(() =>
+    // @ts-expect-error the y are checked at runtime
+    xyArrayMaxY([{ x: [1], y: ['a'] }]),
+  ).toThrow('input must contain numbers');
+});
diff --git a/src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts b/src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts
new file mode 100644
index 000000000..6c798582e
--- /dev/null
+++ b/src/xyArray/__tests__/xyArrayMergeByIntensity.test.ts
@@ -0,0 +1,116 @@
+import type { DataXY } from 'cheminfo-types';
+import { XSadd } from 'ml-xsadd';
+import { expect, test } from 'vitest';
+
+import { xyArrayMergeByIntensity } from '../xyArrayMergeByIntensity.ts';
+
+test('the x of a slot is the intensity weighted average of its peaks', () => {
+  const data = [
+    { x: [100.002, 200.01], y: [10, 20] },
+    { x: [100.001, 200.02], y: [11, 21] },
+  ];
+  const result = xyArrayMergeByIntensity(data, { delta: 0.02 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([
+    (100.002 * 10 + 100.001 * 11) / 21,
+    (200.01 * 20 + 200.02 * 21) / 41,
+  ]);
+  expect(Array.from(result.y)).toStrictEqual([21, 41]);
+  expect(Array.from(result.number)).toStrictEqual([2, 2]);
+  expect(Array.from(result.from)).toStrictEqual([100.001, 200.01]);
+  expect(Array.from(result.to)).toStrictEqual([100.002, 200.02]);
+});
+
+test('the strongest peak takes its neighbours before the weaker ones do', () => {
+  // 5 is closer to 4 than to 10, but 10 is the most intense peak of the set and
+  // takes it first, which leaves 4 alone.
+  const data = [{ x: [4, 5, 10], y: [3, 1, 100] }];
+  const result = xyArrayMergeByIntensity(data, { delta: 5.5 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([4, (5 * 1 + 10 * 100) / 101]);
+  expect(Array.from(result.y)).toStrictEqual([3, 101]);
+  expect(Array.from(result.number)).toStrictEqual([1, 2]);
+});
+
+test('a slot may hold several peaks of the same spectrum', () => {
+  const data = [{ x: [1, 1.1, 1.2], y: [1, 5, 1] }];
+  const result = xyArrayMergeByIntensity(data, { delta: 0.2 });
+
+  expect(Array.from(result.x)).toBeDeepCloseTo([1.1]);
+  expect(Array.from(result.y)).toStrictEqual([7]);
+  expect(Array.from(result.number)).toStrictEqual([3]);
+});
+
+test('no data', () => {
+  const result = xyArrayMergeByIntensity([]);
+
+  expect(Array.from(result.x)).toStrictEqual([]);
+  expect(Array.from(result.y)).toStrictEqual([]);
+});
+
+test('many drifting spectra keep the total intensity and stay within delta', () => {
+  const delta = 0.005;
+  const random = new XSadd(42);
+  const trueMasses = Array.from(
+    { length: 200 },
+    (_, index) => 100 + index * 0.5,
+  );
+  const data: DataXY[] = [];
+  let totalIntensity = 0;
+  for (let spectrum = 0; spectrum < 500; spectrum++) {
+    const x: number[] = [];
+    const y: number[] = [];
+    for (const mass of trueMasses) {
+      if (random.random() < 0.3) continue;
+      x.push(mass + (random.random() - 0.5) * 0.004);
+      const intensity = 1 + random.random() * 1000;
+      y.push(intensity);
+      totalIntensity += intensity;
+    }
+    data.push({ x, y });
+  }
+
+  const result = xyArrayMergeByIntensity(data, { delta });
+
+  expect(result.x).toHaveLength(200);
+
+  let sum = 0;
+  let numberOfPeaks = 0;
+  for (let i = 0; i < result.x.length; i++) {
+    sum += result.y[i];
+    numberOfPeaks += result.number[i];
+
+    expect(result.to[i] - result.from[i]).toBeLessThanOrEqual(2 * delta);
+    expect(result.x[i]).toBeCloseTo(trueMasses[i], 2);
+  }
+
+  for (let i = 1; i < result.x.length; i++) {
+    expect(result.x[i]).toBeGreaterThan(result.x[i - 1]);
+  }
+
+  expect(sum).toBeCloseTo(totalIntensity, 6);
+  expect(numberOfPeaks).toBe(data.reduce((total, s) => total + s.x.length, 0));
+});
+
+test('slots with no intensity keep the x of the peak that opened them', () => {
+  const result = xyArrayMergeByIntensity([{ x: [1, 1.05, 5], y: [0, 0, 0] }], {
+    delta: 0.1,
+  });
+
+  expect(Array.from(result.x)).toStrictEqual([1, 5]);
+  expect(Array.from(result.y)).toStrictEqual([0, 0]);
+  expect(Array.from(result.number)).toStrictEqual([2, 1]);
+});
+
+test('more slots than the initially allocated ones', () => {
+  const length = 5000;
+  const x = Array.from({ length }, (_, index) => index * 10);
+  const y = Array.from({ length }, (_, index) => (index % 7) + 1);
+
+  const result = xyArrayMergeByIntensity([{ x, y }], { delta: 1 });
+
+  expect(result.x).toHaveLength(length);
+  expect(Array.from(result.x)).toStrictEqual(x);
+  expect(Array.from(result.y)).toStrictEqual(y);
+  expect(Array.from(result.number)).toStrictEqual(x.map(() => 1));
+});
diff --git a/src/xyArray/index.ts b/src/xyArray/index.ts
index a084fb3dd..aea413682 100644
--- a/src/xyArray/index.ts
+++ b/src/xyArray/index.ts
@@ -1,4 +1,8 @@
 export * from './xyArrayAlign.ts';
+export * from './xyArrayAlignByIntensity.ts';
 export * from './xyArrayAlignToFirst.ts';
+export * from './xyArrayFilterMinYValue.ts';
+export * from './xyArrayMaxY.ts';
 export * from './xyArrayMerge.ts';
+export * from './xyArrayMergeByIntensity.ts';
 export * from './xyArrayWeightedMerge.ts';
diff --git a/src/xyArray/utils/__tests__/mergeSortedXY.test.ts b/src/xyArray/utils/__tests__/mergeSortedXY.test.ts
new file mode 100644
index 000000000..2685a9a73
--- /dev/null
+++ b/src/xyArray/utils/__tests__/mergeSortedXY.test.ts
@@ -0,0 +1,126 @@
+import { expect, test } from 'vitest';
+
+import { mergeSortedXY } from '../mergeSortedXY.ts';
+
+test('merges the points of every spectrum ascending by x', () => {
+  const merged = mergeSortedXY([
+    { x: [1, 4, 7], y: [10, 40, 70] },
+    { x: [], y: [] },
+    { x: [2, 4.5, 9], y: [20, 45, 90] },
+    { x: [0], y: [0] },
+  ]);
+
+  expect(Array.from(merged.x)).toStrictEqual([0, 1, 2, 4, 4.5, 7, 9]);
+  expect(Array.from(merged.y)).toStrictEqual([0, 10, 20, 40, 45, 70, 90]);
+  expect(Array.from(merged.spectrumIndex)).toStrictEqual([3, 0, 2, 0, 2, 0, 2]);
+});
+
+test('points sharing an x are all kept', () => {
+  const merged = mergeSortedXY([
+    { x: [4], y: [40] },
+    { x: [4], y: [41] },
+  ]);
+
+  expect(Array.from(merged.x)).toStrictEqual([4, 4]);
+  expect(Array.from(merged.y).toSorted((a, b) => a - b)).toStrictEqual([
+    40, 41,
+  ]);
+  expect(
+    Array.from(merged.spectrumIndex).toSorted((a, b) => a - b),
+  ).toStrictEqual([0, 1]);
+});
+
+test('no data', () => {
+  const merged = mergeSortedXY([]);
+
+  expect(Array.from(merged.x)).toStrictEqual([]);
+  expect(Array.from(merged.spectrumIndex)).toStrictEqual([]);
+});
+
+test('x values must be ascending', () => {
+  expect(() => mergeSortedXY([{ x: [3, 1], y: [1, 1] }])).toThrow(
+    'x values of the spectrum 0 must be ascending',
+  );
+});
+
+// Thirty spectra or more are merged by ordering the concatenated points instead of
+// by the heap, so the cases below have to reach that many to exercise it.
+const MANY_SPECTRA = 30;
+
+/**
+ * Builds spectra whose x values interleave, one point of each per unit.
+ * @param numberOfSpectra - how many spectra.
+ * @returns the spectra.
+ */
+function interleavedSpectra(numberOfSpectra: number) {
+  return Array.from({ length: numberOfSpectra }, (_, spectrum) => ({
+    x: [spectrum / numberOfSpectra, 1 + spectrum / numberOfSpectra],
+    y: [spectrum, 100 + spectrum],
+  }));
+}
+
+test('many spectra are merged ascending by x', () => {
+  const merged = mergeSortedXY(interleavedSpectra(MANY_SPECTRA));
+
+  expect(merged.x).toHaveLength(2 * MANY_SPECTRA);
+  expect(Array.from(merged.x)).toStrictEqual(
+    Array.from(merged.x).toSorted((a, b) => a - b),
+  );
+  // Each spectrum contributes its low point, in order, then its high point.
+  expect(Array.from(merged.spectrumIndex)).toStrictEqual([
+    ...Array.from({ length: MANY_SPECTRA }, (_, i) => i),
+    ...Array.from({ length: MANY_SPECTRA }, (_, i) => i),
+  ]);
+  expect(Array.from(merged.y)).toStrictEqual([
+    ...Array.from({ length: MANY_SPECTRA }, (_, i) => i),
+    ...Array.from({ length: MANY_SPECTRA }, (_, i) => 100 + i),
+  ]);
+});
+
+test('both strategies agree on the same points', () => {
+  // Empty spectra contribute no point and keep the indices of the real ones, so
+  // padding to the threshold sends the very same data through the other strategy.
+  const spectra = interleavedSpectra(10);
+  const padded = [
+    ...spectra,
+    ...Array.from({ length: MANY_SPECTRA - spectra.length }, () => ({
+      x: [] as number[],
+      y: [] as number[],
+    })),
+  ];
+
+  expect(mergeSortedXY(padded)).toStrictEqual(mergeSortedXY(spectra));
+});
+
+test('many spectra still reject descending x values', () => {
+  const spectra = interleavedSpectra(MANY_SPECTRA);
+  spectra[7].x = [5, 2];
+
+  expect(() => mergeSortedXY(spectra)).toThrow(
+    'x values of the spectrum 7 must be ascending',
+  );
+});
+
+test('many spectra keep every point sharing an x', () => {
+  const spectra = Array.from({ length: MANY_SPECTRA }, (_, spectrum) => ({
+    x: [4],
+    y: [spectrum],
+  }));
+  const merged = mergeSortedXY(spectra);
+
+  expect(Array.from(merged.x)).toStrictEqual(
+    Array.from({ length: MANY_SPECTRA }, () => 4),
+  );
+  expect(Array.from(merged.y).toSorted((a, b) => a - b)).toStrictEqual(
+    Array.from({ length: MANY_SPECTRA }, (_, i) => i),
+  );
+});
+
+test('many spectra including empty ones', () => {
+  const spectra = interleavedSpectra(MANY_SPECTRA);
+  spectra[3] = { x: [], y: [] };
+  const merged = mergeSortedXY(spectra);
+
+  expect(merged.x).toHaveLength(2 * MANY_SPECTRA - 2);
+  expect(Array.from(merged.spectrumIndex)).not.toContain(3);
+});
diff --git a/src/xyArray/utils/getIntensitySlots.ts b/src/xyArray/utils/getIntensitySlots.ts
new file mode 100644
index 000000000..166adbadc
--- /dev/null
+++ b/src/xyArray/utils/getIntensitySlots.ts
@@ -0,0 +1,183 @@
+import { xDoubleTypedArrayLength } from '../../x/xDoubleTypedArrayLength.ts';
+import { xGetApproximateSortOrder } from '../../x/xGetApproximateSortOrder.ts';
+
+import type { MergedXY } from './mergeSortedXY.ts';
+import { sortIntensitySlots } from './sortIntensitySlots.ts';
+
+const FREE = 0xffffffff;
+
+export interface GetIntensitySlotsOptions {
+  /**
+   * A point joins a slot when its x is within `delta` of the x of the most intense
+   * point of that slot. It may also be a function that allows to change `delta`
+   * depending on the x value.
+   * @default 1
+   */
+  delta?: ((arg: number) => number) | number;
+}
+
+export interface IntensitySlots {
+  /**
+   * The intensity weighted average of the x of the points of each slot, ascending.
+   */
+  x: Float64Array;
+  /**
+   * The sum of the y of the points of each slot.
+   */
+  y: Float64Array;
+  /**
+   * The smallest x of each slot.
+   */
+  from: Float64Array;
+  /**
+   * The largest x of each slot.
+   */
+  to: Float64Array;
+  /**
+   * The number of points of each slot.
+   */
+  number: Uint32Array;
+  /**
+   * For each merged point, the index of the slot it belongs to.
+   */
+  slotIndex: Uint32Array;
+}
+
+/**
+ * Groups points around the most intense ones.
+ * The most intense point of the whole set opens the first slot and takes every
+ * point within `delta` of it, the most intense of the remaining points opens the
+ * next one, and so on until every point belongs to exactly one slot. Unlike a
+ * chain of close neighbours, a slot can therefore never grow wider than `2 * delta`.
+ * @param merged - all the points of all the spectra, ascending by x.
+ * @param options - options.
+ * @returns the slots, ascending by x.
+ */
+export function getIntensitySlots(
+  merged: MergedXY,
+  options: GetIntensitySlotsOptions = {},
+): IntensitySlots {
+  const { delta = 1 } = options;
+  const getDelta = typeof delta === 'function' ? delta : () => delta;
+  const { x, y } = merged;
+  const length = x.length;
+
+  const slotIndex = new Uint32Array(length).fill(FREE);
+  // Union-find over the still free positions, so that claiming a point is done
+  // once and the scan of a window never walks over the points already taken.
+  const nextFree = new Uint32Array(length + 1);
+  const previousFree = new Uint32Array(length + 1);
+  for (let i = 0; i <= length; i++) {
+    nextFree[i] = i;
+    previousFree[i] = i;
+  }
+
+  let capacity = 1024;
+  let sumY = new Float64Array(capacity);
+  let sumXY = new Float64Array(capacity);
+  let seeds = new Float64Array(capacity);
+  let from = new Float64Array(capacity);
+  let to = new Float64Array(capacity);
+  let number = new Uint32Array(capacity);
+  let numberOfSlots = 0;
+
+  const order = xGetApproximateSortOrder(y, { descending: true });
+  for (let i = 0; i < length; i++) {
+    const seed = order[i];
+    if (slotIndex[seed] !== FREE) continue;
+
+    if (numberOfSlots === capacity) {
+      capacity *= 2;
+      sumY = xDoubleTypedArrayLength(sumY);
+      sumXY = xDoubleTypedArrayLength(sumXY);
+      seeds = xDoubleTypedArrayLength(seeds);
+      from = xDoubleTypedArrayLength(from);
+      to = xDoubleTypedArrayLength(to);
+      number = xDoubleTypedArrayLength(number);
+    }
+
+    const slot = numberOfSlots++;
+    const center = x[seed];
+    const currentDelta = getDelta(center);
+    const lowerLimit = center - currentDelta;
+    const upperLimit = center + currentDelta;
+
+    let slotSumY = 0;
+    let slotSumXY = 0;
+    let slotNumber = 0;
+    let slotFrom = center;
+    let slotTo = center;
+
+    let position = seed;
+    while (position < length && x[position] <= upperLimit) {
+      slotIndex[position] = slot;
+      nextFree[position] = position + 1;
+      previousFree[position + 1] = position;
+      slotSumY += y[position];
+      slotSumXY += x[position] * y[position];
+      slotNumber++;
+      slotTo = x[position];
+      position = findNext(nextFree, position + 1);
+    }
+
+    position = findPrevious(previousFree, seed);
+    while (position > 0 && x[position - 1] >= lowerLimit) {
+      const claimed = position - 1;
+      slotIndex[claimed] = slot;
+      nextFree[claimed] = claimed + 1;
+      previousFree[claimed + 1] = claimed;
+      slotSumY += y[claimed];
+      slotSumXY += x[claimed] * y[claimed];
+      slotNumber++;
+      slotFrom = x[claimed];
+      position = findPrevious(previousFree, claimed);
+    }
+
+    sumY[slot] = slotSumY;
+    sumXY[slot] = slotSumXY;
+    seeds[slot] = center;
+    from[slot] = slotFrom;
+    to[slot] = slotTo;
+    number[slot] = slotNumber;
+  }
+
+  return sortIntensitySlots({
+    numberOfSlots,
+    sumY,
+    sumXY,
+    seeds,
+    from,
+    to,
+    number,
+    slotIndex,
+  });
+}
+
+/**
+ * Follows the union-find to the first free position at or after `position`.
+ * @param nextFree - union-find of the free positions.
+ * @param position - where to start looking.
+ * @returns the first free position, at most the length of the data.
+ */
+function findNext(nextFree: Uint32Array, position: number): number {
+  while (nextFree[position] !== position) {
+    nextFree[position] = nextFree[nextFree[position]];
+    position = nextFree[position];
+  }
+  return position;
+}
+
+/**
+ * Follows the union-find to the first boundary at or before `position` whose
+ * preceding point is still free.
+ * @param previousFree - union-find of the free positions.
+ * @param position - where to start looking.
+ * @returns the boundary, 0 when no free point remains on the left.
+ */
+function findPrevious(previousFree: Uint32Array, position: number): number {
+  while (previousFree[position] !== position) {
+    previousFree[position] = previousFree[previousFree[position]];
+    position = previousFree[position];
+  }
+  return position;
+}
diff --git a/src/xyArray/utils/mergeSortedXY.ts b/src/xyArray/utils/mergeSortedXY.ts
new file mode 100644
index 000000000..890897eeb
--- /dev/null
+++ b/src/xyArray/utils/mergeSortedXY.ts
@@ -0,0 +1,177 @@
+import type { DataXY } from 'cheminfo-types';
+
+import { xGetSortOrder } from '../../x/xGetSortOrder.ts';
+
+// Below this many spectra the heap wins, above it the sort does. The heap costs
+// O(log k) per point while ordering the concatenated points is independent of k,
+// and the two meet around thirty spectra.
+const MANY_SPECTRA = 30;
+
+export interface MergedXY {
+  /**
+   * All the x values of all the spectra, ascending.
+   */
+  x: Float64Array;
+  /**
+   * The y value matching each x.
+   */
+  y: Float64Array;
+  /**
+   * The index, in `data`, of the spectrum each point comes from.
+   */
+  spectrumIndex: Uint32Array;
+}
+
+/**
+ * Merges the points of many spectra into one array of points ascending by x.
+ * Each spectrum must already have ascending x values.
+ * @param data - spectra, each one with ascending x values.
+ * @returns every point of every spectrum, ascending by x.
+ */
+export function mergeSortedXY(data: DataXY[]): MergedXY {
+  return data.length < MANY_SPECTRA ? mergeByHeap(data) : mergeBySorting(data);
+}
+
+/**
+ * Merges the spectra by repeatedly taking the smallest x of a binary heap of the
+ * spectra, which reads each spectrum in order and costs O(log k) per point.
+ * @param data - spectra, each one with ascending x values.
+ * @returns every point of every spectrum, ascending by x.
+ */
+function mergeByHeap(data: DataXY[]): MergedXY {
+  const total = getTotalLength(data);
+
+  const x = new Float64Array(total);
+  const y = new Float64Array(total);
+  const spectrumIndex = new Uint32Array(total);
+
+  // Binary heap of the spectra, keyed by the x value each one is currently on.
+  const ids = new Uint32Array(data.length);
+  const keys = new Float64Array(data.length);
+  const positions = new Uint32Array(data.length);
+  let size = 0;
+  for (let i = 0; i < data.length; i++) {
+    const spectrum = data[i];
+    if (spectrum.x.length === 0) continue;
+    const key = spectrum.x[0];
+    let child = size++;
+    while (child > 0) {
+      const parent = (child - 1) >> 1;
+      if (keys[parent] <= key) break;
+      keys[child] = keys[parent];
+      ids[child] = ids[parent];
+      child = parent;
+    }
+    keys[child] = key;
+    ids[child] = i;
+  }
+
+  let at = 0;
+  while (size > 0) {
+    const id = ids[0];
+    const spectrum = data[id];
+    const position = positions[id];
+    const currentX = keys[0];
+    x[at] = currentX;
+    y[at] = spectrum.y[position];
+    spectrumIndex[at] = id;
+    at++;
+
+    const next = position + 1;
+    positions[id] = next;
+
+    let key: number;
+    let newId: number;
+    if (next < spectrum.x.length) {
+      key = spectrum.x[next];
+      if (key < currentX) {
+        throw new Error(`x values of the spectrum ${id} must be ascending`);
+      }
+      newId = id;
+    } else {
+      size--;
+      if (size === 0) break;
+      key = keys[size];
+      newId = ids[size];
+    }
+
+    let parent = 0;
+    for (;;) {
+      const left = 2 * parent + 1;
+      if (left >= size) break;
+      const right = left + 1;
+      const child = right < size && keys[right] < keys[left] ? right : left;
+      if (keys[child] >= key) break;
+      keys[parent] = keys[child];
+      ids[parent] = ids[child];
+      parent = child;
+    }
+    keys[parent] = key;
+    ids[parent] = newId;
+  }
+
+  return { x, y, spectrumIndex };
+}
+
+/**
+ * Merges the spectra by concatenating them and ordering the result by x, which
+ * costs the same whatever the number of spectra.
+ * @param data - spectra, each one with ascending x values.
+ * @returns every point of every spectrum, ascending by x.
+ */
+function mergeBySorting(data: DataXY[]): MergedXY {
+  const total = getTotalLength(data);
+
+  const concatenatedX = new Float64Array(total);
+  const concatenatedY = new Float64Array(total);
+  const concatenatedIndex = new Uint32Array(total);
+  let at = 0;
+  for (let i = 0; i < data.length; i++) {
+    const spectrum = data[i];
+    checkAscending(spectrum.x, i);
+    concatenatedX.set(spectrum.x, at);
+    concatenatedY.set(spectrum.y, at);
+    concatenatedIndex.fill(i, at, at + spectrum.x.length);
+    at += spectrum.x.length;
+  }
+
+  const order = xGetSortOrder(concatenatedX);
+  const x = new Float64Array(total);
+  const y = new Float64Array(total);
+  const spectrumIndex = new Uint32Array(total);
+  for (let i = 0; i < total; i++) {
+    const from = order[i];
+    x[i] = concatenatedX[from];
+    y[i] = concatenatedY[from];
+    spectrumIndex[i] = concatenatedIndex[from];
+  }
+
+  return { x, y, spectrumIndex };
+}
+
+/**
+ * Counts the points of every spectrum.
+ * @param data - spectra.
+ * @returns the number of points of all of them.
+ */
+function getTotalLength(data: DataXY[]): number {
+  let total = 0;
+  for (const spectrum of data) {
+    total += spectrum.x.length;
+  }
+  return total;
+}
+
+/**
+ * Throws when the x values of a spectrum are not ascending, which the heap
+ * notices on its own but ordering the points would silently accept.
+ * @param x - the x values of one spectrum.
+ * @param id - the index of the spectrum, for the message.
+ */
+function checkAscending(x: DataXY['x'], id: number): void {
+  for (let i = 1; i < x.length; i++) {
+    if (x[i] < x[i - 1]) {
+      throw new Error(`x values of the spectrum ${id} must be ascending`);
+    }
+  }
+}
diff --git a/src/xyArray/utils/sortIntensitySlots.ts b/src/xyArray/utils/sortIntensitySlots.ts
new file mode 100644
index 000000000..84aaec432
--- /dev/null
+++ b/src/xyArray/utils/sortIntensitySlots.ts
@@ -0,0 +1,54 @@
+import { xGetSortOrder } from '../../x/xGetSortOrder.ts';
+
+import type { IntensitySlots } from './getIntensitySlots.ts';
+
+export interface RawIntensitySlots {
+  numberOfSlots: number;
+  sumY: Float64Array;
+  sumXY: Float64Array;
+  seeds: Float64Array;
+  from: Float64Array;
+  to: Float64Array;
+  number: Uint32Array;
+  slotIndex: Uint32Array;
+}
+
+/**
+ * Turns the slots built by decreasing intensity into slots ascending by x.
+ * @param raw - the slots, by decreasing intensity of the point that opened them.
+ * @returns the slots, ascending by x, and the slot of each merged point.
+ */
+export function sortIntensitySlots(raw: RawIntensitySlots): IntensitySlots {
+  const { numberOfSlots, sumY, sumXY, seeds, slotIndex } = raw;
+
+  const centers = new Float64Array(numberOfSlots);
+  for (let i = 0; i < numberOfSlots; i++) {
+    // A slot with no intensity at all has no weighted average, it stays where
+    // the point that opened it was.
+    centers[i] = sumY[i] === 0 ? seeds[i] : sumXY[i] / sumY[i];
+  }
+
+  const bySlot = xGetSortOrder(centers);
+
+  const x = new Float64Array(numberOfSlots);
+  const y = new Float64Array(numberOfSlots);
+  const from = new Float64Array(numberOfSlots);
+  const to = new Float64Array(numberOfSlots);
+  const number = new Uint32Array(numberOfSlots);
+  const rankOfSlot = new Uint32Array(numberOfSlots);
+  for (let rank = 0; rank < numberOfSlots; rank++) {
+    const slot = bySlot[rank];
+    rankOfSlot[slot] = rank;
+    x[rank] = centers[slot];
+    y[rank] = sumY[slot];
+    from[rank] = raw.from[slot];
+    to[rank] = raw.to[slot];
+    number[rank] = raw.number[slot];
+  }
+
+  for (let i = 0; i < slotIndex.length; i++) {
+    slotIndex[i] = rankOfSlot[slotIndex[i]];
+  }
+
+  return { x, y, from, to, number, slotIndex };
+}
diff --git a/src/xyArray/xyArrayAlign.ts b/src/xyArray/xyArrayAlign.ts
index 64dbc1a74..e545c844a 100644
--- a/src/xyArray/xyArrayAlign.ts
+++ b/src/xyArray/xyArrayAlign.ts
@@ -1,6 +1,7 @@
 import type { DataXY, NumberArray } from 'cheminfo-types';
 
 import { xyJoinX } from '../xy/index.ts';
+import { xysFilter } from '../xys/index.ts';
 
 import { getSlots } from './utils/getSlots.ts';
 
@@ -53,21 +54,9 @@ export function xyArrayAlign(
     }
   }
 
-  if (requiredY) return filterRequiredY(x, ys);
+  if (requiredY) {
+    return xysFilter({ x, ys }, { requiredY: true });
+  }
 
   return { x, ys };
 }
-
-function filterRequiredY(x: Float64Array, ys: Float64Array[]) {
-  const newX: number[] = [];
-  const newYs: number[][] = Array.from(ys, () => []);
-  for (let i = 0; i < x.length; i++) {
-    if (ys.every((y) => y[i] !== 0)) {
-      newX.push(x[i]);
-      for (let j = 0; j < ys.length; j++) {
-        newYs[j].push(ys[j][i]);
-      }
-    }
-  }
-  return { x: newX, ys: newYs };
-}
diff --git a/src/xyArray/xyArrayAlignByIntensity.ts b/src/xyArray/xyArrayAlignByIntensity.ts
new file mode 100644
index 000000000..beb1b2d04
--- /dev/null
+++ b/src/xyArray/xyArrayAlignByIntensity.ts
@@ -0,0 +1,56 @@
+import type { DataXY, DataXYs } from 'cheminfo-types';
+
+import { getIntensitySlots } from './utils/getIntensitySlots.ts';
+import { mergeSortedXY } from './utils/mergeSortedXY.ts';
+
+export interface XYArrayAlignByIntensityOptions {
+  /**
+   * A peak joins a slot when its x is within `delta` of the x of the most intense
+   * peak of that slot, so a slot is never wider than `2 * delta`. It may also be a
+   * function that allows to change `delta` depending on the x value.
+   * @default 1
+   */
+  delta?: ((arg: number) => number) | number;
+}
+
+/**
+ * Aligns spectra on the most intense peaks, which is meant for centroided data
+ * whose x values drift slightly from one spectrum to the next.
+ *
+ * The most intense peak of the whole set defines the first common x, and takes
+ * every peak within `delta` of it. The most intense of the remaining peaks defines
+ * the next one, and so on until every peak has been placed. The reported x of a
+ * slot is the intensity weighted average of the peaks it holds, and its y is their
+ * sum, so no intensity is created or lost. Because a slot only ever reaches
+ * `delta` away from the peak that opened it, dense data cannot chain into a single
+ * wide slot the way `xyArrayAlign` does.
+ *
+ * The x values of each spectrum must be ascending.
+ *
+ * The result is a dense matrix: it holds `data.length * x.length` numbers, which
+ * is out of reach for tens of thousands of spectra. Use
+ * `xyArrayMergeByIntensity` to get the common peaks alone, and `xysFilter` to
+ * drop the slots too few spectra contain.
+ * @param data - data
+ * @param options - options.
+ * @returns the common `x` axis and one aligned `y` array per input spectrum.
+ */
+export function xyArrayAlignByIntensity(
+  data: DataXY[],
+  options: XYArrayAlignByIntensityOptions = {},
+): DataXYs<Float64Array> {
+  const { delta = 1 } = options;
+
+  const merged = mergeSortedXY(data);
+  const slots = getIntensitySlots(merged, { delta });
+
+  const x = slots.x;
+  const ys = Array.from(data, () => new Float64Array(x.length));
+  const { slotIndex } = slots;
+  const { spectrumIndex, y: mergedY } = merged;
+  for (let i = 0; i < slotIndex.length; i++) {
+    ys[spectrumIndex[i]][slotIndex[i]] += mergedY[i];
+  }
+
+  return { x, ys };
+}
diff --git a/src/xyArray/xyArrayFilterMinYValue.ts b/src/xyArray/xyArrayFilterMinYValue.ts
new file mode 100644
index 000000000..d15a7fe9e
--- /dev/null
+++ b/src/xyArray/xyArrayFilterMinYValue.ts
@@ -0,0 +1,50 @@
+import type { DataXY } from 'cheminfo-types';
+
+import { xyArrayMaxY } from './xyArrayMaxY.ts';
+
+export function xyArrayFilterMinYValue(
+  data: DataXY[],
+  minRelativeYValue: number,
+): Array<DataXY<number[]>>;
+export function xyArrayFilterMinYValue(
+  data: DataXY[],
+  minRelativeYValue?: number,
+): DataXY[];
+
+/**
+ * Removes from every spectrum the points whose y is below a fraction of the
+ * largest y of the whole set.
+ * Unlike `xyFilterMinYValue`, which uses the largest y of the spectrum it filters,
+ * the threshold here is the same for all the spectra, so a weak spectrum may lose
+ * every one of its points.
+ * @param data - data.
+ * @param minRelativeYValue - the minimum relative value compare to the global Y max value.
+ * @returns the spectra without their smallest points, the data itself when there is nothing to filter on.
+ */
+export function xyArrayFilterMinYValue(
+  data: DataXY[],
+  minRelativeYValue?: number,
+): DataXY[] {
+  if (minRelativeYValue === undefined) return data;
+
+  let total = 0;
+  for (const spectrum of data) {
+    total += spectrum.y.length;
+  }
+  const threshold = total === 0 ? 0 : xyArrayMaxY(data) * minRelativeYValue;
+
+  // The threshold is a number and not a callback, so the comparison is inlined
+  // rather than going through `xyFilter`, which costs 25% on 30 million points.
+  return data.map((spectrum) => {
+    const { x, y } = spectrum;
+    const newX: number[] = [];
+    const newY: number[] = [];
+    for (let i = 0; i < y.length; i++) {
+      if (y[i] >= threshold) {
+        newX.push(x[i]);
+        newY.push(y[i]);
+      }
+    }
+    return { x: newX, y: newY };
+  });
+}
diff --git a/src/xyArray/xyArrayMaxY.ts b/src/xyArray/xyArrayMaxY.ts
new file mode 100644
index 000000000..0ef1886c0
--- /dev/null
+++ b/src/xyArray/xyArrayMaxY.ts
@@ -0,0 +1,26 @@
+import type { DataXY } from 'cheminfo-types';
+
+import { xMaxValue } from '../x/index.ts';
+
+/**
+ * Returns the largest y of all the spectra.
+ * @param data - data.
+ * @returns the largest y.
+ */
+export function xyArrayMaxY(data: DataXY[]): number {
+  let max = Number.NEGATIVE_INFINITY;
+  let isEmpty = true;
+
+  for (const spectrum of data) {
+    if (spectrum.y.length === 0) continue;
+    isEmpty = false;
+    const value = xMaxValue(spectrum.y);
+    if (value > max) max = value;
+  }
+
+  if (isEmpty) {
+    throw new Error('can not process empty arrays');
+  }
+
+  return max;
+}
diff --git a/src/xyArray/xyArrayMergeByIntensity.ts b/src/xyArray/xyArrayMergeByIntensity.ts
new file mode 100644
index 000000000..a66d4e0e7
--- /dev/null
+++ b/src/xyArray/xyArrayMergeByIntensity.ts
@@ -0,0 +1,49 @@
+import type { DataXY } from 'cheminfo-types';
+
+import type { IntensitySlots } from './utils/getIntensitySlots.ts';
+import { getIntensitySlots } from './utils/getIntensitySlots.ts';
+import { mergeSortedXY } from './utils/mergeSortedXY.ts';
+
+export interface XYArrayMergeByIntensityOptions {
+  /**
+   * A peak joins a slot when its x is within `delta` of the x of the most intense
+   * peak of that slot, so a slot is never wider than `2 * delta`. It may also be a
+   * function that allows to change `delta` depending on the x value.
+   * @default 1
+   */
+  delta?: ((arg: number) => number) | number;
+}
+
+/**
+ * The merged peaks: every field of a slot except which slot each input point landed in.
+ */
+export type XYArrayMergeByIntensityResult = Omit<IntensitySlots, 'slotIndex'>;
+
+/**
+ * Merges the peaks of many spectra around the most intense ones.
+ *
+ * This is `xyArrayAlignByIntensity` without the per spectrum y arrays: the
+ * most intense peak of the whole set takes every peak within `delta` of it, the
+ * most intense of the remaining peaks opens the next slot, and so on. Every peak
+ * ends up in exactly one slot, so the sum of `y` is the sum of the intensities of
+ * all the spectra.
+ *
+ * The x values of each spectrum must be ascending. Only the merged peaks are held
+ * in memory, so it scales to tens of millions of peaks where a dense alignment
+ * matrix does not.
+ * @param data - data
+ * @param options - options.
+ * @returns the common peaks, ascending by x.
+ */
+export function xyArrayMergeByIntensity(
+  data: DataXY[],
+  options: XYArrayMergeByIntensityOptions = {},
+): XYArrayMergeByIntensityResult {
+  const { delta = 1 } = options;
+
+  const { x, y, from, to, number } = getIntensitySlots(mergeSortedXY(data), {
+    delta,
+  });
+
+  return { x, y, from, to, number };
+}
diff --git a/src/xyObject/__tests__/xyObjectSortX.test.ts b/src/xyObject/__tests__/xyObjectSortX.test.ts
index f2780f27e..e3c0882f6 100644
--- a/src/xyObject/__tests__/xyObjectSortX.test.ts
+++ b/src/xyObject/__tests__/xyObjectSortX.test.ts
@@ -1,3 +1,4 @@
+import { XSadd } from 'ml-xsadd';
 import { expect, test } from 'vitest';
 
 import { xyObjectSortX } from '../xyObjectSortX.ts';
@@ -17,3 +18,21 @@ test('xyObjectSortX', () => {
     { x: 3, y: 3 },
   ]);
 });
+
+// The short array above is ordered by a comparator, the long one by xGetSortOrder: both have
+// to give the same stable order, and both in place.
+test('same stable order whatever the length', () => {
+  for (const length of [100, 1000]) {
+    const { random } = new XSadd(42);
+    const points = Array.from({ length }, (value, index) => ({
+      x: Math.round(random() * 20),
+      y: index,
+    }));
+    const expected = points.toSorted((a, b) => a.x - b.x);
+
+    const result = xyObjectSortX(points);
+
+    expect(result).toBe(points);
+    expect(result).toStrictEqual(expected);
+  }
+});
diff --git a/src/xyObject/xyObjectSortX.ts b/src/xyObject/xyObjectSortX.ts
index 95787951c..07d235934 100644
--- a/src/xyObject/xyObjectSortX.ts
+++ b/src/xyObject/xyObjectSortX.ts
@@ -1,4 +1,6 @@
 import type { Point } from '../types/index.ts';
+import { xGetSortOrder } from '../x/index.ts';
+import { RADIX_ORDER_MIN_LENGTH } from '../x/utils/radixSortOrder.ts';
 
 /**
  * Sorts an array of points in-place.
@@ -6,6 +8,26 @@ import type { Point } from '../types/index.ts';
  * @returns sorted array of points {x,y}.
  */
 export function xyObjectSortX(points: Point[]): Point[] {
-  points.sort((a, b) => a.x - b.x);
+  const { length } = points;
+
+  if (length < RADIX_ORDER_MIN_LENGTH) {
+    points.sort((a, b) => a.x - b.x);
+    return points;
+  }
+
+  const x = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    x[i] = points[i].x;
+  }
+
+  const order = xGetSortOrder(x);
+  const sorted = new Array<Point>(length);
+  for (let i = 0; i < length; i++) {
+    sorted[i] = points[order[i]];
+  }
+  for (let i = 0; i < length; i++) {
+    points[i] = sorted[i];
+  }
+
   return points;
 }
diff --git a/src/xys/__tests__/xysFilter.test.ts b/src/xys/__tests__/xysFilter.test.ts
new file mode 100644
index 000000000..8acd733d4
--- /dev/null
+++ b/src/xys/__tests__/xysFilter.test.ts
@@ -0,0 +1,173 @@
+import { expect, test } from 'vitest';
+
+import { xyArrayAlign } from '../../xyArray/xyArrayAlign.ts';
+import { xysFilter } from '../xysFilter.ts';
+
+/*
+ * Eight spectra, four slots.
+ *          x=10  x=20  x=30  x=40
+ * present    4     3     3     2   spectra
+ * longest    1     3     3     2   run of consecutive spectra
+ */
+function aligned() {
+  return {
+    x: [10, 20, 30, 40],
+    ys: [
+      [1, 0, 0, 1],
+      [0, 0, 1, 1],
+      [1, 0, 1, 0],
+      [0, 0, 1, 0],
+      [1, 0, 0, 0],
+      [0, 1, 0, 0],
+      [1, 1, 0, 0],
+      [0, 1, 0, 0],
+    ],
+  };
+}
+
+test('without options nothing is removed', () => {
+  const data = aligned();
+  const result = xysFilter(data);
+
+  expect(Array.from(result.x)).toStrictEqual(data.x);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual(data.ys);
+});
+
+test('minY on its own removes nothing', () => {
+  const data = aligned();
+  const result = xysFilter(data, { minY: 5 });
+
+  expect(Array.from(result.x)).toStrictEqual(data.x);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual(data.ys);
+});
+
+test('minNumberOfSpectra counts the spectra wherever they are', () => {
+  const result = xysFilter(aligned(), { minNumberOfSpectra: 3 });
+
+  expect(Array.from(result.x)).toStrictEqual([10, 20, 30]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [1, 0, 0],
+    [0, 0, 1],
+    [1, 0, 1],
+    [0, 0, 1],
+    [1, 0, 0],
+    [0, 1, 0],
+    [1, 1, 0],
+    [0, 1, 0],
+  ]);
+});
+
+test('minConsecutiveSpectra rejects the scattered slot the count kept', () => {
+  // x=10 is in four spectra but never in two neighbouring ones.
+  const result = xysFilter(aligned(), { minConsecutiveSpectra: 3 });
+
+  expect(Array.from(result.x)).toStrictEqual([20, 30]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([
+    [0, 0],
+    [0, 1],
+    [0, 1],
+    [0, 1],
+    [0, 0],
+    [1, 0],
+    [1, 0],
+    [1, 0],
+  ]);
+});
+
+test('two consecutive spectra are enough for the narrow slot', () => {
+  const result = xysFilter(aligned(), { minConsecutiveSpectra: 2 });
+
+  expect(Array.from(result.x)).toStrictEqual([20, 30, 40]);
+});
+
+test('the criteria are combined', () => {
+  const result = xysFilter(aligned(), {
+    minNumberOfSpectra: 3,
+    minConsecutiveSpectra: 2,
+  });
+
+  expect(Array.from(result.x)).toStrictEqual([20, 30]);
+});
+
+test('requiredY keeps only the slots every spectrum contains', () => {
+  const result = xysFilter(
+    {
+      x: [1, 2],
+      ys: [
+        [1, 1],
+        [0, 1],
+        [3, 1],
+      ],
+    },
+    { requiredY: true },
+  );
+
+  expect(Array.from(result.x)).toStrictEqual([2]);
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual([[1], [1], [1]]);
+});
+
+test('a negative peak counts as present until minY is set', () => {
+  const data = {
+    x: [1, 2],
+    ys: [
+      [-5, 0],
+      [3, 4],
+    ],
+  };
+
+  expect(Array.from(xysFilter(data, { requiredY: true }).x)).toStrictEqual([1]);
+  // Setting minY compares the y itself, so the negative peak stops counting.
+  expect(
+    Array.from(xysFilter(data, { requiredY: true, minY: 0 }).x),
+  ).toStrictEqual([]);
+});
+
+test('minY sets how intense a peak must be to count', () => {
+  const data = {
+    x: [1, 2],
+    ys: [
+      [5, 1],
+      [5, 1],
+      [5, 1],
+    ],
+  };
+
+  expect(
+    Array.from(xysFilter(data, { minNumberOfSpectra: 3 }).x),
+  ).toStrictEqual([1, 2]);
+  expect(
+    Array.from(xysFilter(data, { minNumberOfSpectra: 3, minY: 2 }).x),
+  ).toStrictEqual([1]);
+});
+
+test('agrees with the requiredY option of xyArrayAlign', () => {
+  const data = [
+    { x: [1, 2, 3], y: [1, 1, 1] },
+    { x: [0.1, 1.1, 2.1, 3.1, 4.1], y: [1, 1, 1, 1, 1] },
+    { x: [2.9, 3.1, 3.9, 4.9], y: [1, 1, 1, 1] },
+  ];
+  const expected = xyArrayAlign(data, { delta: 0.15, requiredY: true });
+  const result = xysFilter(xyArrayAlign(data, { delta: 0.15 }), {
+    requiredY: true,
+  });
+
+  expect(Array.from(result.x)).toStrictEqual(Array.from(expected.x));
+  expect(result.ys.map((y) => Array.from(y))).toStrictEqual(
+    expected.ys.map((y) => Array.from(y)),
+  );
+});
+
+test('no slot survives', () => {
+  const result = xysFilter(aligned(), { minNumberOfSpectra: 8 });
+
+  expect(Array.from(result.x)).toStrictEqual([]);
+  expect(result.ys).toHaveLength(8);
+  expect(Array.from(result.ys[0])).toStrictEqual([]);
+});
+
+test('no spectra', () => {
+  const result = xysFilter({ x: [1, 2], ys: [] }, { requiredY: true });
+
+  expect(Array.from(result.x)).toStrictEqual([1, 2]);
+  expect(result.ys).toStrictEqual([]);
+});
diff --git a/src/xys/index.ts b/src/xys/index.ts
new file mode 100644
index 000000000..b52c468a7
--- /dev/null
+++ b/src/xys/index.ts
@@ -0,0 +1 @@
+export * from './xysFilter.ts';
diff --git a/src/xys/xysFilter.ts b/src/xys/xysFilter.ts
new file mode 100644
index 000000000..0d614ef70
--- /dev/null
+++ b/src/xys/xysFilter.ts
@@ -0,0 +1,120 @@
+import type { DataXYs } from 'cheminfo-types';
+
+export interface XYsFilterOptions {
+  /**
+   * The minimal y a spectrum must have in a slot to contain it, meant to be set
+   * above 0. When it is not set a spectrum contains a slot when its y there is
+   * not 0, so the negative peaks of a signed spectrum count as present.
+   * On its own this never removes a slot, it only changes what the criteria below
+   * count as contained.
+   * @default undefined
+   */
+  minY?: number;
+
+  /**
+   * The smallest number of spectra that must contain a slot. 0 does not filter.
+   * @default 0
+   */
+  minNumberOfSpectra?: number;
+
+  /**
+   * The smallest number of neighbouring spectra that must all contain a slot,
+   * neighbouring meaning next to each other in `ys`. Made for spectra ordered by
+   * retention time, where a real compound is present in a run of consecutive
+   * scans while noise is scattered. 0 does not filter.
+   * @default 0
+   */
+  minConsecutiveSpectra?: number;
+
+  /**
+   * If true, every spectrum must contain the slot. Same as setting
+   * `minNumberOfSpectra` to the number of spectra.
+   * @default false
+   */
+  requiredY?: boolean;
+}
+
+/**
+ * Removes the slots of an aligned set of spectra that too few spectra contain.
+ *
+ * `requiredY`, `minNumberOfSpectra` and `minConsecutiveSpectra` are the criteria,
+ * and `minY` only says what they count as contained. They may be combined, a slot
+ * is kept only when it passes all the ones that are set, so with no criterion set
+ * nothing is removed. They only ever remove slots, the y values are never changed.
+ * @param data - aligned spectra, as returned by the align methods.
+ * @param options - options.
+ * @returns the same spectra with the rejected slots removed.
+ */
+export function xysFilter(
+  data: DataXYs,
+  options: XYsFilterOptions = {},
+): DataXYs<Float64Array> {
+  const {
+    minY,
+    minNumberOfSpectra = 0,
+    minConsecutiveSpectra = 0,
+    requiredY = false,
+  } = options;
+
+  const { x, ys } = data;
+  const length = x.length;
+  const numberOfSpectra = ys.length;
+  const minNumber = requiredY ? numberOfSpectra : minNumberOfSpectra;
+  const hasMinY = minY !== undefined;
+  const threshold = minY ?? 0;
+
+  const numbers = new Uint32Array(length);
+  let longestRuns: Uint32Array | undefined;
+  if (minConsecutiveSpectra > 0) {
+    longestRuns = new Uint32Array(length);
+    const runs = new Uint32Array(length);
+    for (let j = 0; j < numberOfSpectra; j++) {
+      const y = ys[j];
+      for (let i = 0; i < length; i++) {
+        if (hasMinY ? y[i] > threshold : y[i] !== 0) {
+          numbers[i]++;
+          const run = runs[i] + 1;
+          runs[i] = run;
+          if (run > longestRuns[i]) longestRuns[i] = run;
+        } else {
+          runs[i] = 0;
+        }
+      }
+    }
+  } else if (minNumber > 0) {
+    for (let j = 0; j < numberOfSpectra; j++) {
+      const y = ys[j];
+      for (let i = 0; i < length; i++) {
+        if (hasMinY ? y[i] > threshold : y[i] !== 0) numbers[i]++;
+      }
+    }
+  }
+
+  let kept = 0;
+  for (let i = 0; i < length; i++) {
+    if (numbers[i] < minNumber) continue;
+    if (longestRuns && longestRuns[i] < minConsecutiveSpectra) continue;
+    kept++;
+  }
+
+  const newX = new Float64Array(kept);
+  const positions = new Uint32Array(kept);
+  let at = 0;
+  for (let i = 0; i < length; i++) {
+    if (numbers[i] < minNumber) continue;
+    if (longestRuns && longestRuns[i] < minConsecutiveSpectra) continue;
+    newX[at] = x[i];
+    positions[at] = i;
+    at++;
+  }
+
+  const newYs = Array.from(ys, (y) => {
+    const newY = new Float64Array(kept);
+    for (let i = 0; i < kept; i++) {
+      newY[i] = y[positions[i]];
+    }
+    return newY;
+  });
+
+  return { x: newX, ys: newYs };
+}