diff --git a/Makefile b/Makefile
index 58c7bb1f..df38c0d6 100644
--- a/Makefile
+++ b/Makefile
@@ -133,7 +133,7 @@ indent:
 TESTS = $(wildcard tests/*/out/*.json)
 SPACE = $(NULL) $(NULL)
 
-test: tippecanoe tippecanoe-decode $(addsuffix .check,$(TESTS)) raw-tiles-test parallel-test pbf-test join-test enumerate-test decode-test join-filter-test unit json-tool-test allow-existing-test csv-test layer-json-test pmtiles-test decode-pmtiles-test overzoom-test flatgeobuf-test
+test: tippecanoe tippecanoe-decode $(addsuffix .check,$(TESTS)) raw-tiles-test parallel-test radix-sort-test pbf-test join-test enumerate-test decode-test join-filter-test unit json-tool-test allow-existing-test csv-test layer-json-test pmtiles-test decode-pmtiles-test overzoom-test flatgeobuf-test
 	./unit
 
 suffixes = json json.gz
@@ -168,7 +168,7 @@ nogeobuf = tests/overflow/out/-z0.json $(wildcard tests/stringid/out/*.json)
 geobuf-test: tippecanoe-json-tool $(addsuffix .checkbuf,$(filter-out $(nogeobuf),$(TESTS)))
 
 # For quicker address sanitizer build, hope that regular JSON parsing is tested enough by parallel and join tests
-fewer-tests: tippecanoe tippecanoe-decode geobuf-test raw-tiles-test parallel-test pbf-test join-test enumerate-test decode-test join-filter-test unit
+fewer-tests: tippecanoe tippecanoe-decode geobuf-test raw-tiles-test parallel-test radix-sort-test pbf-test join-test enumerate-test decode-test join-filter-test unit
 
 # XXX Use proper makefile rules instead of a for loop
 %.json.checkbuf:
@@ -179,6 +179,42 @@ fewer-tests: tippecanoe tippecanoe-decode geobuf-test raw-tiles-test parallel-te
 	cmp $@.out $(patsubst %.checkbuf,%,$@)
 	rm $@.out $@.mbtiles
 
+# The result of the sort must not depend on how the sort was performed, so instead of
+# checking the sorted output against an expected copy of it, check that sorting by radix
+# produces the same tiles as sorting in memory. --prefer-radix-sort lowers the memory
+# limit to 8K, which radix() then halves again, so the radix subdivision has to recurse
+# until each bucket is under 4K. How deeply that recurses depends on how many files the
+# machine will let us open at once, but the sorted result is the same either way, so
+# this comparison doesn't depend on that.
+#
+# What sends the sort down its rarely-taken paths is the shape of the input rather than
+# the size of it: a feature whose geometry alone is bigger than the memory limit is
+# sorted as a bucket of its own, and features that share a long run of leading index
+# bits have to be subdivided until there are no bits left. The first two inputs are
+# each one of those on purpose -- several separated features too big to sort in memory,
+# and many features at one location -- and the rest are for breadth.
+radix-sort-test: tippecanoe tippecanoe-decode
+	mkdir -p tests/radix-sort
+	perl -e 'for ($$f = 0; $$f < 8; $$f++) { print "{ \"type\": \"Feature\", \"properties\": { \"f\": $$f }, \"geometry\": { \"type\": \"LineString\", \"coordinates\": ["; for ($$i = 0; $$i < 2000; $$i++) { print "," unless $$i == 0; printf "[%f,%f]", $$f * 40 - 175 + $$i * 0.001, $$i % 2 * 0.5 - 20; } print "] } }\n"; }' > tests/radix-sort/bigfeatures.json
+	perl -e 'for ($$i = 0; $$i < 500; $$i++) { print "{ \"type\": \"Feature\", \"properties\": { \"i\": $$i }, \"geometry\": { \"type\": \"Point\", \"coordinates\": [ 17, 42 ] } }\n"; }' > tests/radix-sort/onelocation.json
+	$(MAKE) radix-sort-compare RADIXIN="tests/radix-sort/bigfeatures.json" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/radix-sort/onelocation.json" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/feature-filter/in.json" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/ne_110m_ocean/in.json" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/border/in.json" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/loop/in.json" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/tl_2022_11_tract/in.json.gz" RADIXARGS="-z4"
+	$(MAKE) radix-sort-compare RADIXIN="tests/epsg-3857/in.json" RADIXARGS="-z4 -sEPSG:3857"
+	rm -r tests/radix-sort
+
+radix-sort-compare:
+	./tippecanoe -q -f $(RADIXARGS) -o tests/radix-sort/memory.mbtiles $(RADIXIN)
+	./tippecanoe -q -f $(RADIXARGS) -aR -o tests/radix-sort/radix.mbtiles $(RADIXIN)
+	./tippecanoe-decode -x generator -x generator_options -x name -x description tests/radix-sort/memory.mbtiles > tests/radix-sort/memory.json
+	./tippecanoe-decode -x generator -x generator_options -x name -x description tests/radix-sort/radix.mbtiles > tests/radix-sort/radix.json
+	cmp tests/radix-sort/memory.json tests/radix-sort/radix.json
+	rm tests/radix-sort/memory.mbtiles tests/radix-sort/radix.mbtiles tests/radix-sort/memory.json tests/radix-sort/radix.json
+
 parallel-test: $(eval SHELL:=$(ADVSHELL))
 	mkdir -p tests/parallel
 	perl -e 'for ($$i = 0; $$i < 20; $$i++) { $$lon = rand(360) - 180; $$lat = rand(180) - 90; $$k = rand(1); $$v = rand(1); print "{ \"type\": \"Feature\", \"properties\": { \"yes\": \"no\", \"who\": 1, \"$$k\": \"$$v\" }, \"geometry\": { \"type\": \"Point\", \"coordinates\": [ $$lon, $$lat ] } }\n"; }' > tests/parallel/in1.json
diff --git a/main.cpp b/main.cpp
index b7813995..5842f988 100644
--- a/main.cpp
+++ b/main.cpp
@@ -743,8 +743,20 @@ void start_parsing(int fd, STREAM *fp, long long offset, long long len, std::ato
 }
 
 void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int splits, long long mem, const char *tmpdir, long long *availfiles, FILE *geomfile, FILE *indexfile, std::atomic<long long> *geompos_out, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, int basezoom, double droprate, double gamma, struct drop_state *ds) {
-	// Arranged as bits to facilitate subdividing again if a subdivided file is still huge
-	int splitbits = log(splits) / log(2);
+	// Arranged as bits to facilitate subdividing again if a subdivided file is still huge.
+	//
+	// There must be at least two buckets. With only one, each subdivision would
+	// consume no bits of the index, so it would never reach the maximum prefix
+	// that stops the recursion, and the shift that chooses a feature's bucket
+	// below would be by the full width of the index. That shift is undefined,
+	// and what it does in practice is to mask the shift count down to zero, so
+	// the bucket number comes out as the whole shifted index instead of as 0
+	// and the writes are made through whatever is found beyond the end of the
+	// arrays of buckets.
+	int splitbits = 1;
+	if (splits > 1) {
+		splitbits = log(splits) / log(2);
+	}
 	splits = 1 << splitbits;
 
 	FILE *geomfiles[splits];
@@ -891,7 +903,14 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 		}
 
 		if (indexst.st_size > 0) {
-			if (indexst.st_size + geomst.st_size < mem) {
+			// Subdividing again would only make progress if the next level could
+			// split into at least two buckets, which it can't if there are no longer
+			// enough files left to do it with. In that case, sort in memory instead,
+			// even though this is more memory than we wanted to use at once, since
+			// it is the only way left to get this bucket sorted.
+			bool can_subdivide = *availfiles / 4 > 1;
+
+			if (indexst.st_size + geomst.st_size < mem || !can_subdivide) {
 				std::atomic<long long> indexpos(indexst.st_size);
 				int bytes = sizeof(struct index);
 
@@ -994,7 +1013,11 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 					struct index ix = indexmap[a];
 					long long pos = *geompos_out;
 
-					fwrite_check(geommap + ix.start, ix.end - ix.start, 1, geomfile, geompos_out, "geom");
+					// MAGIC: This knows that the feature minzoom is the last byte of the serialized feature
+					// and is writing one byte less and then adding the byte for the minzoom,
+					// the same as merge() does.
+
+					fwrite_check(geommap + ix.start, 1, ix.end - ix.start - 1, geomfile, geompos_out, "geom");
 					int feature_minzoom = calc_feature_minzoom(&ix, ds, maxzoom, gamma);
 					serialize_byte(geomfile, feature_minzoom, geompos_out, "merge geometry");