Minecraft Java 26.3 · overworld TERRAIN step · RTX 5090

Minecraft terrain generation on CUDA, bit for bit

Biomes, density with every cave noise, aquifers and the material rules with ore veins, reimplemented as CUDA kernels. Every block, biome and density value matches vanilla exactly, checked offline against the game's own values and inside the running game.

chunk 16 × 384 × 16 blocks, y −64 to 319seed 123456789GPU RTX 5090CPU Ryzen 9 9950X, 32 threadsvanilla 26.3 on JDK 25

857×

faster than vanilla on 1 thread

49×

faster than vanilla on all 32 threads

101,352

chunks/s end to end on the GPU (9.9 µs per chunk)

mismatches in 1,782,828,972 values compared with vanilla

Results

Same work everywhere: biomes for the region plus a one-chunk border, then the noise fill with aquifers and the material rules for every chunk. Vanilla runs its own private methods on real ProtoChunks (see Method). End to end includes host preparation and copying every block back to pinned host memory.

configuration	chunks/s	ms/chunk	CV	min	max	n	notes
Vanilla 26.3, 1 thread (region 16²)	118	8.46	4.6%	106	124	15	3 JVM runs, run means 116–120
Vanilla 26.3, 32 threads (region 32²)	2,057	0.49	4.7%	1,910	2,235	15	3 JVM runs, run means 2,018–2,089
CPU port, 1 thread (region 16²)	66	15.16	4.0%	60	70	15	3 runs
CPU port, 32 threads (region 32²)	1,162	0.86	1.1%	1,138	1,191	15	3 runs
RTX 5090 end to end, batch 2048, 2 slots (region 128²)	101,352	0.0099	0.7%	99,857	102,644	15	3 runs, pinned output
RTX 5090 end to end, pageable output	93,730	0.0107	1.7%	92,913	96,499	5	1 run
RTX 5090 kernels only, batch 4096 (region 128²)	109,063	0.0092	1.1%	107,024	110,281	5	sum of kernel event times

Against our own C++ port on 32 threads the GPU is 87× faster. The port is the debugging reference, an interpreter of the density graph written for clarity and verified first, not a tuned CPU implementation: on one thread it runs at 0.56× vanilla's JIT-compiled speed. CV is the coefficient of variation over all timed repetitions; warm-up repetitions are excluded.

Vanilla reaches 74% of linear scaling at 16 threads (one per core); SMT adds another 47% at 32. The GPU line is the end-to-end rate.

Inside the GPU run

CUDA event times per stage for one batch of 4,096 chunks run on its own, without overlap.

Throughput needs a few thousand chunks in flight. With two slots the copy back overlaps the next batch's kernels (v9 in How it got fast).

For a handful of chunks the fixed cost of a call (47 kernel launches including two radix sorts) dominates: one chunk takes 2.63 ms, against 8.5 ms for vanilla on one thread. Median of 20 calls each.

The density core alone

The first milestone ported only the density and the block fill, without aquifers or material rules. Vanilla here is NoiseChunk + sampleVolume(finalDensity) + the substance decision per block; the GPU row uses the first, naive kernels, kept as the baseline implementation cuda-v0. On one thread vanilla spends 2.18 ms per chunk on this core against 8.46 ms for the whole TERRAIN step, so the core is only 26% of the work. Biomes, aquifers and the material rules are the larger part, which is why the full step was the target.

configuration	chunks/s	ms/chunk	CV	min	max	n
Vanilla 26.3 core, 1 thread	459	2.18	1.1%	455	467	5
Vanilla 26.3 core, 32 threads	8,586	0.12	3.8%	8,252	8,946	5
CPU port core, 1 thread	209	4.80	7.4%	192	227	5
CPU port core, 32 threads	3,967	0.25	4.7%	3,664	4,142	5
RTX 5090 core (naive v0 kernels), end to end	40,126	0.0249	0.5%	39,830	40,309	5

Watch it generate

Each implementation generated the same 32×32-chunk region (chunk 462, −296: badlands, a river, beach, ocean and grassland) while logging the moment every chunk was finished. The videos replay those logs at real speed over a top-down map of the generated blocks: highest block per column, water shaded by depth, relief from the height field. The map is the same for every implementation because the blocks are; for this region the GPU's output matches a vanilla dump block for block. These are single recorded runs; the tables in Results are the reference measurements.

All five on one clock. The GPU is done within the first frame (18.5 ms); vanilla on one thread needs 9.1 s.

Vanilla, 1 thread: 9.1 s, 112 chunks/s. Biomes for the region and its border come first, then chunks finish one by one in order.

Vanilla, 32 threads: 610 ms, 1,679 chunks/s. Chunks land out of order as worker threads pick them up.

C++ port, 1 thread: 13.1 s, 78 chunks/s.

C++ port, 32 threads: 856 ms, 1,197 chunks/s.

RTX 5090: 18.5 ms in four batches of 256 chunks, then a replay at 1/160 speed to show each batch land.

RTX 5090 at its benchmark batch size: a 128×128 region (16,384 chunks) in 170 ms, 96,246 chunks/s, eight batches of 2,048, then a replay at 1/25 speed.

In first person

Two recordings fly one 38-second camera path over the same seed: a walk at eye height in the badlands, then a climb over mesas, a forest, a river and the coast. The first is our own renderer generating the terrain live on the GPU; the second is the real game, flying through chunks the GPU generated.

Our CUDA renderer (viz/fly.cu) generates every chunk within 30 chunks of the camera with the bit-exact pipeline and ray traces the view on the same GPU: 2,821 chunks in the first 6 frames (74 ms), then 46 per generating frame in 4.4 ms on average (5.7 ms at most), 10,862 in all. Procedural colours, no game textures.

The Minecraft 26.3 client. A Fabric mod fills every chunk of a 136×136-chunk region with the GPU's blocks and biomes; the game adds its own caves, trees and lighting. 10,683 chunks came from the GPU, and at every frame the client held 0 chunks from anywhere else. 30 fps, captured frame by frame. How the mod works and what it checked: Correctness.

The same frames side by side: the same mesas, terracotta bands and rivers. The ray tracer shows the TERRAIN step alone; the client adds the game's later steps on top of it.

