diff --git a/benchmark/all.js b/benchmark/all.js
new file mode 100644
index 0000000..2402f43
--- /dev/null
+++ b/benchmark/all.js
@@ -0,0 +1,84 @@
+/*
+ * Benchmarks every distance and similarity function so their relative cost can
+ * be compared. Use beforeAfter.js instead to compare two implementations of the
+ * same function.
+ *
+ * Run with: node benchmark/all.js  (or: bun benchmark/all.js)
+ */
+import Benchmark from 'benchmark';
+
+import { distance, similarity } from '../src/index.ts';
+
+const LENGTH = 10000;
+
+function makeVector(seed) {
+  const vector = new Float64Array(LENGTH);
+  let state = seed;
+  for (let i = 0; i < LENGTH; i++) {
+    state = (state * 1103515245 + 12345) % 2147483648;
+    vector[i] = 0.1 + (state / 2147483648) * 2;
+  }
+  return vector;
+}
+
+const a = makeVector(42);
+const b = makeVector(1337);
+
+const entries = [];
+for (const [namespace, functions] of [
+  ['distance', distance],
+  ['similarity', similarity],
+]) {
+  for (const [name, callback] of Object.entries(functions)) {
+    if (typeof callback !== 'function') continue;
+    if (name === 'minkowski') {
+      for (const p of [1, 2, 3]) {
+        entries.push([
+          `${namespace}.${name}(p=${p})`,
+          (x, y) => callback(x, y, p),
+        ]);
+      }
+      continue;
+    }
+    entries.push([`${namespace}.${name}`, callback]);
+  }
+}
+
+function log(message) {
+  // eslint-disable-next-line no-console -- benchmark output
+  console.log(message);
+}
+
+const results = [];
+const suite = new Benchmark.Suite();
+for (const [name, callback] of entries) {
+  suite.add(name, () => callback(a, b), { minSamples: 30 });
+}
+
+suite
+  .on('cycle', (event) => {
+    const { name, hz, stats } = event.target;
+    results.push({
+      name,
+      nanoseconds: 1e9 / hz,
+      rme: stats.rme,
+      samples: stats.sample.length,
+    });
+  })
+  .on('complete', () => {
+    results.sort((first, second) => second.nanoseconds - first.nanoseconds);
+    log(`\nn = ${LENGTH}, sorted slowest first\n`);
+    log(
+      `${'function'.padEnd(32)}${'per call'.padStart(12)}${'per element'.padStart(14)}${'error'.padStart(9)}   result`,
+    );
+    const values = new Map(entries.map(([name, callback]) => [name, callback]));
+    for (const { name, nanoseconds, rme } of results) {
+      // A wide confidence interval means the engine kept re-tiering this one;
+      // treat the number as indicative only.
+      const flag = rme > 10 ? ' (!)' : '';
+      log(
+        `${name.padEnd(32)}${`${(nanoseconds / 1000).toFixed(1)} µs`.padStart(12)}${`${(nanoseconds / LENGTH).toFixed(2)} ns`.padStart(14)}${`±${rme.toFixed(1)}%${flag}`.padStart(13)}   ${values.get(name)(a, b)}`,
+      );
+    }
+  })
+  .run({ async: false });
diff --git a/benchmark/arrayKinds.js b/benchmark/arrayKinds.js
new file mode 100644
index 0000000..5936682
--- /dev/null
+++ b/benchmark/arrayKinds.js
@@ -0,0 +1,105 @@
+/*
+ * Cost of `(ai * bi) / (ai + bi)` depending on the element representation.
+ * One kind per process, so the loads stay monomorphic.
+ *
+ * `cached` copies `x[i]` into a local, `direct` reads it on every use.
+ *
+ * Run with: node benchmark/arrayKinds.js <array|typed|mixed> <cached|direct>
+ */
+import { argv } from 'node:process';
+
+const LENGTH = 10000;
+const WARMUP_MS = 1000;
+const TARGET_MS = 2000; // keep at 1000 or more, shorter runs are too noisy
+
+const kind = argv[2] ?? 'typed';
+const reads = argv[3] ?? 'cached';
+
+function fill(target) {
+  let state = 42;
+  for (let i = 0; i < target.length; i++) {
+    state = (state * 1103515245 + 12345) % 2147483648;
+    target[i] = 0.1 + (state / 2147483648) * 2;
+  }
+  return target;
+}
+
+function plainVector() {
+  return fill(Array.from({ length: LENGTH }, () => 0));
+}
+
+function typedVector() {
+  return fill(new Float64Array(LENGTH));
+}
+
+let a;
+let b;
+if (kind === 'array') {
+  a = plainVector();
+  b = plainVector();
+} else if (kind === 'typed') {
+  a = typedVector();
+  b = typedVector();
+} else if (kind === 'mixed') {
+  a = typedVector();
+  b = plainVector();
+} else {
+  throw new Error(`unknown kind: ${kind}`);
+}
+
+function cachedKernel(x, y) {
+  let sum = 0;
+  for (let i = 0; i < x.length; i++) {
+    const xi = x[i];
+    const yi = y[i];
+    sum += (xi * yi) / (xi + yi);
+  }
+  return sum;
+}
+
+function directKernel(x, y) {
+  let sum = 0;
+  for (let i = 0; i < x.length; i++) {
+    sum += (x[i] * y[i]) / (x[i] + y[i]);
+  }
+  return sum;
+}
+
+const kernel = reads === 'cached' ? cachedKernel : directKernel;
+
+let sink = 0;
+// the arguments are swapped after each call, so `mixed` really sees both orders
+let x = a;
+let y = b;
+
+let start = performance.now();
+while (performance.now() - start < WARMUP_MS) {
+  sink += kernel(x, y);
+  [x, y] = [y, x];
+}
+
+let rounds = 0;
+let elapsed = 0;
+start = performance.now();
+while (elapsed < TARGET_MS) {
+  sink += kernel(x, y);
+  [x, y] = [y, x];
+  rounds++;
+  elapsed = performance.now() - start;
+}
+
+// keeps the loop from being dropped as dead code
+if (!Number.isFinite(sink)) throw new Error('kernel diverged');
+
+const operationsPerSecond = (rounds * LENGTH * 1000) / elapsed;
+// eslint-disable-next-line no-console -- benchmark output
+console.log(
+  [
+    kind,
+    reads,
+    operationsPerSecond.toFixed(0),
+    LENGTH,
+    kernel(a, b).toFixed(10),
+    kernel(b, a).toFixed(10),
+  ].join('\t'),
+);
diff --git a/benchmark/arrayKinds.sh b/benchmark/arrayKinds.sh
new file mode 100755
index 0000000..7f15c1f
--- /dev/null
+++ b/benchmark/arrayKinds.sh
@@ -0,0 +1,41 @@
+#!/bin/bash
+# Runs every element representation against both read styles, one process each,
+# then prints a summary table.
+set -euo pipefail
+
+directory="$(dirname "$0")"
+measurements=""
+
+for kind in array typed mixed; do
+  for reads in cached direct; do
+    measurements="${measurements}$(node "$directory/arrayKinds.js" "$kind" "$reads")
+"
+  done
+done
+
+printf '%s' "$measurements" | awk -F'\t' '
+  {
+    speed[$1 "/" $2] = $3 / 1e6
+    length_ = $4
+    results[$5] = 1
+    results[$6] = 1
+  }
+  END {
+    printf "\n  (ai * bi) / (ai + bi) on %d elements, arguments swapped at each call\n\n", length_
+    printf "  %-8s %14s %14s %10s\n", "kind", "cached", "direct", "ratio"
+    printf "  %-8s %14s %14s %10s\n", "--------", "--------------", "--------------", "----------"
+    split("array typed mixed", kinds, " ")
+    for (i = 1; i <= 3; i++) {
+      cached = speed[kinds[i] "/cached"]
+      direct = speed[kinds[i] "/direct"]
+      printf "  %-8s %11.0f M/s %11.0f M/s %9.2fx\n", kinds[i], cached, direct, cached / direct
+    }
+    distinct = 0
+    for (result in results) { distinct++; sample = result }
+    if (distinct == 1) {
+      printf "\n  all results identical: %s\n\n", sample
+    } else {
+      printf "\n  WARNING: %d different results\n\n", distinct
+    }
+  }
+'
diff --git a/benchmark/beforeAfter.js b/benchmark/beforeAfter.js
new file mode 100644
index 0000000..f80dbab
--- /dev/null
+++ b/benchmark/beforeAfter.js
@@ -0,0 +1,373 @@
+/*
+ * Compares the previous implementation of each optimized function with the
+ * current one. Both versions are copied in here so that they run in the same
+ * process on the same data.
+ *
+ * Run with: node benchmark/beforeAfter.js
+ */
+import Benchmark from 'benchmark';
+
+const LENGTHS = [1000, 10000, 100000];
+const KINDS = ['Array', 'Float64Array'];
+
+function makeArrays(kind, length) {
+  const a =
+    kind === 'Float64Array' ? new Float64Array(length) : new Array(length);
+  const b =
+    kind === 'Float64Array' ? new Float64Array(length) : new Array(length);
+  let state = 42;
+  for (let i = 0; i < length; i++) {
+    state = (state * 1103515245 + 12345) % 2147483648;
+    a[i] = 0.1 + (state / 2147483648) * 2;
+    state = (state * 1103515245 + 12345) % 2147483648;
+    b[i] = 0.1 + (state / 2147483648) * 2;
+  }
+  return [a, b];
+}
+
+/* ------------------------------------------------ distances/pearson */
+function pearsonBefore(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += ((a[i] - b[i]) * (a[i] - b[i])) / b[i];
+  }
+  return d;
+}
+function pearsonAfter(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    const bi = b[i];
+    const diff = a[i] - bi;
+    d += (diff * diff) / bi;
+  }
+  return d;
+}
+
+/* ------------------------------------------------ distances/squared */
+function squaredBefore(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]);
+  }
+  return d;
+}
+function squaredAfter(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const diff = ai - bi;
+    d += (diff * diff) / (ai + bi);
+  }
+  return d;
+}
+
+/* ------------------------------------------------ similarities/cosine */
+function cosineBefore(a, b) {
+  let p = 0;
+  let p2 = 0;
+  let q2 = 0;
+  for (let i = 0; i < a.length; i++) {
+    p += a[i] * b[i];
+    p2 += a[i] * a[i];
+    q2 += b[i] * b[i];
+  }
+  return p / (Math.sqrt(p2) * Math.sqrt(q2));
+}
+function cosineAfter(a, b) {
+  let p = 0;
+  let p2 = 0;
+  let q2 = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    p += ai * bi;
+    p2 += ai * ai;
+    q2 += bi * bi;
+  }
+  return p / (Math.sqrt(p2) * Math.sqrt(q2));
+}
+
+/* ------------------------------------------------ distances/clark */
+function clarkBefore(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += (Math.abs(a[i] - b[i]) / (a[i] + b[i])) ** 2;
+  }
+  return Math.sqrt(d);
+}
+function clarkAfter(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const ratio = (ai - bi) / (ai + bi);
+    d += ratio * ratio;
+  }
+  return Math.sqrt(d);
+}
+
+/* ------------------------------------------------ distances/kumarJohnson */
+function kumarJohnsonBefore(a, b) {
+  let ans = 0;
+  for (let i = 0; i < a.length; i++) {
+    ans += (a[i] * a[i] - b[i] * b[i]) ** 2 / (2 * (a[i] * b[i]) ** 1.5);
+  }
+  return ans;
+}
+function kumarJohnsonAfter(a, b) {
+  let ans = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const numerator = ai * ai - bi * bi;
+    const prod = ai * bi;
+    ans += (numerator * numerator) / (2 * prod * Math.sqrt(prod));
+  }
+  return ans;
+}
+
+/* ------------------------------------------------ distances/minkowski
+ * Each order gets its own copy on purpose. Calling one shared `minkowski` with
+ * p = 1 and p = 2 makes the exponent polymorphic and lets the two orders share
+ * inline caches, which hides the difference the branches are meant to measure.
+ */
+function minkowskiBeforeP1(a, b, p) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += Math.abs(a[i] - b[i]) ** p;
+  }
+  return d ** (1 / p);
+}
+function minkowskiAfterP1(a, b, p) {
+  let d = 0;
+  if (p === 1) {
+    for (let i = 0; i < a.length; i++) {
+      d += Math.abs(a[i] - b[i]);
+    }
+    return d;
+  }
+  if (p === 2) {
+    for (let i = 0; i < a.length; i++) {
+      const diff = a[i] - b[i];
+      d += diff * diff;
+    }
+    return Math.sqrt(d);
+  }
+  for (let i = 0; i < a.length; i++) {
+    d += Math.abs(a[i] - b[i]) ** p;
+  }
+  return d ** (1 / p);
+}
+function minkowskiBeforeP2(a, b, p) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += Math.abs(a[i] - b[i]) ** p;
+  }
+  return d ** (1 / p);
+}
+function minkowskiAfterP2(a, b, p) {
+  let d = 0;
+  if (p === 1) {
+    for (let i = 0; i < a.length; i++) {
+      d += Math.abs(a[i] - b[i]);
+    }
+    return d;
+  }
+  if (p === 2) {
+    for (let i = 0; i < a.length; i++) {
+      const diff = a[i] - b[i];
+      d += diff * diff;
+    }
+    return Math.sqrt(d);
+  }
+  for (let i = 0; i < a.length; i++) {
+    d += Math.abs(a[i] - b[i]) ** p;
+  }
+  return d ** (1 / p);
+}
+
+/* ------------------------------------------------ similarities/pearson */
+function meanOf(input) {
+  let sumValue = 0;
+  for (const value of input) sumValue += value;
+  return sumValue / input.length;
+}
+// A private copy of cosine: sharing cosineBefore with the cosine pair would
+// feed it both Float64Array and Array inputs and make its loads polymorphic.
+function cosineForPearson(a, b) {
+  let p = 0;
+  let p2 = 0;
+  let q2 = 0;
+  for (let i = 0; i < a.length; i++) {
+    p += a[i] * b[i];
+    p2 += a[i] * a[i];
+    q2 += b[i] * b[i];
+  }
+  return p / (Math.sqrt(p2) * Math.sqrt(q2));
+}
+function pearsonSimilarityBefore(a, b) {
+  const avgA = meanOf(a);
+  const avgB = meanOf(b);
+  const newA = new Array(a.length);
+  const newB = new Array(b.length);
+  for (let i = 0; i < newA.length; i++) {
+    newA[i] = a[i] - avgA;
+    newB[i] = b[i] - avgB;
+  }
+  return cosineForPearson(newA, newB);
+}
+function pearsonSimilarityAfter(a, b) {
+  const length = a.length;
+  let sumA = 0;
+  let sumB = 0;
+  for (let i = 0; i < length; i++) {
+    sumA += a[i];
+    sumB += b[i];
+  }
+  const avgA = sumA / length;
+  const avgB = sumB / length;
+  let p = 0;
+  let p2 = 0;
+  let q2 = 0;
+  for (let i = 0; i < length; i++) {
+    const centredA = a[i] - avgA;
+    const centredB = b[i] - avgB;
+    p += centredA * centredB;
+    p2 += centredA * centredA;
+    q2 += centredB * centredB;
+  }
+  return p / (Math.sqrt(p2) * Math.sqrt(q2));
+}
+
+/* ------------------------------------------------ topsoe */
+function topsoeBefore(a, b) {
+  let ans = 0;
+  for (let i = 0; i < a.length; i++) {
+    ans +=
+      a[i] * Math.log((2 * a[i]) / (a[i] + b[i])) +
+      b[i] * Math.log((2 * b[i]) / (a[i] + b[i]));
+  }
+  return ans;
+}
+function topsoeAfter(a, b) {
+  let ans = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const sum = ai + bi;
+    ans += ai * Math.log((2 * ai) / sum) + bi * Math.log((2 * bi) / sum);
+  }
+  return ans;
+}
+
+const PAIRS = [
+  { name: 'distances/pearson', before: pearsonBefore, after: pearsonAfter },
+  { name: 'distances/squared', before: squaredBefore, after: squaredAfter },
+  { name: 'similarities/cosine', before: cosineBefore, after: cosineAfter },
+  { name: 'distances/clark', before: clarkBefore, after: clarkAfter },
+  { name: 'distances/topsoe', before: topsoeBefore, after: topsoeAfter },
+  {
+    name: 'distances/kumarJohnson',
+    before: kumarJohnsonBefore,
+    after: kumarJohnsonAfter,
+  },
+  {
+    name: 'distances/minkowski p=1',
+    before: (a, b) => minkowskiBeforeP1(a, b, 1),
+    after: (a, b) => minkowskiAfterP1(a, b, 1),
+  },
+  {
+    name: 'distances/minkowski p=2',
+    before: (a, b) => minkowskiBeforeP2(a, b, 2),
+    after: (a, b) => minkowskiAfterP2(a, b, 2),
+  },
+  {
+    name: 'similarities/pearson',
+    before: pearsonSimilarityBefore,
+    after: pearsonSimilarityAfter,
+  },
+];
+
+function log(message) {
+  // eslint-disable-next-line no-console -- benchmark output
+  console.log(message);
+}
+
+function runPair({ name, before, after }, a, b) {
+  return new Promise((resolve) => {
+    const beforeValue = before(a, b);
+    const afterValue = after(a, b);
+    const identical = Object.is(beforeValue, afterValue);
+    let beforeHz = 0;
+    let afterHz = 0;
+    let beforeRme = 0;
+    let afterRme = 0;
+    new Benchmark.Suite(name)
+      .add('before', () => before(a, b), { maxTime: 2 })
+      .add('after', () => after(a, b), { maxTime: 2 })
+      .on('cycle', (event) => {
+        const { name: which, hz, stats } = event.target;
+        if (which === 'before') {
+          beforeHz = hz;
+          beforeRme = stats.rme;
+        } else {
+          afterHz = hz;
+          afterRme = stats.rme;
+        }
+        log(
+          `  ${which.padEnd(7)}${(1e3 / hz).toFixed(4).padStart(10)} ms/op  ±${stats.rme.toFixed(2)}%  (${stats.sample.length} samples)`,
+        );
+      })
+      .on('complete', () => {
+        const speedup = afterHz / beforeHz;
+        log(
+          `  => ${speedup.toFixed(2)}x   ${identical ? 'identical result' : `DIFFERS ${beforeValue} -> ${afterValue}`}\n`,
+        );
+        resolve({ speedup, identical, rme: Math.max(beforeRme, afterRme) });
+      })
+      .run({ async: false });
+  });
+}
+
+const summary = [];
+for (const length of LENGTHS) {
+  for (const kind of KINDS) {
+    const [a, b] = makeArrays(kind, length);
+    log(`\n=== ${kind}, ${length} elements ===\n`);
+    for (const pair of PAIRS) {
+      log(pair.name);
+      // eslint-disable-next-line no-await-in-loop -- suites must not overlap
+      const result = await runPair(pair, a, b);
+      summary.push({ ...result, name: pair.name, length, kind });
+    }
+  }
+}
+
+log('\n\n=== speedup summary (after / before) ===\n');
+for (const kind of KINDS) {
+  log(kind);
+  log(
+    `  ${'function'.padEnd(24)}${LENGTHS.map((l) => `${l}`.padStart(10)).join('')}`,
+  );
+  for (const pair of PAIRS) {
+    const cells = LENGTHS.map((length) => {
+      const found = summary.find(
+        (entry) =>
+          entry.name === pair.name &&
+          entry.length === length &&
+          entry.kind === kind,
+      );
+      return `${found.speedup.toFixed(2)}x`.padStart(10);
+    });
+    log(`  ${pair.name.padEnd(24)}${cells.join('')}`);
+  }
+  log('');
+}
+
+const differing = summary.filter((entry) => !entry.identical);
+log(
+  differing.length === 0
+    ? 'all results identical'
+    : `results differ for: ${[...new Set(differing.map((entry) => entry.name))].join(', ')}`,
+);
diff --git a/benchmark/elementReads.js b/benchmark/elementReads.js
new file mode 100644
index 0000000..f2b5907
--- /dev/null
+++ b/benchmark/elementReads.js
@@ -0,0 +1,153 @@
+/*
+ * Why does caching `a[i]` in a local help at all?
+ *
+ * Not because of the `NumberArray` TypeScript type: types are erased and V8
+ * never sees them. What matters is the element representation the load site
+ * actually observes at run time. On a site that has only ever seen one
+ * Float64Array shape, V8 eliminates the repeated loads by itself and the
+ * caching is worth nothing. On plain arrays, or on a site fed more than one
+ * array kind, the loads survive and the caching pays.
+ *
+ * Run with: node benchmark/elementReads.js  (or: bun benchmark/elementReads.js)
+ */
+import Benchmark from 'benchmark';
+
+const LENGTH = 10000;
+
+function fill(target) {
+  let state = 42;
+  for (let i = 0; i < target.length; i++) {
+    state = (state * 1103515245 + 12345) % 2147483648;
+    target[i] = 0.1 + (state / 2147483648) * 2;
+  }
+  return target;
+}
+
+const typedA = fill(new Float64Array(LENGTH));
+const typedB = fill(new Float64Array(LENGTH));
+const plainA = fill(Array.from({ length: LENGTH }, () => 0));
+const plainB = fill(Array.from({ length: LENGTH }, () => 0));
+
+/*
+ * Six copies of the same two loops. They must not be shared: a single copy
+ * called with several array kinds would make every measurement polymorphic,
+ * which is exactly the variable under test.
+ */
+function typedBefore(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]);
+  }
+  return d;
+}
+function typedAfter(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const diff = ai - bi;
+    d += (diff * diff) / (ai + bi);
+  }
+  return d;
+}
+function plainBefore(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]);
+  }
+  return d;
+}
+function plainAfter(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const diff = ai - bi;
+    d += (diff * diff) / (ai + bi);
+  }
+  return d;
+}
+function mixedBefore(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]);
+  }
+  return d;
+}
+function mixedAfter(a, b) {
+  let d = 0;
+  for (let i = 0; i < a.length; i++) {
+    const ai = a[i];
+    const bi = b[i];
+    const diff = ai - bi;
+    d += (diff * diff) / (ai + bi);
+  }
+  return d;
+}
+
+for (let i = 0; i < 5000; i++) {
+  typedBefore(typedA, typedB);
+  typedAfter(typedA, typedB);
+  plainBefore(plainA, plainB);
+  plainAfter(plainA, plainB);
+  mixedBefore(typedA, typedB);
+  mixedBefore(plainA, plainB);
+  mixedAfter(typedA, typedB);
+  mixedAfter(plainA, plainB);
+}
+
+function log(message) {
+  // eslint-disable-next-line no-console -- benchmark output
+  console.log(message);
+}
+
+const times = new Map();
+new Benchmark.Suite()
+  .add('Float64Array only  before', () => typedBefore(typedA, typedB), {
+    minSamples: 30,
+  })
+  .add('Float64Array only  after', () => typedAfter(typedA, typedB), {
+    minSamples: 30,
+  })
+  .add('Array only         before', () => plainBefore(plainA, plainB), {
+    minSamples: 30,
+  })
+  .add('Array only         after', () => plainAfter(plainA, plainB), {
+    minSamples: 30,
+  })
+  .add('both kinds         before', () => mixedBefore(typedA, typedB), {
+    minSamples: 30,
+  })
+  .add('both kinds         after', () => mixedAfter(typedA, typedB), {
+    minSamples: 30,
+  })
+  .on('cycle', (event) => {
+    const { name, hz, stats } = event.target;
+    times.set(name, 1e9 / hz / LENGTH);
+    log(
+      `${name.padEnd(28)}${(1e9 / hz / LENGTH).toFixed(2).padStart(7)} ns/element  ±${stats.rme.toFixed(1)}%`,
+    );
+  })
+  .on('complete', () => {
+    log('\ngain from caching the element reads:');
+    for (const [label, before, after] of [
+      [
+        'Float64Array only',
+        'Float64Array only  before',
+        'Float64Array only  after',
+      ],
+      [
+        'Array only       ',
+        'Array only         before',
+        'Array only         after',
+      ],
+      [
+        'both kinds       ',
+        'both kinds         before',
+        'both kinds         after',
+      ],
+    ]) {
+      log(`  ${label}  ${(times.get(before) / times.get(after)).toFixed(2)}x`);
+    }
+  })
+  .run({ async: false });
diff --git a/benchmark/equivalence.js b/benchmark/equivalence.js
new file mode 100644
index 0000000..000ec69
--- /dev/null
+++ b/benchmark/equivalence.js
@@ -0,0 +1,137 @@
+/*
+ * Differential test: is `x ** 1.5` bit-identical to `x * Math.sqrt(x)`?
+ *
+ * Mathematically x^1.5 = x * sqrt(x), but the two expressions round
+ * differently: the multiplication form rounds twice (sqrt, then *), while `**`
+ * calls the engine's pow. Neither is guaranteed correctly rounded by the spec
+ * (ECMA-262 leaves Math.pow implementation-defined), so this measures the
+ * actual disagreement instead of assuming it away.
+ *
+ * Run with: node benchmark/equivalence.js
+ */
+
+function log(message) {
+  // eslint-disable-next-line no-console -- benchmark output
+  console.log(message);
+}
+
+const buffer = new DataView(new ArrayBuffer(8));
+
+/**
+ * Distance in representable doubles (ULPs) between two finite values.
+ * @param x
+ * @param y
+ */
+function ulpDistance(x, y) {
+  if (x === y) return 0n;
+  if (Number.isNaN(x) || Number.isNaN(y)) return null;
+  return ordinal(y) - ordinal(x);
+}
+
+/**
+ * Maps a double onto a monotonic signed integer, so subtraction counts ULPs.
+ * @param value
+ */
+function ordinal(value) {
+  buffer.setFloat64(0, value);
+  const bits = buffer.getBigUint64(0);
+  return bits & 0x8000000000000000n
+    ? -(bits & 0x7fffffffffffffffn)
+    : BigInt(bits);
+}
+
+log('--- special values ---');
+for (const x of [
+  0,
+  -0,
+  1,
+  -1,
+  Infinity,
+  -Infinity,
+  Number.NaN,
+  Number.MIN_VALUE,
+  Number.MAX_VALUE,
+  Number.EPSILON,
+]) {
+  const pow = x ** 1.5;
+  const mul = x * Math.sqrt(x);
+  const same = Object.is(pow, mul);
+  log(
+    `  x = ${String(x).padEnd(24)} x**1.5 = ${String(pow).padEnd(24)} x*sqrt(x) = ${String(mul).padEnd(24)} ${same ? 'same' : 'DIFFERENT'}`,
+  );
+}
+
+log('\n--- random sweep across magnitudes ---');
+let state = 88172645463325252n;
+const mask = (1n << 64n) - 1n;
+function nextBits() {
+  state ^= (state << 13n) & mask;
+  state ^= state >> 7n;
+  state ^= (state << 17n) & mask;
+  return state;
+}
+
+for (const [label, low, high] of [
+  ['[1e-300, 1e-200)', 1e-300, 1e-200],
+  ['[1e-10, 1e-5)', 1e-10, 1e-5],
+  ['[0.1, 2)', 0.1, 2],
+  ['[1, 1000)', 1, 1000],
+  ['[1e100, 1e200)', 1e100, 1e200],
+]) {
+  const SAMPLES = 2_000_000;
+  const logLow = Math.log(low);
+  const logSpan = Math.log(high) - logLow;
+  let identical = 0;
+  let maxUlp = 0n;
+  let maxRelative = 0;
+  for (let i = 0; i < SAMPLES; i++) {
+    const unit = Number(nextBits() >> 11n) / 2 ** 53;
+    const x = Math.exp(logLow + unit * logSpan);
+    const pow = x ** 1.5;
+    const mul = x * Math.sqrt(x);
+    if (pow === mul) {
+      identical++;
+      continue;
+    }
+    const ulps = ulpDistance(pow, mul);
+    const absolute = ulps < 0n ? -ulps : ulps;
+    if (absolute > maxUlp) maxUlp = absolute;
+    const relative = Math.abs(pow - mul) / Math.abs(pow);
+    if (relative > maxRelative) maxRelative = relative;
+  }
+  const percent = ((identical / SAMPLES) * 100).toFixed(4);
+  log(
+    `  ${label.padEnd(18)} identical: ${percent}%  max |diff|: ${maxUlp} ulp  max relative: ${maxRelative.toExponential(3)}`,
+  );
+}
+
+log('\n--- exhaustive mantissa scan in [1, 2) ---');
+{
+  // Walk consecutive doubles from 1.0 upward: every value tested is adjacent to
+  // the previous one, so this is exhaustive over the scanned prefix.
+  const STEPS = 5_000_000;
+  let identical = 0;
+  let maxUlp = 0n;
+  let x = 1;
+  for (let i = 0; i < STEPS; i++) {
+    const pow = x ** 1.5;
+    const mul = x * Math.sqrt(x);
+    if (pow === mul) {
+      identical++;
+    } else {
+      const ulps = ulpDistance(pow, mul);
+      const absolute = ulps < 0n ? -ulps : ulps;
+      if (absolute > maxUlp) maxUlp = absolute;
+    }
+    x = nextUp(x);
+  }
+  log(
+    `  ${STEPS} consecutive doubles from 1.0: identical ${((identical / STEPS) * 100).toFixed(4)}%  max |diff|: ${maxUlp} ulp`,
+  );
+}
+
+function nextUp(value) {
+  buffer.setFloat64(0, value);
+  buffer.setBigUint64(0, buffer.getBigUint64(0) + 1n);
+  return buffer.getFloat64(0);
+}
diff --git a/benchmark/exponentiation.js b/benchmark/exponentiation.js
new file mode 100644
index 0000000..0f3b542
--- /dev/null
+++ b/benchmark/exponentiation.js
@@ -0,0 +1,184 @@
+/*
+ * The exponentiation operator is the reason minkowski and kumarJohnson were
+ * slow. This isolates which forms of `**` an engine specializes and which fall
+ * back to a generic pow call.
+ *
+ * Run with: node benchmark/exponentiation.js  (or: bun benchmark/exponentiation.js)
+ */
+import Benchmark from 'benchmark';
+
+const LENGTH = 10000;
+const values = new Float64Array(LENGTH);
+let state = 42;
+for (let i = 0; i < LENGTH; i++) {
+  state = (state * 1103515245 + 12345) % 2147483648;
+  values[i] = 0.1 + (state / 2147483648) * 2;
+}
+
+// A variable exponent, opaque to the compiler.
+function exponent(value) {
+  return values.length > 0 ? value : 0;
+}
+const variableOne = exponent(1);
+const variableTwo = exponent(2);
+const variableOneAndAHalf = exponent(1.5);
+const variableHalf = exponent(0.5);
+
+// `power` reaches powerCosine as a parameter with a default, not as a constant.
+function poweredSum(input, power) {
+  let s = 0;
+  for (let i = 0; i < input.length; i++) s += input[i] ** power;
+  return s;
+}
+function sqrtSum(input) {
+  let s = 0;
+  for (let i = 0; i < input.length; i++) s += Math.sqrt(input[i]);
+  return s;
+}
+
+const CASES = [
+  [
+    'x ** 2   (literal)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** 2;
+      return s;
+    },
+  ],
+  [
+    'x * x',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) {
+        const v = values[i];
+        s += v * v;
+      }
+      return s;
+    },
+  ],
+  [
+    'x ** p   (p = 2)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** variableTwo;
+      return s;
+    },
+  ],
+  [
+    'x ** 1   (literal)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** 1;
+      return s;
+    },
+  ],
+  [
+    'x ** p   (p = 1)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** variableOne;
+      return s;
+    },
+  ],
+  [
+    'x        (identity)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i];
+      return s;
+    },
+  ],
+  [
+    'x ** 0.5 (literal)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** 0.5;
+      return s;
+    },
+  ],
+  [
+    'Math.sqrt(x)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += Math.sqrt(values[i]);
+      return s;
+    },
+  ],
+  [
+    'x ** 1.5 (literal)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** 1.5;
+      return s;
+    },
+  ],
+  [
+    'x ** p   (p = 1.5)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** variableOneAndAHalf;
+      return s;
+    },
+  ],
+  [
+    'x * Math.sqrt(x)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) {
+        const v = values[i];
+        s += v * Math.sqrt(v);
+      }
+      return s;
+    },
+  ],
+  [
+    'x ** p   (p = 0.5)',
+    () => {
+      let s = 0;
+      for (let i = 0; i < LENGTH; i++) s += values[i] ** variableHalf;
+      return s;
+    },
+  ],
+  ['x ** p   (p = 0.5 param)', () => poweredSum(values, 0.5)],
+  ['Math.sqrt(x)  (in fn)', () => sqrtSum(values)],
+];
+
+function log(message) {
+  // eslint-disable-next-line no-console -- benchmark output
+  console.log(message);
+}
+
+const suite = new Benchmark.Suite();
+for (const [name, callback] of CASES) {
+  suite.add(name, callback, { maxTime: 2 });
+}
+
+const timings = new Map();
+suite
+  .on('cycle', (event) => {
+    const { name, hz, stats } = event.target;
+    const nanosecondsPerElement = 1e9 / hz / LENGTH;
+    timings.set(name, nanosecondsPerElement);
+    log(
+      `${name.padEnd(22)}${nanosecondsPerElement.toFixed(3).padStart(8)} ns/element  ±${stats.rme.toFixed(2)}%`,
+    );
+  })
+  .on('complete', () => {
+    log('\nrelative cost (vs the multiplication-based equivalent):');
+    for (const [slow, fast] of [
+      ['x ** 2   (literal)', 'x * x'],
+      ['x ** p   (p = 2)', 'x * x'],
+      ['x ** 1   (literal)', 'x        (identity)'],
+      ['x ** p   (p = 1)', 'x        (identity)'],
+      ['x ** 0.5 (literal)', 'Math.sqrt(x)'],
+      ['x ** p   (p = 0.5)', 'Math.sqrt(x)'],
+      ['x ** p   (p = 0.5 param)', 'Math.sqrt(x)  (in fn)'],
+      ['x ** 1.5 (literal)', 'x * Math.sqrt(x)'],
+      ['x ** p   (p = 1.5)', 'x * Math.sqrt(x)'],
+    ]) {
+      log(
+        `  ${slow.padEnd(22)}${(timings.get(slow) / timings.get(fast)).toFixed(2).padStart(7)}x the cost of  ${fast}`,
+      );
+    }
+  })
+  .run({ async: false });
diff --git a/package.json b/package.json
index b458f0c..4752b5d 100644
--- a/package.json
+++ b/package.json
@@ -49,7 +49,6 @@
   "homepage": "https://github.com/mljs/distance",
   "dependencies": {
     "cheminfo-types": "^1.15.0",
-    "ml-array-mean": "^2.0.0",
     "ml-distance-euclidean": "^3.0.1",
     "ml-tree-similarity": "^1.0.0"
   },
