12 KiB
Fixed-seed CUDA FFT experiment (2026-10-07)
This is archived diagnostic code, not part of the active build. The independent implementation follows the bitmap/correlation approach in the public FFT source. It uses our exact Java rejection-loop arithmetic, overlapping tiles, full-result comparison with a direct GPU scan, and exact candidate validation. It handles solid squares only; arbitrary-pattern extensions were deliberately deferred.
RTX 3060 Laptop GPU, seed 61692151016558, origin (46973,415028), 5x5, one warm-up plus five measured repeats; table reports median device times. Setup is separately timed and excluded from these device timings. Error reduction uses warp shuffles, not a global atomic for every pixel.
| anchors | FFT tile | exact direct | FFT pipeline | FFT/direct |
|---|---|---|---|---|
| 67108864 | 2048 | 1.58822 ms | 22.4184 ms | 14.12x slower |
| 67108864 | 4096 | 1.42541 ms | 35.9567 ms | 25.23x slower |
| 268435456 | 2048 | 5.70470 ms | 71.2909 ms | 12.50x slower |
All three scans return exactly the known position (51069,419124). The existing outer-sum scanner on the same 67108864-anchor region measured direct 1.38854 ms and tiled 1.15098 ms: FFT/2048 is about 19.48x slower than tiled scanning. Cold/reused cuFFT setup varied widely (about 7--388 ms), so steady-state pipeline time is deliberately separated from setup cost.
Boundary tests also match the direct scanner: a known match starting at offset (123,123) with tile=128 and step=124 crosses both tile seams; a 257-by-259 window places the match at the final valid anchor. The largest observed rounding error in the final 5x5 runs was 1.33514e-5. A preliminary 8x8, 16777216-anchor run returned zero hits in both paths. These are measured agreement checks, not a proof of numerical completeness for every input. The candidate tolerance is 0.25; excessive observed error (>0.1) aborts.
Conclusion: early-exit exact scanning averages only about 1.11 cell checks per anchor for a random all-slime pattern. Computing whole-grid FFT correlations cannot compete here. This does not rule out FFT for non-slime-heavy patterns, repeated convolutions on a reused bitmap, or other correlation workloads. No FFT dependency or executable is retained in the default research build.
To reproduce, save the fenced source as research/slime_fft_cuda.cu, compile
with nvcc -O3 -std=c++17 -o research/slime_fft_cuda research/slime_fft_cuda.cu -lcufft,
and run with --width 8192 --height 8192 --tile 2048 --repeats 5 (or dimensions
16384 by 16384). All memory capacities fail loudly instead of silently
truncating results.
Archived source
// Independent fixed-seed FFT benchmark. Public approach:
// https://github.com/alchemyyy/UltraSlimeChunkFinder/blob/main/slime_fft.cu
#include "slime_solver_core.hpp"
#include <cuda_runtime.h>
#include <cufft.h>
#include <algorithm>
#include <chrono>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
namespace sr = slime_research;
void ck(cudaError_t e) { if (e != cudaSuccess) throw std::runtime_error(cudaGetErrorString(e)); }
void fk(cufftResult e) { if (e != CUFFT_SUCCESS) throw std::runtime_error("cuFFT error " + std::to_string(e)); }
struct Hit {
int x, z;
bool operator<(const Hit& b) const { return x != b.x ? x < b.x : z < b.z; }
bool operator==(const Hit& b) const { return x == b.x && z == b.z; }
};
constexpr unsigned capacity = 1U << 20;
__device__ bool solid(uint64_t seed, int x, int z, int n) {
for (int dx = 0; dx < n; ++dx)
for (int dz = 0; dz < n; ++dz)
if (!sr::exactIsSlimeFromTerm(seed, sr::coordinateTerm(x + dx, z + dz))) return false;
return true;
}
__device__ void emit(Hit* hits, unsigned* count, int x, int z) {
unsigned index = atomicAdd(count, 1U);
if (index < capacity) hits[index] = {x, z};
}
__global__ void direct(uint64_t seed, int x0, int z0, int w, int h, int n,
Hit* hits, unsigned* count) {
for (uint64_t i = uint64_t(blockIdx.x) * blockDim.x + threadIdx.x;
i < uint64_t(w) * h; i += uint64_t(gridDim.x) * blockDim.x) {
int x = x0 + i / h, z = z0 + i % h;
if (solid(seed, x, z, n)) emit(hits, count, x, z);
}
}
__global__ void bitmap(uint64_t seed, float* image, int tile, int x0, int z0,
int validX, int validZ) {
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < tile * tile;
i += gridDim.x * blockDim.x) {
int dx = i / tile, dz = i % tile;
image[i] = dx < validX && dz < validZ &&
sr::exactIsSlimeFromTerm(seed, sr::coordinateTerm(x0 + dx, z0 + dz)) ? 1.f : 0.f;
}
}
__global__ void multiply(cufftComplex* a, const cufftComplex* b, int count) {
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < count; i += gridDim.x * blockDim.x) {
float ar = a[i].x, ai = a[i].y;
a[i] = {ar * b[i].x - ai * b[i].y, ar * b[i].y + ai * b[i].x};
}
}
__global__ void peaks(uint64_t seed, const float* correlation, int tile,
int x0, int z0, int w, int h, int n, Hit* hits,
unsigned* count, unsigned* errorBits) {
float maxError = 0;
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < w * h; i += gridDim.x * blockDim.x) {
int dx = i / h, dz = i % h;
float value = correlation[dx * tile + dz] / float(tile * tile);
maxError = fmaxf(maxError, fabsf(value - roundf(value)));
if (value > float(n * n) - .25f && solid(seed, x0 + dx, z0 + dz, n))
emit(hits, count, x0 + dx, z0 + dz);
}
for (int offset = 16; offset; offset /= 2)
maxError = fmaxf(maxError, __shfl_down_sync(0xffffffffU, maxError, offset));
if ((threadIdx.x & 31) == 0) atomicMax(errorBits, __float_as_uint(maxError));
}
struct Args {
uint64_t seed = 61692151016558ULL;
int x0 = 46973, z0 = 415028, width = 8192, height = 8192, size = 5, tile = 2048, repeats = 3;
};
int main(int argc, char** argv) {
try {
Args a;
for (int i = 1; i < argc; ++i) {
std::string key = argv[i];
if (key == "--help") {
std::cout << "--seed N --x0 N --z0 N --width N --height N -n N --tile N --repeats N\n";
return 0;
}
if (++i == argc) throw std::runtime_error("missing argument");
std::string value = argv[i];
if (key == "--seed") a.seed = std::stoull(value) & sr::kMask48;
else if (key == "--x0") a.x0 = std::stoi(value);
else if (key == "--z0") a.z0 = std::stoi(value);
else if (key == "--width") a.width = std::stoi(value);
else if (key == "--height") a.height = std::stoi(value);
else if (key == "-n") a.size = std::stoi(value);
else if (key == "--tile") a.tile = std::stoi(value);
else if (key == "--repeats") a.repeats = std::stoi(value);
else throw std::runtime_error("unknown argument " + key);
}
if (a.size < 1 || a.size > 64 || a.tile < a.size || a.tile > 8192 ||
(a.tile & (a.tile - 1)) || a.width < 1 || a.height < 1 || a.repeats < 1 ||
int64_t(a.x0) + a.width + a.size - 2 > INT32_MAX ||
int64_t(a.z0) + a.height + a.size - 2 > INT32_MAX)
throw std::runtime_error("invalid dimensions or int32 coordinate overflow");
ck(cudaFree(nullptr));
auto setupStart = std::chrono::steady_clock::now();
int pixels = a.tile * a.tile, frequencies = a.tile * (a.tile / 2 + 1);
float *image, *correlation;
cufftComplex *spectrum, *pattern;
Hit* hits; unsigned *count, *errorBits;
ck(cudaMalloc(&image, pixels * sizeof(float)));
ck(cudaMalloc(&correlation, pixels * sizeof(float)));
ck(cudaMalloc(&spectrum, frequencies * sizeof(cufftComplex)));
ck(cudaMalloc(&pattern, frequencies * sizeof(cufftComplex)));
ck(cudaMalloc(&hits, capacity * sizeof(Hit)));
ck(cudaMalloc(&count, sizeof(unsigned)));
ck(cudaMalloc(&errorBits, sizeof(unsigned)));
cufftHandle forward, inverse;
fk(cufftPlan2d(&forward, a.tile, a.tile, CUFFT_R2C));
fk(cufftPlan2d(&inverse, a.tile, a.tile, CUFFT_C2R));
std::vector<float> stencil(pixels, 0.f);
for (int dx = 0; dx < a.size; ++dx)
for (int dz = 0; dz < a.size; ++dz)
stencil[((a.tile - dx) % a.tile) * a.tile + (a.tile - dz) % a.tile] = 1.f;
ck(cudaMemcpy(image, stencil.data(), pixels * sizeof(float), cudaMemcpyHostToDevice));
fk(cufftExecR2C(forward, image, pattern));
ck(cudaDeviceSynchronize());
double setupMs = std::chrono::duration<double, std::milli>(std::chrono::steady_clock::now() - setupStart).count();
auto readHits = [&]() {
unsigned total;
ck(cudaMemcpy(&total, count, sizeof(total), cudaMemcpyDeviceToHost));
if (total > capacity) throw std::runtime_error("hit capacity exceeded; reduce window");
std::vector<Hit> result(total);
if (total) ck(cudaMemcpy(result.data(), hits, total * sizeof(Hit), cudaMemcpyDeviceToHost));
std::sort(result.begin(), result.end());
return result;
};
cudaEvent_t begin, end;
ck(cudaEventCreate(&begin)); ck(cudaEventCreate(&end));
std::vector<float> directMs, fftMs;
std::vector<Hit> reference;
float maxError = 0;
for (int repeat = -1; repeat < a.repeats; ++repeat) {
ck(cudaMemset(count, 0, sizeof(unsigned)));
ck(cudaEventRecord(begin));
int blocks = static_cast<int>(std::min<uint64_t>(65535, (uint64_t(a.width) * a.height + 255) / 256));
direct<<<blocks, 256>>>(a.seed, a.x0, a.z0, a.width, a.height, a.size, hits, count);
ck(cudaGetLastError()); ck(cudaEventRecord(end)); ck(cudaEventSynchronize(end));
float elapsed; ck(cudaEventElapsedTime(&elapsed, begin, end));
if (repeat >= 0) directMs.push_back(elapsed);
reference = readHits();
ck(cudaMemset(count, 0, sizeof(unsigned)));
ck(cudaMemset(errorBits, 0, sizeof(unsigned)));
ck(cudaEventRecord(begin));
int step = a.tile - a.size + 1;
for (int x = 0; x < a.width; x += step) for (int z = 0; z < a.height; z += step) {
int w = std::min(step, a.width - x), h = std::min(step, a.height - z);
bitmap<<<std::min(4096, (pixels + 255) / 256), 256>>>(a.seed, image, a.tile, a.x0 + x, a.z0 + z, w + a.size - 1, h + a.size - 1);
fk(cufftExecR2C(forward, image, spectrum));
multiply<<<std::min(4096, (frequencies + 255) / 256), 256>>>(spectrum, pattern, frequencies);
fk(cufftExecC2R(inverse, spectrum, correlation));
peaks<<<std::min(4096, (w * h + 255) / 256), 256>>>(a.seed, correlation, a.tile, a.x0 + x, a.z0 + z, w, h, a.size, hits, count, errorBits);
}
ck(cudaGetLastError()); ck(cudaEventRecord(end)); ck(cudaEventSynchronize(end));
ck(cudaEventElapsedTime(&elapsed, begin, end));
if (repeat >= 0) fftMs.push_back(elapsed);
float error; ck(cudaMemcpy(&error, errorBits, sizeof(error), cudaMemcpyDeviceToHost));
maxError = std::max(maxError, error);
if (error > .1f) throw std::runtime_error("FFT numeric error too large");
if (readHits() != reference) throw std::runtime_error("FFT disagrees with exact direct scan");
}
std::sort(directMs.begin(), directMs.end()); std::sort(fftMs.begin(), fftMs.end());
float d = directMs[directMs.size() / 2], f = fftMs[fftMs.size() / 2];
std::cout << "positions=" << uint64_t(a.width) * a.height << " size=" << a.size
<< " tile=" << a.tile << " hits=" << reference.size()
<< " direct_ms=" << d << " fft_ms=" << f << " speedup=" << d / f
<< " fft_setup_ms=" << setupMs << " max_rounding_error=" << maxError << '\n';
for (auto hit : reference) std::cout << "hit " << hit.x << ' ' << hit.z << '\n';
cufftDestroy(forward); cufftDestroy(inverse);
cudaEventDestroy(begin); cudaEventDestroy(end);
cudaFree(image); cudaFree(correlation); cudaFree(spectrum); cudaFree(pattern);
cudaFree(hits); cudaFree(count); cudaFree(errorBits);
return 0;
} catch (const std::exception& error) {
std::cerr << "error: " << error.what() << '\n'; return 1;
}
}