Correctness

Bit-exactness is the definition of done: every block state, every biome and every density value is compared for equality with vanilla, with no tolerance. Java is strict IEEE with no fused multiply-add, so device code is compiled with -fmad=false -prec-div=true -prec-sqrt=true -ftz=false, never fast-math. Every optimization keeps the exact IEEE operations or replaces them by integer work with a written proof that the result is identical.

CMakeLists.txtlines 24–27

24add_compile_options($<$<COMPILE_LANGUAGE:CUDA>:-fmad=false> $<$<COMPILE_LANGUAGE:CUDA>:-ftz=false>
25                    $<$<COMPILE_LANGUAGE:CUDA>:-prec-div=true> $<$<COMPILE_LANGUAGE:CUDA>:-prec-sqrt=true>
26                    $<$<COMPILE_LANGUAGE:CUDA>:-Xcompiler=/fp:precise> $<$<COMPILE_LANGUAGE:CUDA>:-Xcompiler=/utf-8> $<$<COMPILE_LANGUAGE:CUDA>:-Xcompiler=/Zc:preprocessor>
27                    $<$<COMPILE_LANGUAGE:CUDA>:--expt-relaxed-constexpr>)

Device flags for every CUDA file: no fused multiply-add, IEEE division and square root, denormals kept. The host side uses /fp:precise, and a unit test checks that the host compiler does not contract.

The evidence comes in three layers, from single values to a running game.

Against values from the game

The reference oracle runs the real game classes (Minecraft 26.3 server jar via Fabric Loom, on this machine only) and dumps values for 32 seeds (8 fixed, 24 random) at 18 chunks each: a 3×3 block at the origin, chunks near and far, chunks at the world border (±1,874,999) and at the 32-bit wrap points of the noise lattice, over the full build height (y = −64 to 319). mccuda_tests --set full (RESULT: PASS (bit-exact)):

level / check	values		notes
unit	72	PASS	Java arithmetic semantics, gradient table, MD5, host FP contraction
seedFromHashOf (MD5)	12	PASS
Xoroshiro128++ source	1,112	PASS
legacy LCG source	968	PASS
noise tables + samples	506,240	PASS
zoom fast path	2,147,483,648	PASS	GPU self-test against the plain FP64 BiomeManager zoom (not a vanilla comparison); 50492 near ties took the FP64 path
L1 named functions	29,491,200	PASS
L2 final density (cpu)	56,623,104	PASS
L3 block fill (cpu)	56,623,104	PASS
L2 final density (gen-host)	56,623,104	PASS
L3 block fill (gen-host)	56,623,104	PASS
L2 final density (cuda-v0)	56,623,104	PASS
L3 block fill (cuda-v0)	56,623,104	PASS
biomes (cpu)	884,736	PASS
aquifer fill (cpu)	56,623,104	PASS
terrain (cpu)	56,623,104	PASS
biomes (cuda-terrain)	884,736	PASS
aquifer fill (cuda-terrain)	56,623,104	PASS
terrain (cuda-terrain)	56,623,104	PASS

Whole regions, block for block

The benchmark harness dumps three contiguous regions from vanilla and compares them block for block with both generators:

generator	region	blocks
GPU	origin	25,165,824
CPU port	origin	25,165,824
GPU	negative	25,165,824
CPU port	negative	25,165,824
GPU	border	6,291,456
CPU port	border	6,291,456

Inside the running game

The in-game recording (Videos) uses a Fabric mod in a development client. Inside a 136×136-chunk region it puts the GPU's output in at vanilla's own seams: the noise fill writes the GPU's finished blocks, the surface pass is skipped because its result is already there, and biomes are swapped in the same way at doCreateBiomes. The client reads a dump the GPU wrote beforehand (1.8 GB); it does not call CUDA itself.

client/src/main/java/mccuda/client/mixin/NoiseGeneratorMixin.javalines 27–37

27@Inject(method = "doFill", at = @At("HEAD"), cancellable = true)
28private void mccuda$fill(NoiseChunk noiseChunk, ChunkAccess chunk, CallbackInfo ci) {
29    GpuDump dump = GpuDump.get();
30    if (dump == null) return;
31    if (dump.serves(chunk.getPos())) {
32        dump.fillBlocks(chunk);
33        ci.cancel();
34    } else {
35        dump.vanillaChunks.incrementAndGet();
36    }
37}

Inside the region the chunk is filled from the dump and vanilla's fill is cancelled; outside it the game generates as usual.

For control takes the same mod leaves generation to the game and compares each chunk with the GPU's dump as the server generates it: the biomes right after doCreateBiomes, and every block right after buildSurface, the end of the TERRAIN step.

client/src/main/java/mccuda/client/GpuDump.javalines 212–234

212public void verifyBlocks(ChunkAccess chunk) {
213    long base = (long) index(chunk.getPos()) * chunkBytes, bad = 0;
214    LevelChunkSection[] sections = chunk.getSections();
215    for (int y = 0; y < height; y++) {
216        LevelChunkSection sec = sections[y >> 4];
217        int ly = y & 15;
218        for (int z = 0; z < 16; z++)
219            for (int x = 0; x < 16; x++)
220                if (sec.getBlockState(x, ly, z) != palette[terrain.get((int) (base + ((long) y * 16 + z) * 16 + x)) & 0xFF]) bad++;
221    }
222    boolean adapted = chunk.getAllReferences().entrySet().stream().anyMatch(e -> e.getKey().terrainAdaptation() != TerrainAdjustment.NONE && !e.getValue().isEmpty());
223    checkedBlocks.addAndGet(256L * height);
224    blockMismatches.addAndGet(bad);
225    checkedBlockChunks.incrementAndGet();
226    if (adapted) {
227        adaptedChunks.incrementAndGet();
228        adaptedMismatches.addAndGet(bad);
229        if (bad > 0) adaptedBadChunks.incrementAndGet();
230    } else if (bad > 0) {
231        plainBadChunks.incrementAndGet();
232        System.out.printf("[mccuda] mismatch without structure adaptation: chunk %s, %d blocks%n", chunk.getPos(), bad);
233    }
234}