@@ -57,6 +56,7 @@
     "@types/node": "^26.1.1",
     "@vitest/coverage-v8": "^4.1.10",
     "@zakodium/tsconfig": "^1.0.5",
+    "benchmark": "^2.1.4",
     "eslint": "^9.39.5",
     "eslint-config-cheminfo-typescript": "^22.1.0",
     "prettier": "^3.9.6",
diff --git a/src/distances/__tests__/kumarJohnson.test.ts b/src/distances/__tests__/kumarJohnson.test.ts
index c9bb5de..9dd489e 100644
--- a/src/distances/__tests__/kumarJohnson.test.ts
+++ b/src/distances/__tests__/kumarJohnson.test.ts
@@ -6,5 +6,5 @@ const v1 = [0.2, 0.4, 0.3, 0.1];
 const v2 = [0.3, 0.2, 0.3, 0.2];
 
 test('should be correct', () => {
-  expect(distance.kumarJohnson(v1, v2)).toBe(0.5623488044808911);
+  expect(distance.kumarJohnson(v1, v2)).toBeCloseTo(0.5623488044808911, 15);
 });
diff --git a/src/distances/kumarJohnson.ts b/src/distances/kumarJohnson.ts
index 208309f..8435c1a 100644
--- a/src/distances/kumarJohnson.ts
+++ b/src/distances/kumarJohnson.ts
@@ -8,7 +8,11 @@ import type { NumberArray } from 'cheminfo-types';
 export function kumarJohnson(a: NumberArray, b: NumberArray): number {
   let ans = 0;
   for (let i = 0; i < a.length; i++) {
-    ans += (a[i] * a[i] - b[i] * b[i]) ** 2 / (2 * (a[i] * b[i]) ** 1.5);
+    const numerator = a[i] * a[i] - b[i] * b[i];
+    // `prod * Math.sqrt(prod)` is ~6x faster than `prod ** 1.5`, which no
+    // engine specializes; it costs at most 1 ulp of accuracy
+    const prod = a[i] * b[i];
+    ans += (numerator * numerator) / (2 * prod * Math.sqrt(prod));
   }
   return ans;
 }
