diff --git a/BENCHMARKS.md b/BENCHMARKS.md
new file mode 100644
index 0000000..3fe1121
--- /dev/null
+++ b/BENCHMARKS.md
@@ -0,0 +1,186 @@
+# Radix sort benchmarks
+
+Measured on Apple M3, darwin/arm64, Go 1.27.1, on 2026-09-22.
+
+The baseline is `radixSort` from commit `e882f60`. Both implementations were
+measured in the same test binary using a temporary side-by-side benchmark,
+100 ms per sample and six samples per case. The table reports medians.
+Input generation was excluded; each timed iteration copied the same deterministic
+input into a reusable slice before sorting, matching the current
+`BenchmarkSortDistributions` methodology. The temporary copy of the old sort
+was removed after measurement. These are local measurements, not performance
+guarantees for other CPUs or distributions.
+
+| Elements | Distribution | Before (µs) | After (µs) | Speedup |
+| ---: | --- | ---: | ---: | ---: |
+| 1,000 | duplicates | 8.84 | 7.25 | 1.22× |
+| 1,000 | random | 8.84 | 7.11 | 1.24× |
+| 1,000 | wide | 9.82 | 4.74 | 2.07× |
+| 1,000 | sorted | 0.51 | 0.26 | 1.97× |
+| 1,000 | response_time | 26.25 | 7.83 | 3.35× |
+| 10,000 | duplicates | 82.10 | 54.20 | 1.51× |
+| 10,000 | random | 123.31 | 82.97 | 1.49× |
+| 10,000 | wide | 94.58 | 65.91 | 1.43× |
+| 10,000 | sorted | 5.83 | 3.38 | 1.72× |
+| 10,000 | response_time | 277.15 | 72.95 | 3.80× |
+| 100,000 | duplicates | 828.36 | 499.02 | 1.66× |
+| 100,000 | random | 1,197.41 | 832.73 | 1.44× |
+| 100,000 | wide | 900.32 | 634.49 | 1.42× |
+| 100,000 | sorted | 52.58 | 26.94 | 1.95× |
+| 100,000 | response_time | 2,780.05 | 657.68 | 4.23× |
+| 1,000,000 | duplicates | 8,120.89 | 4,872.78 | 1.67× |
+| 1,000,000 | random | 12,314.01 | 8,304.43 | 1.48× |
+| 1,000,000 | wide | 9,431.60 | 6,346.96 | 1.49× |
+| 1,000,000 | sorted | 599.04 | 334.09 | 1.79× |
+| 1,000,000 | response_time | 27,772.10 | 6,137.81 | 4.52× |
+
+Sorting uses one scratch allocation (802,816 B/op for 100,000 elements), as
+before. Sorted inputs use zero allocations. The SIMD scan itself uses zero
+allocations. Go may also grow its goroutine stack on the first call; the LSD
+histograms occupy 96 KiB on 64-bit systems.
+
+## `slices.Sort` vs `radixSort`
+
+Measured on the same host and harness as above, after adding a reverse-order
+fast path to `radixSort`. The table compares `slices.Sort` (Go 1.27's pdqsort)
+with the package's `radixSort`. Each timed iteration copied the same
+deterministic input into a reusable slice before sorting, so the copy cost is
+included for both algorithms. Values are medians of six 100 ms samples.
+
+| Elements | Distribution | `slices.Sort` (µs) | `radixSort` (µs) | Speedup |
+| ---: | --- | ---: | ---: | ---: |
+| 512 | duplicates | 6.63 | 4.35 | 1.52× |
+| 512 | equal | 0.41 | 0.13 | 3.15× |
+| 512 | random | 5.72 | 3.29 | 1.74× |
+| 512 | response_time | 5.59 | 4.86 | 1.15× |
+| 512 | reverse | 0.54 | 0.27 | 2.00× |
+| 512 | reverse_duplicates | 2.07 | 0.27 | 7.67× |
+| 512 | sorted | 0.41 | 0.13 | 3.15× |
+| 512 | wide | 5.52 | 2.61 | 2.12× |
+| 1,000 | duplicates | 8.80 | 7.24 | 1.22× |
+| 1,000 | equal | 0.79 | 0.26 | 3.04× |
+| 1,000 | random | 12.05 | 7.13 | 1.69× |
+| 1,000 | response_time | 11.93 | 7.87 | 1.52× |
+| 1,000 | reverse | 1.02 | 0.52 | 1.96× |
+| 1,000 | reverse_duplicates | 3.12 | 0.52 | 6.00× |
+| 1,000 | sorted | 0.79 | 0.26 | 3.04× |
+| 1,000 | wide | 12.06 | 4.82 | 2.50× |
+| 2,048 | duplicates | 19.89 | 13.14 | 1.51× |
+| 2,048 | equal | 1.61 | 0.53 | 3.04× |
+| 2,048 | random | 27.30 | 22.36 | 1.22× |
+| 2,048 | response_time | 26.56 | 15.36 | 1.73× |
+| 2,048 | reverse | 2.06 | 1.05 | 1.96× |
+| 2,048 | reverse_duplicates | 6.91 | 1.05 | 6.58× |
+| 2,048 | sorted | 1.62 | 0.52 | 3.12× |
+| 2,048 | wide | 27.20 | 18.43 | 1.48× |
+| 4,096 | duplicates | 36.15 | 24.99 | 1.45× |
+| 4,096 | equal | 3.27 | 1.06 | 3.08× |
+| 4,096 | random | 60.78 | 38.51 | 1.58× |
+| 4,096 | response_time | 57.36 | 29.53 | 1.94× |
+| 4,096 | reverse | 4.62 | 2.12 | 2.18× |
+| 4,096 | reverse_duplicates | 14.53 | 2.10 | 6.92× |
+| 4,096 | sorted | 3.35 | 1.06 | 3.16× |
+| 4,096 | wide | 62.34 | 29.63 | 2.10× |
+| 10,000 | duplicates | 93.18 | 54.81 | 1.70× |
+| 10,000 | equal | 8.79 | 3.60 | 2.44× |
+| 10,000 | random | 243.34 | 83.06 | 2.93× |
+| 10,000 | response_time | 146.77 | 75.65 | 1.94× |
+| 10,000 | reverse | 10.60 | 5.77 | 1.84× |
+| 10,000 | reverse_duplicates | 32.21 | 6.03 | 5.34× |
+| 10,000 | sorted | 9.15 | 3.55 | 2.58× |
+| 10,000 | wide | 260.16 | 64.66 | 4.02× |
+| 100,000 | duplicates | 1,465.99 | 509.70 | 2.88× |
+| 100,000 | equal | 81.71 | 28.06 | 2.91× |
+| 100,000 | random | 6,053.73 | 830.42 | 7.29× |
+| 100,000 | response_time | 3,276.17 | 653.52 | 5.01× |
+| 100,000 | reverse | 103.17 | 57.21 | 1.80× |
+| 100,000 | reverse_duplicates | 305.79 | 56.72 | 5.39× |
+| 100,000 | sorted | 81.46 | 27.63 | 2.95× |
+| 100,000 | wide | 6,050.62 | 628.90 | 9.62× |
+| 1,000,000 | duplicates | 15,794.87 | 4,853.88 | 3.25× |
+| 1,000,000 | equal | 843.87 | 313.70 | 2.69× |
+| 1,000,000 | random | 74,769.90 | 8,246.56 | 9.07× |
+| 1,000,000 | response_time | 32,746.58 | 5,941.56 | 5.51× |
+| 1,000,000 | reverse | 1,071.76 | 608.60 | 1.76× |
+| 1,000,000 | reverse_duplicates | 3,118.17 | 610.98 | 5.10× |
+| 1,000,000 | sorted | 843.92 | 313.21 | 2.69× |
+| 1,000,000 | wide | 73,764.41 | 6,209.09 | 11.88× |
+
+`slices.Sort` is in-place and allocation-free. `radixSort` uses one scratch
+allocation for non-trivial inputs, but is substantially faster on random,
+wide, duplicates, response-time, reverse, and reverse-duplicate distributions
+once the element count exceeds a few thousand. For already-sorted or equal
+inputs the radix scan overhead makes it slower than pdqsort's fast paths at
+the smallest sizes, though it still wins by roughly 2.7–3× at 1,000,000
+elements.
+
+## SIMD scan in isolation
+
+These measurements exclude the reset copy and sorting. Medians of three
+100 ms samples, 100,000 elements:
+
+| Distribution | Portable Go (µs) | ARM64 NEON (µs) | Speedup |
+| --- | ---: | ---: | ---: |
+| wide | 54.27 | 16.93 | 3.21× |
+| sorted | 42.23 | 16.93 | 2.49× |
+
+## Reproducing current measurements
+
+```sh
+go test -run='^$' -bench=BenchmarkSortDistributions -benchmem -benchtime=100ms -count=6 .
+go test -run='^$' -bench=BenchmarkRadixScan -benchmem -benchtime=100ms -count=3 .
+go test -run='^$' -bench=BenchmarkSampdoAll -benchmem -benchtime=100ms -count=3 .
+```
+
+When comparing revisions, use the same benchmark harness in both versions.
+The previous `BenchmarkSortDistributions` excluded the reset copy and regenerated
+input on every iteration, so its historical timings are not directly comparable
+with the current harness.
+
+## Validation
+
+- `go test ./...`, `go test -race ./...`, and `go test -tags=purego ./...`.
+- `golangci-lint run --timeout 5m ./...`.
+- 30 seconds of `FuzzRadixSort`: 334,306 executions without a mismatch.
+- Differential scan checks against portable Go, every adjacent inversion around
+  SIMD boundaries, every non-sign bit, unaligned slices, and scalar tails.
+- Sorting checks against `slices.Sort`, including subnormals, positive infinity,
+  odd/even radix pass counts, both radix algorithms, and in-place aliasing.
+- Linux/amd64 test binary cross-build and `go vet`; Linux/386 cross-build.
+
+Only ARM64 was executed locally. The amd64 SSE2 implementation was cross-built
+and vetted, but could not be executed on this host. The existing Linux/amd64 CI
+test job exercises it. Other architectures and race builds use the portable scan.
+
+## Descending-input follow-up
+
+The follow-up to `7e6b298` detects a fully non-increasing input before the radix
+scan and reverses it in place. The entire input is checked; descending endpoints
+alone never cause an early return. This also supports descending inputs with
+duplicate values, positive infinity, subnormal values, and zero.
+
+Same Apple M3 / Go 1.27.1 environment, medians of three 100 ms samples, including
+the same reset copy for both algorithms:
+
+| Elements | slices.Sort (µs) | radixSort (µs) | Speedup vs slices.Sort |
+| ---: | ---: | ---: | ---: |
+| 512 | 0.53 | 0.27 | 1.98× |
+| 1,000 | 1.01 | 0.51 | 1.97× |
+| 2,048 | 2.15 | 1.05 | 2.04× |
+| 4,096 | 4.07 | 2.08 | 1.96× |
+| 10,000 | 10.56 | 5.73 | 1.84× |
+| 100,000 | 104.93 | 57.15 | 1.84× |
+| 1,000,000 | 1071.92 | 607.19 | 1.77× |
+
+All descending cases allocate 0 B/op and 0 allocs/op. Before this follow-up,
+100,000 descending elements took 593.38 µs and allocated 802,816 bytes; afterward
+they take 57.15 µs with no allocation. Other existing distributions were also
+remeasured at 1,000 through 1,000,000 elements, with median changes within ±5%.
+
+Tests cover cutoff boundaries, duplicates, extreme values, in-place aliasing,
+zero allocations on every reset-and-sort iteration, and every possible single
+ascending pair within a 513-element descending input.
+
+```sh
+go test -run='^$' -bench='BenchmarkSortDistributions/(slices|radix)/.*/reverse$' -benchmem -benchtime=100ms -count=3 .
+```
diff --git a/Makefile b/Makefile
index 60d0933..3a9420d 100644
--- a/Makefile
+++ b/Makefile
@@ -6,6 +6,7 @@ all: check lint bench
 check: *.go
 	go test -v ./...
 	go test -race ./...
+	go test -tags=purego ./...
 
 lint:
 	golangci-lint run --timeout 5m ./...
diff --git a/README.md b/README.md
index 1f1c0c7..745cfec 100644
--- a/README.md
+++ b/README.md
@@ -7,11 +7,41 @@
 ## Features / 特徴
 
 - **Sort-based aggregation / ソートベースの集計**: Compute several statistics from a single sort. サンプルを 1 回ソートするだけで複数の統計量を取得できます。
-- **Radix sort for non-negative numbers / 非負数向け基数ソート**: Falls back to `slices.Sort` only when negative or `-0` values are present. 負数や `-0` がある場合のみ `slices.Sort` にフォールバックします。
+- **Radix sort for non-negative numbers / 非負数向け基数ソート**: Uses `slices.Sort` for small inputs, negative values, or `-0`. 小さい入力、負数や `-0` を含む入力では `slices.Sort` を使用します。
 - **NaN rejection / NaN を排除**: `Append` rejects `NaN` so results stay stable. `Append` は `NaN` を弾き、集計結果を安定させます。
 - **Immutable sorted view / イミュータブルなソート結果**: Appending after `Sorted()` does not mutate the returned `*Sorted`. `Sorted()` 後に追加しても返された `*Sorted` は影響を受けません。
 - **Infinity support / 無限大対応**: Handles `+Inf` and `-Inf` in percentile calculations. `+Inf`、`-Inf` を含むデータのパーセンタイル計算にも対応しています。
 
+### Radix sort implementation / 基数ソートの実装
+
+The key scan checks ordering and finds varying key bits using ARM64 NEON or
+amd64 SSE2. Other architectures, `-tags=purego`, and race builds use the portable
+Go implementation. No additional dependencies or experimental build flags are required.
+
+事前走査では ARM64 NEON / amd64 SSE2 を使い、整列状態と値が異なるビットをまとめて調べます。
+その他のアーキテクチャ、`-tags=purego`、race ビルドでは純 Go 実装を使います。
+追加の依存パッケージや実験的なビルドフラグは不要です。
+
+Large inputs normally use six 11-bit LSD passes, with all histograms collected in
+one traversal and constant digits skipped. Small inputs and sampled distributions
+with many repeated full-width keys use an MSD partitioning path, which can stop
+at equal buckets. Sampling only selects the algorithm; all keys are sorted exactly.
+Both paths reuse one scratch slice and write the result into the original input.
+Already ordered inputs of at least 512 elements return without a scratch allocation.
+Descending inputs (including equal neighbors) are confirmed in a linear scan and
+reversed in place, avoiding radix passes and scratch allocation. The endpoint
+check skips this extra scan for ascending and equal inputs, and unordered inputs
+stop it at the first ascending pair.
+
+大きい入力は通常、全桁の度数を一度に集計する最大 6 パスの 11 ビット LSD 方式で処理し、
+値が変わらない桁は省略します。小さい入力や、広いビット範囲にわたる値の重複が多い入力では
+MSD 方式を使い、同じ値だけの区画の処理を打ち切ります。サンプリングは方式選択のみに使い、
+ソート結果は全要素について厳密です。作業配列は 1 個を再利用し、結果は元の配列に書き戻します。
+512 要素以上の整列済み入力では作業配列を確保しません。
+降順の入力は、同じ値の連続も許して全体を確認し、その場で反転します。
+基数ソートや作業配列の確保は不要です。昇順・全要素同値の入力では端点の比較で
+降順チェックを省略し、未整列の入力では昇順の隣接ペアが見つかった時点で打ち切ります。
+
 ## Installation / インストール
 
 ```bash
@@ -124,11 +154,18 @@ Run benchmarks with:
 
 ```bash
 go test -bench BenchmarkSampdoAll -benchmem -benchtime=100ms -count 10 -run '^$' .
+go test -bench BenchmarkSortDistributions -benchmem -benchtime=100ms -count 6 -run '^$' .
+go test -bench BenchmarkRadixScan -benchmem -benchtime=100ms -run '^$' .
 ```
 
+`BenchmarkSortDistributions` includes an identical input-reset copy for each algorithm;
+input generation is excluded. See [BENCHMARKS.md](BENCHMARKS.md) for before/after results.
+
+`BenchmarkSortDistributions` は各方式に同じ入力復元コピーを含め、入力生成は計測から除外します。
+変更前後の実測値は [BENCHMARKS.md](BENCHMARKS.md) を参照してください。
+
 ## License / ライセンス
 
 This project is licensed under the [MIT License](LICENSE).
 
 本プロジェクトは [MIT ライセンス](LICENSE) の下で提供されています。
-
diff --git a/benchmark_test.go b/benchmark_test.go
index 81738fb..5c575bc 100644
--- a/benchmark_test.go
+++ b/benchmark_test.go
@@ -55,15 +55,18 @@ func BenchmarkSortDistributions(b *testing.B) {
 	}{
 		{"slices", aliasSlicesSort}, {"radix", radixSort},
 	} {
-		for _, n := range []int{1000, 10_000, 100_000} {
-			for _, distribution := range []string{"duplicates", "random", "wide", "sorted", "response_time"} {
+		for _, n := range []int{512, 1000, 2048, 4096, 10_000, 100_000, 1_000_000} {
+			for _, distribution := range []string{"duplicates", "random", "wide", "sorted", "reverse", "reverse_duplicates", "equal", "response_time"} {
 				{
 					b.Run(fmt.Sprintf("%s/%d/%s", algorithm.name, n, distribution), func(b *testing.B) {
+						input := radixInput(n, distribution)
+						points := make([]float64, n)
 						b.ReportAllocs()
+						b.ResetTimer()
+						// Include the same reset copy for both algorithms. Reuse
+						// storage so input generation/GC does not distort timing.
 						for b.Loop() {
-							b.StopTimer()
-							points := radixInput(n, distribution)
-							b.StartTimer()
+							copy(points, input)
 							algorithm.sort(points)
 						}
 					})
@@ -72,3 +75,23 @@ func BenchmarkSortDistributions(b *testing.B) {
 		}
 	}
 }
+
+func BenchmarkRadixScan(b *testing.B) {
+	for _, n := range []int{1000, 100_000} {
+		for _, distribution := range []string{"wide", "sorted"} {
+			points := radixInput(n, distribution)
+			for _, scan := range []struct {
+				name string
+				scan func([]float64) (uint64, bool)
+			}{{"native", radixScan}, {"generic", radixScanGeneric}} {
+				b.Run(fmt.Sprintf("%s/%d/%s", scan.name, n, distribution), func(b *testing.B) {
+					b.ReportAllocs()
+					b.SetBytes(int64(n) * 8)
+					for b.Loop() {
+						scan.scan(points)
+					}
+				})
+			}
+		}
+	}
+}
diff --git a/radix.go b/radix.go
new file mode 100644
index 0000000..e28743d
--- /dev/null
+++ b/radix.go
@@ -0,0 +1,163 @@
+package sampdo
+
+import (
+	"math"
+	"math/bits"
+	"slices"
+)
+
+// radixSort accepts non-negative float64 values, excluding NaN and negative zero.
+// Sampdo.Append and Sampdo.Sorted enforce these preconditions. The returned slice
+// always shares the input's backing array, including an odd number of radix passes.
+func radixSort(points []float64) []float64 {
+	if len(points) < 512 {
+		slices.Sort(points)
+		return points
+	}
+	// A descending input only needs a reversal. Check the endpoints first so
+	// ascending/equal inputs avoid this scan; unordered inputs exit it early.
+	if points[0] > points[len(points)-1] && radixDescending(points) {
+		slices.Reverse(points)
+		return points
+	}
+	varying, ordered := radixScan(points)
+	if ordered {
+		return points
+	}
+	scratch := make([]float64, len(points))
+	// Small partitions amortize the compact MSD histogram. Repeated full-width
+	// keys also benefit: equal buckets finish without visiting the remaining bits.
+	// Narrow keys need few LSD passes, so don't sample them.
+	if len(points) < 2048 || bits.Len64(varying)-bits.TrailingZeros64(varying) > 22 && radixRepeatedSample(points) {
+		radixSortMSDInto(points, scratch, varying)
+	} else {
+		radixSortLSD(points, scratch, varying)
+	}
+	return points
+}
+
+// Equal neighbors are allowed. Confirm the entire input before reversing it:
+// a descending prefix or descending endpoints alone are not sufficient.
+func radixDescending(points []float64) bool {
+	for i := 1; i < len(points); i++ {
+		if points[i] > points[i-1] {
+			return false
+		}
+	}
+	return true
+}
+
+// Sample 512 evenly spaced keys without allocating (len(points) >= 2048).
+// This only chooses
+// an algorithm; neither correctness nor the worst-case radix depth depends on
+// the sample being representative. Stop once at least 1/8 of the sample repeats.
+func radixRepeatedSample(points []float64) bool {
+	var seen [1024]uint64
+	repeats := 0
+	step := len(points) / 512
+	for i := range 512 {
+		key := math.Float64bits(points[i*step]) + 1 // zero is the empty slot sentinel
+		slot := key * 0x9e3779b97f4a7c15 >> 54
+		for seen[slot] != 0 && seen[slot] != key {
+			slot = (slot + 1) & 1023
+		}
+		if seen[slot] == key {
+			repeats++
+			if repeats == 64 {
+				return true
+			}
+		}
+		seen[slot] = key
+	}
+	return false
+}
+
+// Build all six histograms in one traversal. Eleven bits reduce the number of
+// scatter passes from eight to six while keeping each histogram at 16 KiB on
+// 64-bit systems. Histograms are independent of the order of previous passes.
+func radixSortLSD(points, scratch []float64, varying uint64) {
+	var counts [6][2048]int
+	for _, v := range points {
+		k := math.Float64bits(v)
+		counts[0][k&2047]++
+		counts[1][(k>>11)&2047]++
+		counts[2][(k>>22)&2047]++
+		counts[3][(k>>33)&2047]++
+		counts[4][(k>>44)&2047]++
+		counts[5][k>>55]++
+	}
+	src, dst := points, scratch
+	for pass := range 6 {
+		shift := uint(pass * 11)
+		if (varying>>shift)&2047 == 0 {
+			continue
+		}
+		c := &counts[pass]
+		sum := 0
+		for i, n := range c {
+			c[i] = sum
+			sum += n
+		}
+		for _, v := range src {
+			idx := (math.Float64bits(v) >> shift) & 2047
+			dst[c[idx]] = v
+			c[idx]++
+		}
+		src, dst = dst, src
+	}
+	if &src[0] != &points[0] {
+		copy(points, src)
+	}
+}
+
+// Partition from the most significant varying byte. Recursion skips shared
+// prefixes and stops on ordered/equal buckets; at most eight levels are needed.
+func radixSortMSDInto(points, scratch []float64, v uint64) {
+	shift := max(0, bits.Len64(v)-8)
+	count := radixPartition(points, scratch, uint(shift))
+	copy(points, scratch)
+	if shift == 0 {
+		return
+	}
+	start := 0
+	for _, end := range count {
+		if end-start > 128 {
+			vary, sorted := radixScan(points[start:end])
+			if !sorted {
+				radixSortMSDInto(points[start:end], scratch[start:end], vary)
+			}
+		} else if end-start > 1 {
+			slices.Sort(points[start:end])
+		}
+		start = end
+	}
+}
+
+// Four independent counters break the histogram dependency chain when many
+// adjacent keys share a byte. After scattering, offsets are bucket ends.
+func radixPartition(points, scratch []float64, shift uint) [256]int {
+	var counts [4][256]int
+	i := 0
+	for ; i+4 <= len(points); i += 4 {
+		counts[0][(math.Float64bits(points[i])>>shift)&255]++
+		counts[1][(math.Float64bits(points[i+1])>>shift)&255]++
+		counts[2][(math.Float64bits(points[i+2])>>shift)&255]++
+		counts[3][(math.Float64bits(points[i+3])>>shift)&255]++
+	}
+	for ; i < len(points); i++ {
+		counts[0][(math.Float64bits(points[i])>>shift)&255]++
+	}
+	count := &counts[0]
+	sum := 0
+	for i := range count {
+		n := counts[0][i] + counts[1][i] + counts[2][i] + counts[3][i]
+		count[i] = sum
+		sum += n
+	}
+	for _, val := range points {
+		idx := (math.Float64bits(val) >> shift) & 255
+		scratch[count[idx]] = val
+		count[idx]++
+	}
+	return *count
+}
diff --git a/radix_scan.go b/radix_scan.go
new file mode 100644
index 0000000..44dde29
--- /dev/null
+++ b/radix_scan.go
@@ -0,0 +1,21 @@
+package sampdo
+
+import "math"
+
+// radixScanGeneric is the portable reference for the SIMD scan. Non-negative,
+// non-NaN IEEE 754 bit patterns have the same ordering as their numeric values.
+func radixScanGeneric(points []float64) (varying uint64, ordered bool) {
+	if len(points) == 0 {
+		return 0, true
+	}
+	first := math.Float64bits(points[0])
+	previous := first
+	ordered = true
+	for _, v := range points {
+		key := math.Float64bits(v)
+		varying |= key ^ first
+		ordered = ordered && previous <= key
+		previous = key
+	}
+	return
+}
diff --git a/radix_scan_amd64.go b/radix_scan_amd64.go
new file mode 100644
index 0000000..98a7a68
--- /dev/null
+++ b/radix_scan_amd64.go
@@ -0,0 +1,8 @@
+//go:build !purego && !race
+
+package sampdo
+
+// radixScan uses SSE2, available on every amd64 CPU. No AVX feature check is needed.
+//
+//go:noescape
+func radixScan(points []float64) (varying uint64, ordered bool)
diff --git a/radix_scan_amd64.s b/radix_scan_amd64.s
new file mode 100644
index 0000000..53f66a1
--- /dev/null
+++ b/radix_scan_amd64.s
@@ -0,0 +1,85 @@
+//go:build !purego && !race
+
+#include "textflag.h"
+
+// The keys are non-negative, non-NaN float64s, so CMPPD has the same
+// ordering as an unsigned comparison of their bit patterns.
+TEXT ·radixScan(SB), NOSPLIT, $0-33
+	MOVQ points_base+0(FP), SI
+	MOVQ points_len+8(FP), CX
+	XORQ AX, AX
+	XORQ DX, DX
+	TESTQ CX, CX
+	JZ done
+	MOVQ (SI), R8
+	MOVQ R8, X0
+	PUNPCKLQDQ X0, X0
+	MOVOU X0, X1
+	PXOR X2, X2
+	PXOR X3, X3
+	CMPQ CX, $8
+	JB tail
+loop:
+	MOVOU (SI), X4
+	MOVOU 16(SI), X5
+	MOVOU 32(SI), X6
+	MOVOU 48(SI), X7
+	SHUFPD $1, X4, X1
+	MOVOU X4, X8
+	SHUFPD $1, X5, X8
+	MOVOU X5, X9
+	SHUFPD $1, X6, X9
+	MOVOU X6, X10
+	SHUFPD $1, X7, X10
+	CMPPD X4, X1, $6
+	CMPPD X5, X8, $6
+	CMPPD X6, X9, $6
+	CMPPD X7, X10, $6
+	POR X1, X8
+	POR X9, X10
+	POR X8, X10
+	POR X10, X3
+	MOVOU X7, X1
+	PXOR X0, X4
+	PXOR X0, X5
+	PXOR X0, X6
+	PXOR X0, X7
+	POR X4, X5
+	POR X6, X7
+	POR X5, X7
+	POR X7, X2
+	ADDQ $64, SI
+	SUBQ $8, CX
+	CMPQ CX, $8
+	JAE loop
+	MOVQ X2, AX
+	PSRLDQ $8, X2
+	MOVQ X2, R9
+	ORQ R9, AX
+	MOVMSKPD X3, DX
+	PSRLDQ $8, X1
+	MOVQ X1, R9
+	JMP tailCheck
+tail:
+	MOVQ R8, R9
+tailCheck:
+	TESTQ CX, CX
+	JZ done
+tailLoop:
+	MOVQ (SI), R10
+	MOVQ R10, R11
+	XORQ R8, R11
+	ORQ R11, AX
+	CMPQ R10, R9
+	SETCS R11
+	MOVBQZX R11, R11
+	ORQ R11, DX
+	MOVQ R10, R9
+	ADDQ $8, SI
+	DECQ CX
+	JNZ tailLoop
+done:
+	MOVQ AX, varying+24(FP)
+	TESTQ DX, DX
+	SETEQ ordered+32(FP)
+	RET
diff --git a/radix_scan_arm64.go b/radix_scan_arm64.go
new file mode 100644
index 0000000..0a35989
--- /dev/null
+++ b/radix_scan_arm64.go
@@ -0,0 +1,8 @@
+//go:build !purego && !race
+
+package sampdo
+
+// radixScan uses baseline ARM64 NEON instructions. It reads only within points.
+//
+//go:noescape
+func radixScan(points []float64) (varying uint64, ordered bool)
diff --git a/radix_scan_arm64.s b/radix_scan_arm64.s
new file mode 100644
index 0000000..8493a88
--- /dev/null
+++ b/radix_scan_arm64.s
@@ -0,0 +1,72 @@
+//go:build !purego && !race
+
+#include "textflag.h"
+
+// Reduce key differences and adjacent inversions, eight keys at a time.
+TEXT ·radixScan(SB), NOSPLIT, $0-33
+	MOVD points_base+0(FP), R0
+	MOVD points_len+8(FP), R1
+	MOVD $0, R2
+	MOVD $0, R3
+	CBZ R1, done
+	MOVD (R0), R4
+	VDUP R4, V0.D2
+	VMOV V0.B16, V1.B16
+	VEOR V2.B16, V2.B16, V2.B16
+	VEOR V3.B16, V3.B16, V3.B16
+	CMP $8, R1
+	BLT tail
+loop:
+	VLD1.P 64(R0), [V4.D2, V5.D2, V6.D2, V7.D2]
+	VEXT $8, V4.B16, V1.B16, V8.B16
+	VEXT $8, V5.B16, V4.B16, V9.B16
+	VEXT $8, V6.B16, V5.B16, V10.B16
+	VEXT $8, V7.B16, V6.B16, V11.B16
+	VCMHI V4.D2, V8.D2, V8.D2
+	VCMHI V5.D2, V9.D2, V9.D2
+	VCMHI V6.D2, V10.D2, V10.D2
+	VCMHI V7.D2, V11.D2, V11.D2
+	VORR V8.B16, V9.B16, V8.B16
+	VORR V10.B16, V11.B16, V10.B16
+	VORR V8.B16, V10.B16, V8.B16
+	VORR V8.B16, V3.B16, V3.B16
+	VMOV V7.B16, V1.B16
+	VEOR V0.B16, V4.B16, V4.B16
+	VEOR V0.B16, V5.B16, V5.B16
+	VEOR V0.B16, V6.B16, V6.B16
+	VEOR V0.B16, V7.B16, V7.B16
+	VORR V4.B16, V5.B16, V4.B16
+	VORR V6.B16, V7.B16, V6.B16
+	VORR V4.B16, V6.B16, V4.B16
+	VORR V4.B16, V2.B16, V2.B16
+	SUB $8, R1
+	CMP $8, R1
+	BGE loop
+	VMOV V2.D[0], R2
+	VMOV V2.D[1], R5
+	ORR R5, R2, R2
+	VMOV V3.D[0], R3
+	VMOV V3.D[1], R5
+	ORR R5, R3, R3
+	VMOV V1.D[1], R6
+	B tailCheck
+tail:
+	MOVD R4, R6
+tailCheck:
+	CBZ R1, done
+tailLoop:
+	MOVD.P 8(R0), R5
+	EOR R4, R5, R7
+	ORR R7, R2, R2
+	CMP R6, R5
+	CSET LO, R7
+	ORR R7, R3, R3
+	MOVD R5, R6
+	SUB $1, R1
+	CBNZ R1, tailLoop
+done:
+	MOVD R2, varying+24(FP)
+	CMP $0, R3
+	CSET EQ, R3
+	MOVB R3, ordered+32(FP)
+	RET
diff --git a/radix_scan_generic.go b/radix_scan_generic.go
new file mode 100644
index 0000000..ee4c997
--- /dev/null
+++ b/radix_scan_generic.go
@@ -0,0 +1,7 @@
+//go:build (!arm64 && !amd64) || purego || race
+
+package sampdo
+
+func radixScan(points []float64) (uint64, bool) {
+	return radixScanGeneric(points)
+}
diff --git a/radix_test.go b/radix_test.go
new file mode 100644
index 0000000..32e850f
--- /dev/null
+++ b/radix_test.go
@@ -0,0 +1,186 @@
+package sampdo
+
+import (
+	"encoding/binary"
+	"fmt"
+	"math"
+	"math/rand/v2"
+	"slices"
+	"testing"
+
+	"github.com/stretchr/testify/require"
+)
+
+func TestRadixScan(t *testing.T) {
+	for _, n := range []int{0, 1, 2, 3, 7, 8, 9, 15, 16, 17, 31, 32, 33, 127, 128, 129, 511, 512, 513, 2047, 2048, 2049} {
+		for _, d := range []string{"random", "wide", "sorted", "reverse", "equal", "response_time"} {
+			// Offset the slice to exercise unaligned SIMD loads as well as tails.
+			storage := make([]float64, n+2)
+			points := storage[1 : n+1]
+			copy(points, radixInput(n, d))
+			wantV, wantO := radixScanGeneric(points)
+			gotV, gotO := radixScan(points)
+			require.Equal(t, wantV, gotV, "varying %s/%d", d, n)
+			require.Equal(t, wantO, gotO, "ordered %s/%d", d, n)
+		}
+	}
+}
+
+func TestRadixScanEveryInversion(t *testing.T) {
+	// Check both SIMD lanes, vector boundaries, and the scalar tail, including
+	// an inversion that is the only difference from an otherwise ordered slice.
+	for n := 2; n <= 65; n++ {
+		for i := 0; i < n-1; i++ {
+			points := radixInput(n, "sorted")
+			points[i], points[i+1] = points[i+1], points[i]
+			want, _ := radixScanGeneric(points)
+			got, ordered := radixScan(points)
+			require.Equal(t, want, got, "n=%d inversion=%d", n, i)
+			require.False(t, ordered, "n=%d inversion=%d", n, i)
+		}
+	}
+	for bit := range 63 {
+		for i := range 33 {
+			points := make([]float64, 33)
+			points[i] = math.Float64frombits(uint64(1) << bit)
+			got, ordered := radixScan(points)
+			require.Equal(t, uint64(1)<<bit, got)
+			require.Equal(t, i == 32, ordered)
+		}
+	}
+}
+
+func checkRadixInPlace(t *testing.T, input []float64) {
+	t.Helper()
+	want := slices.Clone(input)
+	slices.Sort(want)
+	storage := make([]float64, len(input)+2)
+	storage[0], storage[len(storage)-1] = -123, -456
+	points := storage[1 : len(storage)-1]
+	copy(points, input)
+	got := radixSort(points)
+	require.Equal(t, want, got)
+	require.Equal(t, want, points, "input must be sorted even if return value is ignored")
+	if len(points) > 0 {
+		require.Same(t, &points[0], &got[0])
+	}
+	require.Equal(t, -123.0, storage[0])
+	require.Equal(t, -456.0, storage[len(storage)-1])
+}
+
+func TestRadixSortInPlace(t *testing.T) {
+	for _, n := range []int{511, 512, 513, 1023, 2047, 2048, 2049, 8193, 100000} {
+		for _, d := range []string{"duplicates", "random", "wide", "sorted", "reverse", "equal", "response_time", "nearly_sorted"} {
+			t.Run(fmt.Sprintf("%s/%d", d, n), func(t *testing.T) { checkRadixInPlace(t, radixInput(n, d)) })
+		}
+	}
+	// Force every LSD digit individually, and odd/even numbers of active digits.
+	// These also exercise all-subnormal inputs and skipping common middle digits.
+	for passes := 1; passes <= 6; passes++ {
+		for start := 0; start+passes <= 6; start++ {
+			points := make([]float64, 4097)
+			for i := range points {
+				var key uint64
+				for p := start; p < start+passes; p++ {
+					key |= uint64((i*17+p)%8) << uint(p*11)
+				}
+				points[i] = math.Float64frombits(key)
+			}
+			checkRadixInPlace(t, points)
+		}
+	}
+}
+
+func TestRadixSortPaths(t *testing.T) {
+	for _, d := range []string{"duplicates", "random", "wide", "response_time", "reverse", "equal"} {
+		original := radixInput(8193, d)
+		want := slices.Clone(original)
+		slices.Sort(want)
+		varying, _ := radixScanGeneric(original)
+		for _, sort := range []func([]float64, []float64, uint64){radixSortLSD, radixSortMSDInto} {
+			points := slices.Clone(original)
+			sort(points, make([]float64, len(points)), varying)
+			require.Equal(t, want, points, d)
+		}
+	}
+	require.True(t, radixRepeatedSample(radixInput(10000, "response_time")))
+	require.False(t, radixRepeatedSample(radixInput(10000, "wide")))
+}
+
+func TestRadixSortBitPatterns(t *testing.T) {
+	r := rand.New(rand.NewPCG(19, 37))
+	for trial := range 30 {
+		points := make([]float64, 2048+trial)
+		mask := r.Uint64() & 0x7fefffffffffffff
+		for i := range points {
+			points[i] = math.Float64frombits(r.Uint64() & mask)
+		}
+		points[len(points)/2] = math.Inf(1)
+		checkRadixInPlace(t, points)
+	}
+}
+
+func FuzzRadixSort(f *testing.F) {
+	f.Add([]byte{0, 0, 0, 0, 0, 0, 0xf0, 0x7f, 1, 0, 0, 0, 0, 0, 0, 0})
+	f.Add([]byte{1, 7, 2, 11, 3, 19, 5, 23, 8})
+	f.Fuzz(func(t *testing.T, data []byte) {
+		if len(data) == 0 {
+			return
+		}
+		n := min(4097, max(513, len(data)))
+		points := make([]float64, n)
+		for i := range points {
+			var raw [8]byte
+			for j := range raw {
+				raw[j] = data[(i*8+j)%len(data)]
+			}
+			key := binary.LittleEndian.Uint64(raw[:]) & 0x7fffffffffffffff
+			key = min(key, 0x7ff0000000000000)
+			points[i] = math.Float64frombits(key)
+		}
+		checkRadixInPlace(t, points)
+	})
+}
+
+func TestRadixSortOrderedAllocs(t *testing.T) {
+	for _, distribution := range []string{"sorted", "equal"} {
+		points := radixInput(8193, distribution)
+		require.Zero(t, testing.AllocsPerRun(10, func() { radixSort(points) }), distribution)
+	}
+}
+
+func TestRadixSortDescending(t *testing.T) {
+	for _, n := range []int{0, 1, 511, 512, 513, 2048, 2049, 100000} {
+		for _, distribution := range []string{"reverse", "reverse_duplicates"} {
+			points := radixInput(n, distribution)
+			checkRadixInPlace(t, points)
+			work := make([]float64, n)
+			require.Zero(t, testing.AllocsPerRun(10, func() {
+				copy(work, points)
+				radixSort(work)
+			}), "%s/%d", distribution, n)
+		}
+	}
+	points := make([]float64, 1025)
+	values := []float64{math.Inf(1), math.MaxFloat64, 1, math.SmallestNonzeroFloat64, 0}
+	for i := range points {
+		points[i] = values[i%len(values)]
+	}
+	slices.Sort(points)
+	slices.Reverse(points)
+	checkRadixInPlace(t, points)
+}
+
+func TestRadixSortAlmostDescending(t *testing.T) {
+	// Every possible ascending pair within a descending input, including a pair
+	// at the end, must prevent the reversal-only path from returning early.
+	for i := range 512 {
+		points := radixInput(513, "reverse")
+		points[i], points[i+1] = points[i+1], points[i]
+		checkRadixInPlace(t, points)
+	}
+	// Equal endpoints can also enclose an unordered interior.
+	points := radixInput(2049, "reverse")
+	points[0] = points[len(points)-1]
+	checkRadixInPlace(t, points)
+}
diff --git a/sorted.go b/sorted.go
index 7a1c376..4258724 100644
--- a/sorted.go
+++ b/sorted.go
@@ -2,9 +2,7 @@ package sampdo
 
 import (
 	"errors"
-
 	"math"
-	"slices"
 )
 
 var ErrNoPoints = errors.New("no points available")
@@ -16,57 +14,6 @@ type Sorted struct {
 	sumValid bool
 }
 
-// radixSort accepts non-negative float64 values, excluding NaN and negative zero.
-// Stats.Append and Stats.Sorted enforce these preconditions.
-// It sorts the provided slice in-place and returns it.
-func radixSort(points []float64) []float64 {
-	// Small inputs do not amortize radix sort's histogram and scratch space.
-	if len(points) < 512 {
-		slices.Sort(points)
-		return points
-	}
-
-	first := math.Float64bits(points[0])
-	var varying uint64
-	ordered := true
-	previous := first
-	for _, v := range points {
-		key := math.Float64bits(v)
-		varying |= key ^ first
-		ordered = ordered && previous <= key
-		previous = key
-	}
-	if ordered {
-		return points
-	}
-
-	sorted := make([]float64, len(points))
-	copy(sorted, points)
-	for shift := uint(0); shift < 64; shift += 8 {
-		// A byte shared by every key cannot change their order.
-		if (varying>>shift)&0xff == 0 {
-			continue
-		}
-		var count [256]int
-		for _, v := range sorted {
-			count[(math.Float64bits(v)>>shift)&0xff]++
-		}
-		sum := 0
-		for i, n := range count {
-			count[i] = sum
-			sum += n
-		}
-		for _, v := range sorted {
-			idx := (math.Float64bits(v) >> shift) & 0xff
-			points[count[idx]] = v
-			count[idx]++
-		}
-		sorted, points = points, sorted
-	}
-	copy(points, sorted)
-	return points
-}
-
 func (s *Sorted) Count() int {
 	return len(s.points)
 }
diff --git a/sorted_test.go b/sorted_test.go
index 05ddf9c..934275b 100644
--- a/sorted_test.go
+++ b/sorted_test.go
@@ -70,7 +70,7 @@ func TestSortedMedian(t *testing.T) {
 		t.Run(tt.name, func(t *testing.T) {
 			s := &Sorted{points: tt.pts}
 			median, err := s.Median()
-			if tt.pts == nil || len(tt.pts) == 0 {
+			if len(tt.pts) == 0 {
 				require.Error(t, err)
 				return
 			}
@@ -239,7 +239,7 @@ func TestRadixSortExtremeValues(t *testing.T) {
 	for i := range points {
 		points[i] = extremes[i%len(extremes)]
 	}
-	require.Equal(t, aliasSlicesSort(points), radixSort(points))
+	checkRadixInPlace(t, points)
 }
 
 func radixInput(n int, distribution string) []float64 {
@@ -253,6 +253,12 @@ func radixInput(n int, distribution string) []float64 {
 			points[i] = r.Float64() * 10000
 		case "wide":
 			points[i] = math.Float64frombits(r.Uint64() & 0x7fefffffffffffff)
+		case "reverse_duplicates":
+			points[i] = float64((n - 1 - i) / 8)
+		case "reverse":
+			points[i] = float64(n - i)
+		case "nearly_sorted":
+			points[i] = float64(i ^ 1)
 		case "sorted":
 			points[i] = float64(i)
 		case "equal":