The game's own chunk after its surface pass, compared state by state with the GPU's palette index at the same position. Block states are canonical objects, so identity is equality.

take	structures	chunks	blocks compared	mismatches	biome values
vanilla	off	10,727	1,054,507,008	0	17,075,712
vanilla, second run	off	10,646	1,046,544,384	0	16,952,832
vanilla, structures on	on	10,663	1,048,215,552	826,034	16,977,408

With structures on, every mismatch is in a chunk whose density a structure reshapes (the Beardifier, for villages, trail ruins and the like): 524 of the 1,142 such chunks differ, 826,034 blocks in all, and the other 9,521 chunks match exactly. That adaptation is outside this project's scope, so the recordings use worlds without structures. One more thing the dump does not carry: vanilla's doFill also marks aquifer fluid blocks that should start flowing once the chunk is finished. In GPU-filled chunks those fluids stay still until a neighbouring block changes; the blocks themselves are identical.

The recorded frames still differ, and are expected to. The GPU take differs from the first vanilla take in 5.1% of pixels per frame on average (by more than 8 of 255), and the two vanilla takes differ from each other by about as much, 5.6%. The differences are trees. A chunk's trees are placed while its neighbours may already be decorated, so where trees compete for space the outcome depends on the order in which the server's worker threads finish chunks.

The same forest in three takes: GPU chunks, vanilla, and vanilla again — Frame 785 of three takes, cropped where the GPU take and the first vanilla take differ most. The trees change between the two vanilla runs as much as between GPU and vanilla; the ground under them does not.

What exactness takes

These are the behaviours of 26.3's terrain generation that are easy to get almost right. They were found by reading public signatures and bytecode of the locally prepared game jar and by comparing against the oracle; no game source is reproduced. The longer write-up is docs/ENGINE_NOTES.md.

Two evaluation semantics

26.3's density functions are float-valued and have two evaluation paths that do not produce identical bits. The volume path (sampleVolume) generates a chunk's density, the climate for biomes, ore-vein density and richness, and the aquifer's surface-level prefill. The point path (sampleValue) answers the aquifer's barrier, floodedness, spread and lava queries, the ore-vein filler gap and the find_top_surface probes.

	volume path	point path
noise coordinates	`block * (scale * frequency)`	`(block * scale) * frequency`
gradient dot product	`(gxx + gzz) + gy*y`, xz part reused over y	`(gxx + gyy) + gz*z`
interpolation	z, then x, then incremental steps in y	lerp3 at the point
min / max	keep the first argument unless the second is strictly smaller or larger	`Math.min` / `max`, with range-based short cuts
mul / div	plain	zero short-circuit on the first operand

The code generator (cuda/codegen/gen_router.py) emits a separate evaluator per context, each with fixed semantics. Mixing them up changes the last bit at a small fraction of positions, which is enough to flip blocks. The gradient dot product shows how small the difference is:

cpu/include/mccuda/noise_eval.hlines 74–87

