233 lines
12 KiB
Markdown
233 lines
12 KiB
Markdown
# 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](https://github.com/alchemyyy/UltraSlimeChunkFinder/blob/main/slime_fft.cu).
|
|
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
|
|
|
|
```cpp
|
|
// 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;
|
|
}
|
|
}
|
|
```
|