diff --git a/src/distances/minkowski.ts b/src/distances/minkowski.ts
index e0639c4..e241504 100644
--- a/src/distances/minkowski.ts
+++ b/src/distances/minkowski.ts
@@ -8,6 +8,21 @@ import type { NumberArray } from 'cheminfo-types';
  */
 export function minkowski(a: NumberArray, b: NumberArray, p: number) {
   let d = 0;
+  // `x ** p` is far slower than the equivalent multiplication: ~9x for p = 1
+  // and ~1.5x for p = 2, the two orders that are used in practice.
+  if (p === 1) {
+    for (let i = 0; i < a.length; i++) {
+      d += Math.abs(a[i] - b[i]);
+    }
+    return d;
+  }
+  if (p === 2) {
+    for (let i = 0; i < a.length; i++) {
+      const diff = a[i] - b[i];
+      d += diff * diff;
+    }
+    return Math.sqrt(d);
+  }
   for (let i = 0; i < a.length; i++) {
     d += Math.abs(a[i] - b[i]) ** p;
   }
diff --git a/src/similarities/pearson.ts b/src/similarities/pearson.ts
index 59d7e4c..442d3e0 100644
--- a/src/similarities/pearson.ts
+++ b/src/similarities/pearson.ts
@@ -1,7 +1,4 @@
 import type { NumberArray } from 'cheminfo-types';
-import mean from 'ml-array-mean';
-
-import { cosine } from './cosine.ts';
 
 /**
  * Returns the Pearson correlation between vectors a and b, i.e. the cosine
@@ -10,15 +7,25 @@ import { cosine } from './cosine.ts';
  * @param b - second vector
  */
 export function pearson(a: NumberArray, b: NumberArray): number {
-  const avgA = mean(a);
-  const avgB = mean(b);
-
-  const newA = new Array(a.length);
-  const newB = new Array(b.length);
-  for (let i = 0; i < newA.length; i++) {
-    newA[i] = a[i] - avgA;
-    newB[i] = b[i] - avgB;
+  const length = a.length;
+  let sumA = 0;
+  let sumB = 0;
+  for (let i = 0; i < length; i++) {
+    sumA += a[i];
+    sumB += b[i];
   }
+  const avgA = sumA / length;
+  const avgB = sumB / length;
 
-  return cosine(newA, newB);
+  let p = 0;
+  let p2 = 0;
+  let q2 = 0;
+  for (let i = 0; i < length; i++) {
+    const centredA = a[i] - avgA;
+    const centredB = b[i] - avgB;
+    p += centredA * centredB;
+    p2 += centredA * centredA;
+    q2 += centredB * centredB;
+  }
+  return p / (Math.sqrt(p2) * Math.sqrt(q2));
 }