74// GRADIENT[hash & 15].dot(x, y, z) with float arguments: gx * x + gy * y + gz * z
75MCC_HD MCC_FI float grad_dot3(int hash, float x, float y, float z) {
76    float a = static_cast<float>(grad_c(hash, 0)) * x;
77    float b = static_cast<float>(grad_c(hash, 1)) * y;
78    float c = static_cast<float>(grad_c(hash, 2)) * z;
79    float s = a + b;
80    return s + c;
81}
82
83// Gradient.dotXz(x, z) = this.x * x + this.z * z
84MCC_HD MCC_FI float grad_dot_xz(int hash, float x, float z) {
85    float a = static_cast<float>(grad_c(hash, 0)) * x;
86    float c = static_cast<float>(grad_c(hash, 2)) * z;
87    return a + c;

Point samples add y before z; volume samples compute the xz part once and reuse it over y. In float arithmetic the two orders round differently.

The compiled density graph

Vanilla rewrites the density graph before sampling: references are inlined, structurally equal cache nodes are merged, and functions that do not depend on an axis are evaluated on slices. The IR exporter (oracle/.../IrExporter.java) mirrors these rewrites, and the binary operations carry vanilla's compile-time decisions: constant detection, division by a constant as a multiply by the float reciprocal, and min/max pruning from value ranges. Those decisions change results, so they are made once, at export time, exactly as vanilla makes them.

Random numbers

Xoroshiro128++ with positional factories: fromHashOf seeds from the MD5 of a name, at(x, y, z) from Mth.getSeed. The legacy 48-bit LCG still drives the biome temperature noise. nextDouble is (bits >>> 11) * 2^-53 computed in double: a decompiler rendering suggested float arithmetic, and the bytecode (l2d, dmul) settled it. Noise construction sums the octave amplitudes with a compensated sum, so plain summation is wrong.

Biomes and the R-tree tie rule

The six climate parameters are sampled at every quart (4 blocks) with the volume path and quantized as (long)(value * 10000f). The lookup is an R-tree over 7,594 parameter boxes (19 children per node, built with stable sorts): the search descends into a child only if its box is strictly closer than the best leaf so far, and accepts a leaf only if it is strictly closer. BiomeManager then picks a quart per block from the 8 surrounding quarts with a seeded offset; that seed is the SHA-256 of the world seed, read little-endian.

The tree keeps a last-result candidate and starts each search from it. On an exact distance tie the candidate wins if it reaches the minimal distance; otherwise the first minimal leaf in DFS order wins. This is not a corner case: over 5 seeds (28,372,992 quart lookups each, region 128×128 plus borders) the candidate changed the chosen leaf at 547–703 ppm of lookups and the biome at 329–461 ppm, because quantized climate values regularly land exactly on a box boundary. Returning the plain nearest leaf would put a wrong biome in roughly one quart in 2,600. Both generators reproduce the rule along each chunk's chain of lookups. Vanilla's only other history dependence is the first lookup of a chunk, which starts from whatever the worker thread looked up last; the oracle and our generators start it empty, so a live server can differ where that first lookup is an exact tie. Inside the running game none of the 17,075,712 biome values compared above differed.

The GPU finds the first minimal leaf of every quart in parallel (with a bounded search, v8 below). The rule then needs the previous result only where the previous leaf also reaches the minimum:

cuda/src/terrain.culines 358–367

358// tie[p]: the previous position's L* is also at minimal distance here (only then can the
359// R-tree's last-result candidate make the result differ from L*).
360__global__ void k_biome_ties(int nb, const int64_t* targets, const int32_t* first_min, const int64_t* min_dist, uint8_t* tie) {
361    const int per = 16 * (c_tt.height / 4);
362    int idx = blockIdx.x * blockDim.x + threadIdx.x;
363    if (idx >= nb * per) return;
364    int b = idx / per, p = idx % per;
365    if (p == 0) { tie[idx] = 0; return; }
366    size_t i = static_cast<size_t>(b) * per + chain_quart(p), ip = static_cast<size_t>(b) * per + chain_quart(p - 1);
367    tie[idx] = tree_distance(c_tt.leaf_node[first_min[ip]], targets + i * 7) == min_dist[i] ? 1 : 0;

Only where this flag is set can the last-result candidate change the answer, and only there does the chain kernel take its slow path; one warp per chunk resolves those windows.

Aquifers

A 16×12×16 grid. Cell locations are random per grid cell; fluid status comes from the surface level around the cell, the floodedness and spread noises, and a lava-noise check for fluid levels at or below −10. The surface level comes from find_top_surface: positions inside the aquifer constructor's prefill window (quarts 4c−4 to 4c+6) use volume semantics, other positions are sampled lazily with point semantics, and the GPU must use the same split. Blocks above a per-chunk skip height use the global fluid picker directly.

Material rules

A column is walked top-down. Stone depth above and below, water height and surface depth feed the rule tree: about 150 conditions, 110 block results, 43 noise thresholds, 2 vertical gradients, 2 ore veins (new in 26.3's material rules) and clay bands. 26.3 checks a cave biome for almost every stone block below the preliminary surface, so the biome zoom runs for most blocks of a chunk.

The heightmap is updated on every block placement, and the slope condition reads neighbour heights. Vanilla processes columns in x-major, z-minor order, so a neighbour's height is from before or after both surface extensions. The GPU walks all columns in parallel and reconstructs what each column would have seen from three height snapshots, with the frozen-ocean extension dry-run in a first pass:

cuda/src/terrain.culines 905–914

905// Height of another column as seen when vanilla processes (x, z): columns earlier in
906// the x-major, z-minor loop have finished both extensions; later ones none.
907auto H = [&](int ox, int oz) {
908    int j = ci * 256 + ox * 16 + oz;
909    int v;
910    if (ox == x && oz == z) v = h1[j];
911    else if (ox < x || (ox == x && oz < z)) v = h2[j];
912    else v = h0[j];
913    return v - 1;                                        // ChunkAccess.getHeight
914};

Columns before (x, z) in the x-major, z-minor loop have finished both surface extensions (h2), later ones neither (h0), and the column itself has finished the eroded-badlands extension (h1). The slope condition reads these heights.

How it got fast

Starting from a straight port of the reference code (v0), each step was chosen from profiling data, measured right after the change on the same 64×64 region (4,096 chunks), and re-verified on the full oracle set before it was committed. The headline rate uses the 128×128 region, where the two-slot pipeline runs eight batches instead of two.

v9 is measured with two pipeline slots and batch 2048; its stage bars come from a one-slot run of the same code.

version	change	chunks/s	step	vs v0	why the output stays identical
v0	naive port	20,450		1.0×	reference code, ported directly
v1	FP64-free floor/wrap	23,748	+16%	1.2×	the same results from integer work on the IEEE bits
v2	axis tables	40,774	+72%	2.0×	the same IEEE operations, computed once per distinct coordinate
v3	zoom table	52,958	+30%	2.6×	the same expression, tabulated over its 8 × 1,024 inputs
v4	grid dedupe	60,369	+14%	3.0×	every value is a pure function of position and semantics
v5	column dedupe	64,398	+7%	3.1×	every value is a pure function of position
v6	zoom FP32 filter	65,019	+1%	3.2×	proven error bound; near ties take the FP64 comparison
v7	inline rules	70,995	+9%	3.5×	code placement only
v8	biome search	76,489	+8%	3.7×	any bound above the minimum yields the same first minimal leaf
v9	pipelined, 2 slots	92,559	+21%	4.5×	scheduling only

What mattered, with the code

FP64 is the bottleneck on a GeForce part

The RTX 5090 runs FP64 at 1/64 of its FP32 rate. Minecraft 26.3's gradient noise is FP32 inside, but each octave sets up its three coordinates in FP64, about 29 FP64 operations per octave. A static SASS count of the v0 corner kernel showed about 7,000 FP64 arithmetic, compare and conversion instructions against about 16,000 FP32 ones, so at 1:64 the FP64 work dominates by an order of magnitude.

cpu/include/mccuda/noise_eval.hlines 162–178

162// PerlinNoise.addToVolume: the value added for one sample (before the amplitude multiply).
163MCC_HD MCC_FI float perlin_volume(const PerlinOctave& o, int32_t bx, int32_t by, int32_t bz, double xzScale, double yScale) {
164    double z = gn_wrap(static_cast<double>(bz) * xzScale) + o.oz;
165    int32_t fz = mth_floor(z);
166    float rz = static_cast<float>(z - static_cast<double>(fz));
167    float az = mth_smoothstep(rz);
168    double x = gn_wrap(static_cast<double>(bx) * xzScale) + o.ox;
169    int32_t fx = mth_floor(x);
170    float rx = static_cast<float>(x - static_cast<double>(fx));
171    float ax = mth_smoothstep(rx);
172    double y = gn_wrap(static_cast<double>(by) * yScale) + o.oy;
173    int32_t fy = mth_floor(y);
174    float ry = static_cast<float>(y - static_cast<double>(fy));
175    float ay = mth_smoothstep(ry);
176    VolumeCorners c = volume_corners(o, fx, fy, fz, rx, rz);
177    return volume_lerp(c, ax, ay, az, ry);
178}

Every double line is FP64 work: a multiply, the wrap test, an add, a floor and a subtraction per axis. Everything after volume_corners is FP32.

v1: floor and wrap as integer work

Mth.floor is (int) Math.floor(v), and Java's cast saturates with NaN going to 0. On the device it becomes one conversion instruction with the NaN case decided on the bit pattern, and the wrap range test becomes integer compares on the IEEE bits:

cpu/include/mccuda/javamath.hlines 68–79

68// Mth.floor(double) = (int) Math.floor(v)
69MCC_HD MCC_FI int32_t mth_floor(double v) {
70#if defined(__CUDA_ARCH__)
71    // One cvt.rmi.s32.f64: floors and saturates exactly like Java's (int) cast.
72    // NaN -> 0 is decided on the bit pattern, which keeps the FP64 pipe (1/64 rate
73    // on GeForce parts) out of it.
74    if ((d2bits(v) & 0x7FFFFFFFFFFFFFFFull) > 0x7FF0000000000000ull) return 0;
75    return __double2int_rd(v);
76#else
77    return java_d2i(std::floor(v));
78#endif
79}

cvt.rmi.s32.f64 floors and saturates exactly like Java's (int) cast. The host keeps the reference form.

cpu/include/mccuda/noise_eval.hlines 56–64

56MCC_HD MCC_FI double gn_wrap(double x) {
57#if defined(__CUDA_ARCH__)
58    // -H <= x < H decided on the bit pattern (IEEE magnitudes order like integers):
59    // integer ops instead of two FP64 compares. -0.0 and -H pass, +H and NaN fail.
60    const uint64_t b = d2bits(x), mag = b & 0x7FFFFFFFFFFFFFFFull, hb = 0x416FFFFFFFFFFFFFull;   // bits of kHalfRoundOff
61    if ((b >> 63) ? mag <= hb : mag < hb) return x;
62#else
63    if (x >= -kHalfRoundOff && x < kHalfRoundOff) return x;
64#endif

IEEE magnitudes order like unsigned integers, so -H <= x < H needs no FP64 compare; the sign bit picks the inclusive or exclusive bound.

v2: per-axis coordinate tables

Every corner-level noise is unshifted and at a fixed scale, so each axis's coordinate setup depends on that axis alone. Its exact result is computed once per distinct coordinate and read back:

cpu/include/mccuda/noise_eval.hlines 223–238

223// In addToVolume each axis's coordinate setup (scale, wrap, offset, floor, fraction,
224// smoothstep) depends on that axis's block coordinate only, and it is the only FP64
225// work in a sample. AxisRec holds its exact result, so it can be computed once per
226// distinct coordinate and reused (identical arithmetic, so identical bits).
227struct AxisRec {
228    int32_t f;   // lattice cell (floor)
229    float r;     // relative coordinate handed to the gradients / lerp (fudged for smeared y)
230    float a;     // smoothstep weight
231};
232
233// x or z of PerlinNoise/SmearedPerlinNoise.addToVolume, and y of PerlinNoise.addToVolume.
234MCC_HD MCC_FI AxisRec axis_plain(int32_t b, double scale, double offset) {
235    double v = gn_wrap(static_cast<double>(b) * scale) + offset;
236    int32_t f = mth_floor(v);
237    float r = static_cast<float>(v - static_cast<double>(f));
238    return {f, r, mth_smoothstep(r)};

The record holds what one axis contributes to a sample. axis_plain is the same arithmetic as the axis setup in perlin_volume above, so the bits are identical.

cpu/include/mccuda/noise_eval.hlines 252–264

252// stack_add_volume with each layer's coordinate setup read from per-axis records.
253MCC_HD MCC_FI float stack_volume_tab(float acc, const StackLayer* layers, int count, const PerlinOctave* octaves,
254                                     const AxisRec* ax, const AxisRec* ay, const AxisRec* az, float amplitude) {
255    for (int i = 0; i < count; ++i) {
256        const StackLayer& L = layers[i];
257        float amp = amplitude * L.amplitude;
258        const AxisRec X = ax[i], Y = ay[i], Z = az[i];
259        VolumeCorners c = volume_corners(octaves[L.octave], X.f, Y.f, Z.f, X.r, Z.r);
260        float v = volume_lerp(c, X.a, Y.a, Z.a, Y.r);
261        acc = acc + amp * v;
262    }
263    return acc;
264}

The sample itself is now FP32 only: lattice hashes, gradients and lerps, with the coordinate setup read from three records.

cuda/src/terrain.culines 496–513

496__global__ void k_corners(int upper, const uint64_t* ucols, const int* count, const float* col, float* corners, int rows, const AxisRec* axis_x,
497                          const AxisRec* axis_y, const AxisRec* axis_z, int sl) {
498    int idx = blockIdx.x * blockDim.x + threadIdx.x;
499    if (idx >= upper * rows) return;
500    int u = idx / rows, row = idx % rows;
501    if (u >= *count) return;
502    int x, z;
503    col_xz(ucols[u], x, z);
504    int y = c_tt.min_y + row * 8;
505    const float* cv = col + static_cast<size_t>(u) * gen::kNumColumns;
506    const AxisRec* ax = axis_x + static_cast<size_t>(u) * sl;
507    const AxisRec* az = axis_z + static_cast<size_t>(u) * sl;
508    const AxisRec* ay = axis_y + static_cast<size_t>(row) * sl;
509    float out[kAllChannels];
510    gen::eval_corners_tab(c_tt.rt, x, y, z, cv, ax, ay, az, out);
511    if (gen::kNumOreChannels > 0) gen::eval_ore_corners_tab(c_tt.rt, x, y, z, cv, ax, ay, az, out + gen::kNumChannels);
512    for (int ch = 0; ch < kAllChannels; ++ch) corners[(static_cast<size_t>(u) * kAllChannels + ch) * rows + row] = out[ch];
513}

Each distinct corner column u owns its x and z records; the y records are shared by every column at the same corner row.

v3 and v6: the biome zoom

The material rules check a cave biome for nearly every stone block, and the BiomeManager zoom was 79% of the surface kernel in an ablation build. Its squared offsets take only 8 × 1,024 distinct values (v3), and an FP32 pre-filter with a proven error bound leaves the FP64 comparison to near ties (v6):

cuda/include/mccuda_cuda/terrain_device.hlines 107–116

107// Squared fiddled offset of one axis, (d + fiddle)^2, for the 8 possible offsets d
108// (fraction (a & 3) / 4, minus 1 on the far side) and the 1024 fiddle indices:
109// table index ((a & 3) + 4 * far) * 1024 + m. Same expression as fiddled_distance.
110MCC_HD MCC_FI double zoom_square(int t) {
111    int k = t >> 10, m = t & 1023;
112    double f = static_cast<double>(k & 3) / 4.0;
113    double d = k < 4 ? f : f - 1.0;
114    double v = fiddle_value(m);
115    return (d + v) * (d + v);
116}

The table is built on the device by this function, with the same expression as the per-lookup code it replaces.

cuda/include/mccuda_cuda/terrain_device.hlines 164–207

164// zoom_quart_tab with an exact FP32 pre-filter (sqf = the squares rounded to float).
165//
166// The 8 distances are first computed in FP32. Every squared term lies in [0, 2.1025]
167// ((1 + 0.45)^2), so a distance S is at most 6.31. Three table roundings and two FP32
168// additions give |f - S| <= 3u * S <= 1.13e-6 (u = 2^-24), and vanilla's FP64 value d is
169// within 1.4e-15 of S, so |f - d| < 1.14e-6. A corner j with f_j > lim = fl(fmin + 8e-6)
170// (that rounding is <= 3.8e-7) therefore has
171//   d_j > fmin + 8e-6 - 3.8e-7 - 1.14e-6 > fmin + 1.14e-6 + 5e-6 > d_k + 5e-6
172// for k = argmin f: it can neither be the minimum nor tie it. If one corner remains it
173// is vanilla's choice; otherwise the remaining corners are compared with vanilla's FP64
174// arithmetic, in corner order, exactly as zoom_quart_tab does.
175MCC_HD MCC_FI void zoom_quart_fast(int64_t seed, const double* sq, const float* sqf, int x, int y, int z, int& bx, int& by, int& bz) {
176    int ax = x - 2, ay = y - 2, az = z - 2;
177    int px = ax >> 2, py = ay >> 2, pz = az >> 2;
178    int tx = (ax & 3) << 10, ty = (ay & 3) << 10, tz = (az & 3) << 10;
179    int ix[8], iy[8], iz[8];
180    float f[8];
181    float fmin = bits2f(0x7F800000u);   // +infinity
182#pragma unroll
183    for (int i = 0; i < 8; ++i) {
184        bool xe = (i & 4) == 0, ye = (i & 2) == 0, ze = (i & 1) == 0;
185        int m = zoom_corner_m(seed, xe ? px : px + 1, ye ? py : py + 1, ze ? pz : pz + 1);
186        ix[i] = tx + (xe ? 0 : 4096) + (m & 1023);
187        iy[i] = ty + (ye ? 0 : 4096) + ((m >> 10) & 1023);
188        iz[i] = tz + (ze ? 0 : 4096) + (m >> 20);
189        f[i] = (sqf[iz[i]] + sqf[iy[i]]) + sqf[ix[i]];
190        fmin = f[i] < fmin ? f[i] : fmin;
191    }
192    const float lim = fmin + 8.0e-6f;
193    int candidates = 0, min_i = 0;
194#pragma unroll
195    for (int i = 7; i >= 0; --i)
196        if (f[i] <= lim) { ++candidates; min_i = i; }
197    if (candidates > 1) {
198        uint64_t min_bits = 0x7FF0000000000000ull;   // +infinity
199#pragma unroll
200        for (int i = 0; i < 8; ++i) {
201            if (f[i] <= lim) {
202                uint64_t d = d2bits(sq[iz[i]] + sq[iy[i]] + sq[ix[i]]);
203                if (min_bits > d) { min_i = i; min_bits = d; }
204            }
205        }
206    }
207    bx = (min_i & 4) == 0 ? px : px + 1;

The comment holds the error bound. One FP32 candidate is the answer; several fall back to vanilla's FP64 arithmetic in corner order.

bench/results/final/tests_full.txtlines 11–11

11zoom fast path                     PASS  (2147483648 checked, 0 mismatches; 50492 near ties took the FP64 path)

The self-test compares the fast path with the plain FP64 zoom on 2^31 random blocks up to the world border. The near ties exercise the FP64 fallback.

v4 and v5: deduplication across chunks

Neighbouring chunks need the same positions. For the 64×64 region, 2,117,632 surface-level lookups are 144,778 distinct ones (14.6×), and 102,400 corner columns are 66,049 (1.55×). The dedupe runs on the device for any chunk list:

cuda/src/dedupe.cuhlines 64–77

64// unique[0 .. *count) = distinct keys in increasing order; map[i] = index of keys[i] in unique.
65void run(const uint64_t* keys, int n, uint64_t* unique, int* map, int* count, cudaStream_t s = 0) {
66    if (n > cap_) throw std::runtime_error("dedupe capacity exceeded");
67    if (n == 0) { chk(cudaMemsetAsync(count, 0, sizeof(int), s), "memset"); return; }
68    const int tb = 256, g = (n + tb - 1) / tb;
69    dedupe_detail::k_iota<<<g, tb, 0, s>>>(n, idx_);
70    size_t t = temp_bytes_;
71    chk(cub::DeviceRadixSort::SortPairs(temp_, t, keys, sorted_, idx_, sorted_idx_, n, 0, bits_, s), "radix sort");
72    dedupe_detail::k_heads<<<g, tb, 0, s>>>(n, sorted_, flags_);
73    t = temp_bytes_;
74    chk(cub::DeviceScan::InclusiveSum(temp_, t, flags_, incl_, n, s), "scan");
75    dedupe_detail::k_scatter<<<g, tb, 0, s>>>(n, sorted_, sorted_idx_, flags_, incl_, unique, map, count);
76    chk(cudaGetLastError(), "dedupe kernels");
77}

Radix sort, mark the first key of each run, scan, scatter. The unique count stays on the device, so the evaluation kernels read it without a host round trip.

cuda/src/terrain.culines 526–540

526__global__ void k_grid_keys(const ChunkXZ* chunks, int n, uint64_t* keys) {
527    int idx = blockIdx.x * blockDim.x + threadIdx.x;
528    if (idx >= n * kGridKeys) return;
529    int c = idx / kGridKeys, e = idx % kGridKeys;
530    int cx = chunks[c].x, cz = chunks[c].z;
531    if (e < kGridCells) {
532        int qx = 4 * cx - 16 + e % kGridX, qz = 4 * cz - 8 + e / kGridX;
533        // Positions inside NoiseBasedAquifer's constructor prefill use the volume path.
534        bool prefill = qx >= 4 * cx - 4 && qx <= 4 * cx + 6 && qz >= 4 * cz - 4 && qz <= 4 * cz + 6;
535        keys[idx] = grid_key(prefill, qx, qz);
536    } else {
537        int k = e - kGridCells;   // chunk_surface_level corner: volume path at (16cx + 16(k & 1), 16cz + 16(k >> 1))
538        keys[idx] = grid_key(true, 4 * cx + 4 * (k & 1), 4 * cz + 4 * (k >> 1));
539    }
540}

The key carries the evaluation semantics: a position inside a chunk's aquifer prefill window is a volume-path sample, outside it a point-path sample (see Exactness).

v7: inline the generated rule tree

The material-rule functions were noinline, which kept the whole rule context in local memory. Inlining them cut the stage by 38%:

cuda/src/terrain.culines 32–36

32#if defined(__CUDACC__)
33#define MCC_RULE_FN __device__ __forceinline__
34#else
35#define MCC_RULE_FN inline
36#endif

The generator emits 469 condition and rule functions with this qualifier, one per node of the material-rule tree; the change is this one line.

v8: a bounded R-tree search and a warp-parallel tie chain

A search seeded with a neighbouring quart's leaf distance plus one returns the same first minimal leaf and prunes most of the tree:

cuda/src/terrain.culines 257–268

257// Climate.RTree.SubTree.search (iterative DFS with an explicit stack), returning the
258// first leaf in DFS order at minimal distance (L*), which is what vanilla returns for an
259// empty candidate. bound must exceed the minimal distance; then the result is L* for any
260// bound: any ancestor of L* is reached while min_d > m >= its box distance (min_d == m
261// would mean an earlier leaf at distance m), so it is entered, and L* is accepted. A tight
262// bound (the distance of a nearby quart's leaf, plus one) prunes most of the tree.
263__device__ int tree_search_first_min(const int64_t* t, int64_t bound, int64_t* out_dist) {
264    struct Frame { int node; int child; };
265    Frame stack[16];
266    int sp = 0;
267    int64_t min_d = bound;
268    int closest = -1;

The proof sits in the comment; the only change from an unbounded search is the initial min_d.

cuda/src/terrain.culines 388–392

388const unsigned any_tie = __ballot_sync(0xFFFFFFFFu, t);
389const int last = min(32, per - p0) - 1;
390if (any_tie == 0 && prev == prev_lstar) {
391    if (valid) biomes[i] = static_cast<uint8_t>(c_tt.leaf_biome[lstar]);
392    prev = prev_lstar = __shfl_sync(0xFFFFFFFFu, lstar, last);

The common case in the chain kernel: a window of 32 positions with no tie and in sync writes all 32 results at once.

v9: overlap the copy with the next batch

The first pipelined version overlapped nothing: copies run in submission order, so a download enqueued before the next batch's upload held that upload back. Downloads are now enqueued one batch later:

cuda/src/terrain.culines 1395–1414

1395Slot* prev = nullptr;
1396I.batch_done_ms.assign(static_cast<size_t>((count + I.max_batch - 1) / I.max_batch), 0.0);
1397check(cudaEventRecord(I.t0_event, I.slots[0]->stream), "event");
1398int batch = 0;
1399for (int off = 0; off < count; off += I.max_batch, ++batch) {
1400    int n = std::min(I.max_batch, count - off);
1401    Slot& S = *I.slots[I.next_slot];
1402    I.next_slot = (I.next_slot + 1) % I.slots.size();
1403    if (prev == &S) {   // single slot: its batch's download goes first
1404        I.download(S);
1405        prev = nullptr;
1406    }
1407    I.collect(S, timing);   // the slot's previous batch must be finished before its buffers are reused
1408    S.batch_index = batch;
1409    I.run_batch(S, chunks + off, n, terrain_out + static_cast<size_t>(off) * I.per, biomes_out ? biomes_out + static_cast<size_t>(off) * I.qper : nullptr,
1410                aquifer_out ? aquifer_out + static_cast<size_t>(off) * I.per : nullptr);
1411    if (prev) I.download(*prev);   // after this batch's uploads (see download())
1412    prev = &S;
1413}
1414if (prev) I.download(*prev);

Batch k's download is enqueued after batch k+1's uploads and kernels, so it overlaps that batch's kernels. Found from an Nsight Systems timeline.

Ablations

Throwaway builds that remove one piece of work, to find what a kernel spends its time on (results not bit-exact by design):

experiment	code	kernel	baseline ms	ablated ms	share
E1 ore veins disabled	v2	`k_surface2`	39.96	39.45	1%
E2 biome zoom bypassed	v2	`k_surface2`	39.96	8.23	79%
E3 R-tree search skipped	v7	`k_biome_quarts`	11.90	0.41	97%

Tried and reverted

Per-column memo of the 2D rule noises: no measurable change.
Per-column cache of zoom corner hashes, and a per-chunk hash table: both slower (local-memory traffic and scattered loads beat recomputing).
Coalesced column mapping in the surface kernels: neutral on its own (kept, it is the cleaner access pattern).

Final kernel profile

Nsight Systems, region 64×64 in one batch of 4,096 chunks, average over repetitions:

kernel	ms	share	launches	what it does
`k_blocks`	10.218	24.9%	1	interpolation, final density and aquifer substance per block
`k_corners`	9.903	24.1%	1	density and ore channels at the 4x8x4 cell corners
`k_biome_quarts`	8.043	19.6%	1	climate at every quart + bounded R-tree search (first minimal leaf)
`k_surface2`	7.990	19.4%	1	material rules (incl. ore veins) and frozen-ocean extension
`k_biome_chain`	1.553	3.8%	1	R-tree last-result rule (one warp per biome chunk)
`k_aquifer_cells`	0.998	2.4%	1	aquifer cell locations and fluid status
`k_fts_unique`	0.562	1.4%	1	find_top_surface per distinct surface-level cell
`k_heights`	0.418	1.0%	1	initial heightmap
14 more kernels	1.423	3.5%	39	y-independent inputs and axis records per distinct column, R-tree tie flags, the eroded-badlands pass, aquifer skip heights, and the dedupe (keys, CUB radix sort and scan, run heads, scatters)

Method and scope

Method

Machine: AMD Ryzen 9 9950X (16 cores / 32 threads), 128 GB RAM, NVIDIA GeForce RTX 5090 (170 SMs, 32 GB, driver 591.86, WDDM), Windows 11. CUDA 13.1, MSVC 14.44, JDK 25.0.4.1. GPU state sampled right after the GPU runs: 2910 MHz SM (max 3090 MHz), 44 °C, 265.28 W, P1.
Clocks unlocked (locking needs administrator rights); variance is reported instead: CV over repetitions and, for the headline configurations, the spread of independent process runs.
Workload: seed 123456789, square regions anchored at chunk (0, 0). Vanilla and the CPU port use 16×16 (1 to 4 threads) or 32×32 chunks; the GPU uses 128×128 (it needs the larger region to fill the device). Warm-up repetitions run first (JIT, allocations) and are not timed.
Fairness: vanilla computes biomes once per chunk of the region plus border, exactly what its BIOMES step does; the GPU recomputes the biomes of tile borders shared between batches, which is extra work on its side. Vanilla's worker threads share one RandomState like a server; its R-tree last-result cache is left alone (live-server behaviour) in throughput runs and reset per chunk in the dump runs.
Profiling: Nsight Systems for kernel times and timelines, cuobjdump for register/stack use and SASS instruction mixes, CUDA events per stage, and targeted ablation builds. Nsight Compute could not collect counters: it reports ERR_NVGPUCTRPERM (counter access is admin-only on this machine and was not changed).

Vanilla is timed through its own code path. The inner loop of the vanilla benchmark makes the calls that buildTerrain makes before carvers, on real ProtoChunks whose biomes were generated in phase 1:

oracle/src/main/java/mccuda/oracle/TerrainBench.javalines 142–159

142// phase 2: noise fill + surface for the region
143AtomicInteger next2 = new AtomicInteger();
144futures.clear();
145for (int t = 0; t < threads; t++) {
146    futures.add(pool.submit(() -> {
147        int i;
148        while ((i = next2.getAndIncrement()) < interior) {
149            int cx = x0 + i % region, cz = z0 + i / region;
150            ProtoChunk center = chunks.get(ChunkPos.pack(cx, cz));
151            BiomeManager biomeManager = new BiomeManager(
152                (qx, qy, qz) -> chunks.get(ChunkPos.pack(qx >> 2, qz >> 2)).getNoiseBiome(qx, qy, qz), zoomSeed);
153            Set<Holder<Biome>> possible = new ObjectOpenHashSet<>();
154            for (int dz = -1; dz <= 1; dz++)
155                for (int dx = -1; dx <= 1; dx++) chunks.get(ChunkPos.pack(cx + dx, cz + dz)).collectBiomesInPalette(possible);
156            DensityVolume vol = new DensityVolume(16, height, 16, cx * 16, minY, cz * 16);
157            try (NoiseChunk nc = new NoiseChunk(rs, Beardifier.EMPTY, stock, picker, Blender.empty(), vol)) {
158                doFill.invoke(generator, nc, center);
159                buildSurface.invoke(generator, center, nc, rs, biomeManager, possible, materialRule);

One NoiseChunk per chunk, doFill with aquifers on, then buildSurface with a BiomeManager over the 3×3 neighbourhood. The GPU does exactly this work; nothing is skipped on either side.

Scope

Included: biomes (multi-noise climate, exact R-tree search including its last-result tie rule), final density with all cave noises, NoiseBasedAquifer, doFill, the material rules with ore veins, clay bands, eroded-badlands and frozen-ocean extensions, and the heightmap updates they depend on.
Not included: carvers (they were not requested), the terrain adaptation around structures (measured in the game, see Correctness), blending with old-world chunks, features and lighting. A live server's first biome lookup per chunk can differ on an exact tie (see Exactness).
Decisions taken without the author: Minecraft 26.3 on native Windows; game-derived data (IR, terrain model, oracle dumps) is extracted locally into the git-ignored data/ and never committed; clocks unlocked; overworld only. A Fabric mod (client/) exists only to record the in-game videos from a precomputed dump; running CUDA inside the game (JNI) was not attempted.

Reproduce

powershell -File build.ps1               # CMake + Ninja, Release
powershell -File oracle\regen.ps1        # reference data from the local game jar (data/, git-ignored)
build\tests\mccuda_tests.exe --set full  # bit-exactness, all levels
python bench\run_all.py                  # every measurement in this report
python bench\make_videos.py record       # top-down videos: run every implementation on the video region
python bench\make_videos.py render       # top-down videos: render into docs/videos/
python bench\make_flythrough.py          # first-person GPU flythrough
cd client; .\gradlew runClient           # in-game take; control takes in bench/make_ingame.py
python bench\make_ingame.py              # in-game videos, frame comparison, in-game checks
python docs\make_report.py               # this page and docs/img/*.png

oracle/regen.ps1 writes, under the git-ignored data/, the optimized density graph (overworld_ir.json: floats as raw bits, doubles as hex, with value ranges and axis dependence per node), the material-rule tree, block palette and biome parameters (overworld_terrain.json), and the per-seed oracle dumps: named-function volumes (L1), final density (L2), core block fill (L3), biomes, aquifer fill and final terrain, in ((y − minY) * 16 + z) * 16 + x order with SHA-256 manifests. The in-game region comes from build\bench\mccuda_bench.exe --mode terrain-gpu --region 136 --x0 404 --z0 -324 --reps 1 --warmup 0 --dump data\videos\ingame_dump.

Generated 2026-10-09 by docs/make_report.py from the measurements in bench/results/. Source: github.com/bonquifo/minecraft-terrain-cuda. Not an official Minecraft product. Not approved by or associated with Mojang or Microsoft.