// ============================================================================
// SHARED-NOISE PARITY — is the .inc the same function on both sides?
//
// src/shared/blender_hash.inc and src/shared/blender_perlin.inc are compiled
// TWICE: as C++ through blender_noise_cpp_shim.h, and as GLSL through the
// generated headers the engine hands the driver. There is only one
// implementation of each algorithm, which is the entire point -- but "one
// text" is not "one function" until something compiles both and compares.
// What survives having one text is the set of places the two LANGUAGES read
// the same characters differently: arithmetic vs logical shift, signed
// overflow, the bit cast, truncated vs floored modulo, and float literals
// that are double in C++.
//
// NEITHER SIDE IS WRITTEN FOR THE TEST. The CPU calls the shim; the GLSL is
// EXTRACTED from the generated header at run time, never retyped. The erosion
// filter's lesson was that a transcription used AS the reference can only ever
// agree with itself.
//
// TWO CLASSES OF CASE, AND THEY ARE NOT COMPARED THE SAME WAY.
//
// EXACT the hash layer. Integer arithmetic, which GLSL requires to wrap and
// C++ defines as modulo 2^32, so the two languages compute the same
// integers or they do not. Its float results come from an integer
// conversion, a division and a multiply, all correctly rounded by
// IEEE 754 on both sides. No epsilon, and a near-miss is a finding.
//
// NEAR the perlin layer. Up to eight octaves of multiply-add, where a GPU
// may contract a*b+c into one FMA and the CPU rounds twice. That is a
// difference in the last bits, not in the function, so these carry a
// relative tolerance -- and the run PRINTS the worst error it saw and
// how many were bit-identical anyway, because a tolerance nobody
// measures against is where a real divergence hides.
//
// ---------------------------------------------------------------------------
// ⚠ WHY THIS IS DATA-DRIVEN, AND WHY THE MEMORY GUARD EXISTS.
//
// The first version of this gate did what tests/test_emitted_op_parity.cpp
// does: emit one line of GLSL per case into a single main(), on the reasoning
// that "246 pipelines to answer 246 questions is a slow gate nobody runs".
// That is sound for 658 cheap Math cases. It does not survive 1300 perlin
// cases, each an eight-octave chain: the shader becomes a combinatorial
// blow-up that is inlined and register-allocated IN THIS PROCESS, and on
// 2026-09-05 it twice took a 31 GB machine to global OOM -- 13.9 GB, then
// 16.1 GB resident -- killing the session around it. With the inputs as
// literals glslc also tried to constant-fold every chain, which is worse than
// slow: a folded case is the COMPILER's arithmetic, not the GPU's, so the gate
// quietly stops involving the machine it exists to test.
//
// So: ONE record per case in a buffer, ONE dispatcher in the shader, cases run
// across invocations. Shader size is O(number of node families), not O(number
// of cases), and every input is opaque to the optimiser. Adding a thousand
// cases now costs a thousand records and not one line of GLSL.
//
// The opcode constants are EMITTED INTO THE SHADER from the C++ enum below, so
// the two dispatchers cannot drift apart by editing one of them.
//
// And the process caps its own address space (see MEM_LIMIT_BYTES): if this
// ever blows up again it dies alone with a clear message instead of taking the
// machine with it.
// ---------------------------------------------------------------------------
//
// Build (no engine, no godot-cpp, no scons):
// g++ -O2 -std=c++17 -I src/shared \
// -o /tmp/noise_parity tests/test_shared_noise_parity.cpp -lvulkan
//
// Run:
// /tmp/noise_parity # default subjects
// /tmp/noise_parity <file> <marker> [...] # override, in dependency order
//
// Exit status is 0 when every case agrees.
// ============================================================================
#include <vulkan/vulkan.h>
#include <sys/resource.h>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
#include "blender_noise_cpp_shim.h"
// NOT `using`d. dcerosion has its own vec2/vec3/vec4, floor, max, dot and
// mix, and so does bnoisegl -- two shims for two .inc families, deliberately
// separate. Every erosion call below is qualified, and the types are built
// explicitly, so nothing is resolved by luck.
#include "erosion_cpp_shim.h"
using namespace godot::bnoisegl;
// hash FIRST: perlin's lattice gradients are that hash, and GLSL has no
// #include to say so. Same order the C++ shim uses and the emitter emits.
static const char *SUBJECT_FILES[] = {
"src/gpu/blender_hash_glsl.gen.h",
"src/gpu/blender_perlin_glsl.gen.h",
"src/gpu/blender_voronoi_glsl.gen.h",
"src/gpu/blender_ramp_blend_glsl.gen.h",
"src/gpu/proc_simplex_glsl.gen.h",
// THE SIXTH, added 2026-09-06. erosion_filter.inc was the one dual-compiled
// source not among the subjects -- the oldest of them, the one three other
// systems evaluate, and the only one written before the f-suffix rule. It
// carries 119 bare float literals and not one suffixed, so on the C++ side
// every expression containing one is evaluated in DOUBLE and narrowed at the
// end, where the shader computes the whole chain in float. `t*t*t*(t*(t*6.0
// - 15.0) + 10.0)` is that hazard in one line, and it is the quintic fade at
// the bottom of every octave of every erosion evaluation.
//
// It is self-contained: pn_max_octaves appears only in its comments, and the
// pn_ero_* globals belong to erosion_wrapper_glsl.h, not to this text.
//
// ⚠ THE HAZARD'S EFFECT WAS MEASURED, AND IT DOES NOT JUSTIFY SUFFIXING.
// hyper-games #65 asked for the effect, not the hazard. Both texts were
// compiled at once -- the shipped .inc and an f-suffixed twin, in two
// namespaces -- and asked the same questions over a dense sweep:
//
// erosion_hash 29282 cases 0 differ
// grad_noised_2d 43923 cases 0 differ
// fractal_noise_2d 43923 cases 0 differ
// phacelle_noise 58564 cases 4276 differ (7.3%) worst 9.2e-07
// erosion_filter 8405 cases 1713 differ (20.4%) worst 3.0e-07
// erosion_node_eval_full 3362 cases 581 differ (17.3%) worst 2.4e-07
//
// The three zeros are the finding. `vec2 u = f*f*f*(f*(f*6.0-15.0)+10.0)` --
// the quintic fade, which I predicted would be the sharpest site -- is
// IMMUNE, because `f*6.0` resolves to the shim's `operator*(vec2, float)`
// and the literal narrows AT THE CALL. Only genuinely scalar chains
// (`t - 0.5 * smoothing`) promote, which is why the three that differ are
// the three with scalar arithmetic in them.
//
// And suffixing does not move the CPU toward the GPU. With the twin in
// place the whole gate went from 3862/5371 bit-identical to 3861 -- one
// case WORSE. At this magnitude the GPU's FMA contraction dominates and it
// is not systematically on either side of the double promotion. So the
// literals are left alone: rewriting 119 of them in a filter that
// dc-terrain and hyper-tps both evaluate, to move numbers by 3e-07 in a
// direction the measurement does not endorse, is churn with a surface
// change attached. The class is not removed by construction here; it is
// covered by this gate instead.
"src/gpu/erosion_filter_glsl.gen.h",
};
static const char *SUBJECT_MARKS[] = { "BLENDER_HASH_GLSL", "BLENDER_PERLIN_GLSL",
"BLENDER_VORONOI_GLSL", "BLENDER_RAMP_BLEND_GLSL", "PROC_SIMPLEX_GLSL",
"EROSION_FILTER_GLSL" };
static const int SUBJECT_N = 6;
static const double TOL_REL = 1e-5;
static const double TOL_ABS = 1e-6;
// 4 GiB. A healthy run of this gate peaks around a tenth of that; the headroom
// is for the Vulkan driver's own allocations, which are not small and are not
// ours to predict. The point is not to be tight, it is to be BOUNDED.
static const rlim_t MEM_LIMIT_BYTES = (rlim_t)4 * 1024 * 1024 * 1024;
static void cap_own_memory() {
struct rlimit rl;
if (getrlimit(RLIMIT_AS, &rl) != 0) {
fprintf(stderr, "WARNING: getrlimit failed; running with NO memory cap\n");
return;
}
if (rl.rlim_cur != RLIM_INFINITY && rl.rlim_cur <= MEM_LIMIT_BYTES) {
return; // already at least as strict as we want
}
const rlim_t want = (rl.rlim_max == RLIM_INFINITY || rl.rlim_max > MEM_LIMIT_BYTES)
? MEM_LIMIT_BYTES
: rl.rlim_max;
rl.rlim_cur = want;
if (setrlimit(RLIMIT_AS, &rl) != 0) {
fprintf(stderr, "WARNING: setrlimit failed; running with NO memory cap\n");
return;
}
printf("memcap : %.1f GiB address space (this gate OOM-killed the machine "
"twice before it was data-driven)\n", (double)want / (1024.0 * 1024.0 * 1024.0));
}
static void die(const char *what, VkResult r) {
fprintf(stderr, "FATAL: %s -> VkResult %d\n", what, (int)r);
// ALWAYS name the cap, not only on VK_ERROR_OUT_OF_HOST_MEMORY. Starving
// this process shows up as VK_ERROR_INITIALIZATION_FAILED (-3) from
// vkEnumeratePhysicalDevices long before anyone gets to allocate anything,
// and a reader who does not know a cap is in force will go looking at the
// driver. Measured: at `ulimit -v 200000` that is exactly what happens.
struct rlimit rl;
if (getrlimit(RLIMIT_AS, &rl) == 0 && rl.rlim_cur != RLIM_INFINITY) {
fprintf(stderr, " This process is capped at %.1f GiB of address space.\n"
" The cap is deliberate: without it this gate twice took a\n"
" 31 GB machine to global OOM. If it is being hit, check\n"
" whether the shader has grown back into something whose\n"
" size depends on the case count -- fix that, do not raise\n"
" the cap.\n",
(double)rl.rlim_cur / (1024.0 * 1024.0 * 1024.0));
}
exit(1);
}
#define VK(x) do { VkResult _r = (x); if (_r != VK_SUCCESS) die(#x, _r); } while (0)
static std::string read_file(const std::string &path) {
std::ifstream f(path);
if (!f) { fprintf(stderr, "FATAL: cannot read %s\n", path.c_str()); exit(1); }
std::stringstream ss;
ss << f.rdbuf();
return ss.str();
}
// Pull the GLSL out of a C++ raw string literal that follows `marker`, reading
// the delimiter off the text rather than assuming R"( -- the generated headers
// carry one because their payload is GLSL full of parentheses.
static std::string extract_raw_string(const std::string &src, const std::string &marker) {
size_t m = src.find(marker);
if (m == std::string::npos) {
fprintf(stderr, "FATAL: marker '%s' not found\n", marker.c_str());
exit(1);
}
size_t rq = src.find("R\"", m);
if (rq == std::string::npos) {
fprintf(stderr, "FATAL: no raw string after '%s'\n", marker.c_str());
exit(1);
}
size_t lp = src.find('(', rq);
if (lp == std::string::npos) {
fprintf(stderr, "FATAL: malformed raw string after '%s'\n", marker.c_str());
exit(1);
}
const std::string delim = src.substr(rq + 2, lp - (rq + 2));
const std::string term = ")" + delim + "\"";
size_t close = src.find(term, lp + 1);
if (close == std::string::npos) {
fprintf(stderr, "FATAL: unterminated raw string (delimiter '%s') after '%s'\n",
delim.c_str(), marker.c_str());
exit(1);
}
return src.substr(lp + 1, close - (lp + 1));
}
static std::vector<uint32_t> compile_glsl(const std::string &src) {
char in_path[] = "/tmp/proc_noise_parity_XXXXXX";
int fd = mkstemp(in_path);
if (fd < 0) { perror("mkstemp"); exit(1); }
FILE *f = fdopen(fd, "w");
fwrite(src.data(), 1, src.size(), f);
fclose(f);
std::string spv = std::string(in_path) + ".spv";
std::string cmd = "glslc -fshader-stage=comp -o " + spv + " " + in_path;
if (system(cmd.c_str()) != 0) {
fprintf(stderr, "FATAL: glslc failed (source left at %s)\n", in_path);
exit(1);
}
FILE *sf = fopen(spv.c_str(), "rb");
fseek(sf, 0, SEEK_END); long n = ftell(sf); fseek(sf, 0, SEEK_SET);
std::vector<uint32_t> out(n / 4);
if (fread(out.data(), 1, n, sf) != (size_t)n) { fprintf(stderr, "FATAL: short read\n"); exit(1); }
fclose(sf);
remove(in_path); remove(spv.c_str());
return out;
}
struct Ctx {
VkInstance instance = VK_NULL_HANDLE;
VkPhysicalDevice phys = VK_NULL_HANDLE;
VkDevice dev = VK_NULL_HANDLE;
VkQueue queue = VK_NULL_HANDLE;
uint32_t qfam = 0;
VkCommandPool cmdpool = VK_NULL_HANDLE;
VkPhysicalDeviceMemoryProperties memprops{};
char devname[256] = {0};
};
static uint32_t find_mem(Ctx &c, uint32_t bits, VkMemoryPropertyFlags want) {
for (uint32_t i = 0; i < c.memprops.memoryTypeCount; i++) {
if ((bits & (1u << i)) && (c.memprops.memoryTypes[i].propertyFlags & want) == want) {
return i;
}
}
fprintf(stderr, "FATAL: no memory type with flags 0x%x\n", want);
exit(1);
}
static void init_vulkan(Ctx &c) {
VkApplicationInfo app{VK_STRUCTURE_TYPE_APPLICATION_INFO};
app.pApplicationName = "proc-nodes shared-noise parity";
app.apiVersion = VK_API_VERSION_1_1;
VkInstanceCreateInfo ii{VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO};
ii.pApplicationInfo = &app;
VK(vkCreateInstance(&ii, nullptr, &c.instance));
uint32_t n = 0;
VK(vkEnumeratePhysicalDevices(c.instance, &n, nullptr));
if (!n) { fprintf(stderr, "FATAL: no Vulkan device\n"); exit(1); }
std::vector<VkPhysicalDevice> devs(n);
VK(vkEnumeratePhysicalDevices(c.instance, &n, devs.data()));
c.phys = devs[0];
for (size_t i = 0; i < devs.size(); i++) {
VkPhysicalDeviceProperties p;
vkGetPhysicalDeviceProperties(devs[i], &p);
if (p.deviceType == VK_PHYSICAL_DEVICE_TYPE_DISCRETE_GPU) { c.phys = devs[i]; break; }
}
VkPhysicalDeviceProperties props;
vkGetPhysicalDeviceProperties(c.phys, &props);
snprintf(c.devname, sizeof(c.devname), "%s", props.deviceName);
vkGetPhysicalDeviceMemoryProperties(c.phys, &c.memprops);
uint32_t qn = 0;
vkGetPhysicalDeviceQueueFamilyProperties(c.phys, &qn, nullptr);
std::vector<VkQueueFamilyProperties> qs(qn);
vkGetPhysicalDeviceQueueFamilyProperties(c.phys, &qn, qs.data());
bool found = false;
for (uint32_t i = 0; i < qn; i++) {
if (qs[i].queueFlags & VK_QUEUE_COMPUTE_BIT) { c.qfam = i; found = true; break; }
}
if (!found) { fprintf(stderr, "FATAL: no compute queue\n"); exit(1); }
float prio = 1.0f;
VkDeviceQueueCreateInfo qi{VK_STRUCTURE_TYPE_DEVICE_QUEUE_CREATE_INFO};
qi.queueFamilyIndex = c.qfam; qi.queueCount = 1; qi.pQueuePriorities = &prio;
VkDeviceCreateInfo di{VK_STRUCTURE_TYPE_DEVICE_CREATE_INFO};
di.queueCreateInfoCount = 1; di.pQueueCreateInfos = &qi;
VK(vkCreateDevice(c.phys, &di, nullptr, &c.dev));
vkGetDeviceQueue(c.dev, c.qfam, 0, &c.queue);
VkCommandPoolCreateInfo ci{VK_STRUCTURE_TYPE_COMMAND_POOL_CREATE_INFO};
ci.queueFamilyIndex = c.qfam;
ci.flags = VK_COMMAND_POOL_CREATE_RESET_COMMAND_BUFFER_BIT;
VK(vkCreateCommandPool(c.dev, &ci, nullptr, &c.cmdpool));
}
// ============================================================================
// THE CASE RECORD
//
// One fixed-size record per case, read by both dispatchers. Every numeric
// input travels here, so the shader compiler cannot fold it and the GPU has to
// actually run. Only the opcode and the dimension pick a code path.
// ============================================================================
// 24 rather than 16 since Voronoi arrived: its parameter block is eleven
// fields and every one of them changes the answer. 56 since the erosion filter
// arrived: erosion_node_eval_full takes NINETEEN parameters, thirty floats once
// the vec2s and vec4s are counted out, and every one of them changes the
// surface. Widening the record costs 1300 x 32 more words -- 166 KB -- and the
// alternative is reusing slots that already mean something else, which is how a
// gate comes to compare two things neither side meant.
static const uint32_t STRIDE = 56;
enum Op {
OP_HASH1 = 1, OP_HASH2, OP_HASH3, OP_HASH4,
OP_UINT_TO_FLOAT_01,
OP_HASH_TO_FLOAT1, OP_HASH_TO_FLOAT2, OP_HASH_TO_FLOAT3, OP_HASH_TO_FLOAT4,
OP_HASH_FLOAT, // dim 1..4
OP_HFF, // dim 1..4
OP_HFF2, // dim 2..4, comp 0..1
OP_HFF3, // dim 1..4, comp 0..2
OP_HFF4, // comp 0..3
OP_WRAP_MUL, OP_WRAP_ADD,
OP_INT_TO_FLOAT_01,
OP_PCG2D, OP_PCG3D, OP_PCG4D, // comp
OP_INT_TO_FLOAT2, OP_INT_TO_FLOAT3_I3, OP_INT_TO_FLOAT4,
OP_INT_TO_FLOAT3_I2, OP_INT_TO_FLOAT3_I4,
OP_PERLIN_NOISE, // dim
OP_PERLIN_SIGNED, // dim
OP_PERLIN_SELECT, // dim, type, flags
OP_PERLIN_DISTORTED, // dim, type, flags
OP_PERLIN_FLOAT3, // dim, comp, type, flags
OP_VORONOI_XFX, // dim, comp (0=dist,1..3=color,4..7=position)
OP_VORONOI_EDGE, // dim
OP_VORONOI_RADIUS, // dim
OP_RAMP_BLEND, // type = blend mode, comp = channel
OP_SIMPLEX, // the raw Ashima lattice
OP_SIMPLEX_FBM, // type = octaves
// erosion_filter.inc, bottom-up. The layering is the one the Blender port
// was built in, and for the same reason: a divergence in erosion_hash is
// unreadable once it has been through eight octaves and a fade.
OP_ERO_HASH, // comp 0..1 -- integer lattice, EXACT
OP_ERO_GRAD_NOISED, // comp 0..2 -- height + ANALYTIC derivative
OP_ERO_FRACTAL, // comp 0..2 -- the octave loop
OP_ERO_CLAMP01,
OP_ERO_POW_INV,
OP_ERO_EASE_OUT,
OP_ERO_SMOOTH_START,
OP_ERO_NORMALIZE_SAFE, // comp 0..1
OP_ERO_PHACELLE, // comp 0..3 -- the fade-approach primitive
OP_ERO_FILTER, // comp 0..4 -- delta.xyz, magnitude, ridge
OP_ERO_EVAL_FULL, // comp 0..1 -- the wrapper ErosionNoiseNode calls
OP_ERO_EVAL_GRAD, // comp 0..3 -- same, keeping the slope
OP__COUNT
};
// Emitted into the shader so the two dispatchers cannot drift.
struct OpName { Op op; const char *name; };
static const OpName OP_NAMES[] = {
{ OP_HASH1, "OP_HASH1" }, { OP_HASH2, "OP_HASH2" }, { OP_HASH3, "OP_HASH3" },
{ OP_HASH4, "OP_HASH4" }, { OP_UINT_TO_FLOAT_01, "OP_UINT_TO_FLOAT_01" },
{ OP_HASH_TO_FLOAT1, "OP_HASH_TO_FLOAT1" }, { OP_HASH_TO_FLOAT2, "OP_HASH_TO_FLOAT2" },
{ OP_HASH_TO_FLOAT3, "OP_HASH_TO_FLOAT3" }, { OP_HASH_TO_FLOAT4, "OP_HASH_TO_FLOAT4" },
{ OP_HASH_FLOAT, "OP_HASH_FLOAT" }, { OP_HFF, "OP_HFF" }, { OP_HFF2, "OP_HFF2" },
{ OP_HFF3, "OP_HFF3" }, { OP_HFF4, "OP_HFF4" },
{ OP_WRAP_MUL, "OP_WRAP_MUL" }, { OP_WRAP_ADD, "OP_WRAP_ADD" },
{ OP_INT_TO_FLOAT_01, "OP_INT_TO_FLOAT_01" },
{ OP_PCG2D, "OP_PCG2D" }, { OP_PCG3D, "OP_PCG3D" }, { OP_PCG4D, "OP_PCG4D" },
{ OP_INT_TO_FLOAT2, "OP_INT_TO_FLOAT2" }, { OP_INT_TO_FLOAT3_I3, "OP_INT_TO_FLOAT3_I3" },
{ OP_INT_TO_FLOAT4, "OP_INT_TO_FLOAT4" }, { OP_INT_TO_FLOAT3_I2, "OP_INT_TO_FLOAT3_I2" },
{ OP_INT_TO_FLOAT3_I4, "OP_INT_TO_FLOAT3_I4" },
{ OP_PERLIN_NOISE, "OP_PERLIN_NOISE" }, { OP_PERLIN_SIGNED, "OP_PERLIN_SIGNED" },
{ OP_PERLIN_SELECT, "OP_PERLIN_SELECT" }, { OP_PERLIN_DISTORTED, "OP_PERLIN_DISTORTED" },
{ OP_PERLIN_FLOAT3, "OP_PERLIN_FLOAT3" },
{ OP_VORONOI_XFX, "OP_VORONOI_XFX" }, { OP_VORONOI_EDGE, "OP_VORONOI_EDGE" },
{ OP_VORONOI_RADIUS, "OP_VORONOI_RADIUS" }, { OP_RAMP_BLEND, "OP_RAMP_BLEND" },
{ OP_SIMPLEX, "OP_SIMPLEX" }, { OP_SIMPLEX_FBM, "OP_SIMPLEX_FBM" },
{ OP_ERO_HASH, "OP_ERO_HASH" }, { OP_ERO_GRAD_NOISED, "OP_ERO_GRAD_NOISED" },
{ OP_ERO_FRACTAL, "OP_ERO_FRACTAL" }, { OP_ERO_CLAMP01, "OP_ERO_CLAMP01" },
{ OP_ERO_POW_INV, "OP_ERO_POW_INV" }, { OP_ERO_EASE_OUT, "OP_ERO_EASE_OUT" },
{ OP_ERO_SMOOTH_START, "OP_ERO_SMOOTH_START" },
{ OP_ERO_NORMALIZE_SAFE, "OP_ERO_NORMALIZE_SAFE" },
{ OP_ERO_PHACELLE, "OP_ERO_PHACELLE" }, { OP_ERO_FILTER, "OP_ERO_FILTER" },
{ OP_ERO_EVAL_FULL, "OP_ERO_EVAL_FULL" }, { OP_ERO_EVAL_GRAD, "OP_ERO_EVAL_GRAD" },
};
struct Rec {
uint32_t w[STRIDE];
std::string name;
bool exact;
double tol_abs = 0.0;
uint32_t cpu;
};
static uint32_t bits_of(float f) { uint32_t u; memcpy(&u, &f, sizeof(u)); return u; }
static float float_of(uint32_t u) { float f; memcpy(&f, &u, sizeof(f)); return f; }
// Word slots, named once so both sides agree by construction.
enum { W_OP = 0, W_DIM = 1, W_TYPE = 2, W_FLAGS = 3, W_COMP = 4,
W_A = 5, W_B = 6, W_C = 7, W_D = 8,
W_DETAIL = 9, W_ROUGH = 10, W_LAC = 11, W_OFF = 12, W_GAIN = 13, W_DIST = 14,
// Voronoi's parameter block
W_SCALE = 15, W_SMOOTH = 16, W_EXPONENT = 17, W_RANDOMNESS = 18, W_MAXDIST = 19,
W_METRIC = 20,
// The erosion block. Voronoi's slots are spelled `gF(b, 15u)` by hand in the
// GLSL dispatcher, which is the drift this file warns about everywhere else,
// so these are EMITTED into the shader from the enum (see SLOT_NAMES) and
// both sides name them.
W_E_PX = 24, W_E_PY = 25,
W_E_NX = 26, W_E_NY = 27, W_E_NZ = 28,
W_E_S0 = 29, // 16 scalars, W_E_S0 + k
W_E_R0 = 45, // rounding vec4
W_E_O0 = 49, // onset vec4
W_E_A0 = 53, // assumed slope vec2
W_E_OCT = 55 };
// Emitted into the shader beside the opcodes, same reason.
struct SlotName { int slot; const char *name; };
static const SlotName SLOT_NAMES[] = {
{ W_E_PX, "W_E_PX" }, { W_E_PY, "W_E_PY" },
{ W_E_NX, "W_E_NX" }, { W_E_NY, "W_E_NY" }, { W_E_NZ, "W_E_NZ" },
{ W_E_S0, "W_E_S0" }, { W_E_R0, "W_E_R0" }, { W_E_O0, "W_E_O0" },
{ W_E_A0, "W_E_A0" }, { W_E_OCT, "W_E_OCT" },
};
// Filled from the record, all eleven fields -- a GLSL struct has no defaults
// and a field left unset is garbage, not a sensible value.
static BvParams params_from(const uint32_t *w) {
BvParams p;
p.scale = float_of(w[W_SCALE]);
p.detail = float_of(w[W_DETAIL]);
p.roughness = float_of(w[W_ROUGH]);
p.lacunarity = float_of(w[W_LAC]);
p.smoothness = float_of(w[W_SMOOTH]);
p.exponent = float_of(w[W_EXPONENT]);
p.randomness = float_of(w[W_RANDOMNESS]);
p.max_distance = float_of(w[W_MAXDIST]);
p.normalize = (w[W_FLAGS] & 1u) != 0u;
p.feature = (int)w[W_TYPE];
p.metric = (int)w[W_METRIC];
return p;
}
// The erosion block's accessors. es(w, k) is scalar k; e2/e4 build the vector
// bundles at a named base. Spelled once so the arms below read like the .inc's
// own signatures.
static float es(const uint32_t *w, int k) { return float_of(w[W_E_S0 + k]); }
static dcerosion::vec2 e2(const uint32_t *w, int base) {
return dcerosion::vec2(float_of(w[base]), float_of(w[base + 1]));
}
static dcerosion::vec4 e4(const uint32_t *w, int base) {
return dcerosion::vec4(float_of(w[base]), float_of(w[base + 1]),
float_of(w[base + 2]), float_of(w[base + 3]));
}
// --- the CPU half of the dispatcher -----------------------------------------
static uint32_t eval_cpu(const uint32_t *w) {
const uint32_t op = w[W_OP], dim = w[W_DIM], comp = w[W_COMP];
const int type = (int)w[W_TYPE];
const bool nrm = (w[W_FLAGS] & 1u) != 0u;
const uint32_t ua = w[W_A], ub = w[W_B], uc = w[W_C], ud = w[W_D];
const int32_t ia = (int32_t)ua, ib = (int32_t)ub, ic = (int32_t)uc, id = (int32_t)ud;
const float x = float_of(ua), y = float_of(ub), z = float_of(uc), wv = float_of(ud);
const float detail = float_of(w[W_DETAIL]), rough = float_of(w[W_ROUGH]);
const float lac = float_of(w[W_LAC]), off = float_of(w[W_OFF]);
const float gain = float_of(w[W_GAIN]), dist = float_of(w[W_DIST]);
const vec2 v2(x, y); const vec3 v3(x, y, z); const vec4 v4(x, y, z, wv);
const ivec2 i2(ia, ib); const ivec3 i3(ia, ib, ic); const ivec4 i4(ia, ib, ic, id);
// The erosion block, read by slot name on both sides.
const float ex = float_of(w[W_E_PX]), ey = float_of(w[W_E_PY]);
const float enx = float_of(w[W_E_NX]), eny = float_of(w[W_E_NY]),
enz = float_of(w[W_E_NZ]);
switch (op) {
case OP_HASH1: return bn_hash(ua);
case OP_HASH2: return bn_hash(ua, ub);
case OP_HASH3: return bn_hash(ua, ub, uc);
case OP_HASH4: return bn_hash(ua, ub, uc, ud);
case OP_UINT_TO_FLOAT_01: return bits_of(bn_uint_to_float_01(ua));
case OP_HASH_TO_FLOAT1: return bits_of(bn_hash_to_float(ua));
case OP_HASH_TO_FLOAT2: return bits_of(bn_hash_to_float(ua, ub));
case OP_HASH_TO_FLOAT3: return bits_of(bn_hash_to_float(ua, ub, uc));
case OP_HASH_TO_FLOAT4: return bits_of(bn_hash_to_float(ua, ub, uc, ud));
case OP_HASH_FLOAT:
if (dim == 1) return bn_hash_float(x);
if (dim == 2) return bn_hash_float(v2);
if (dim == 3) return bn_hash_float(v3);
return bn_hash_float(v4);
case OP_HFF:
if (dim == 1) return bits_of(bn_hash_float_to_float(x));
if (dim == 2) return bits_of(bn_hash_float_to_float(v2));
if (dim == 3) return bits_of(bn_hash_float_to_float(v3));
return bits_of(bn_hash_float_to_float(v4));
case OP_HFF2: {
vec2 r = (dim == 2) ? bn_hash_float_to_float2(v2)
: (dim == 3) ? bn_hash_float_to_float2(v3) : bn_hash_float_to_float2(v4);
return bits_of(comp == 0 ? r.x : r.y);
}
case OP_HFF3: {
vec3 r = (dim == 1) ? bn_hash_float_to_float3(x)
: (dim == 2) ? bn_hash_float_to_float3(v2)
: (dim == 3) ? bn_hash_float_to_float3(v3) : bn_hash_float_to_float3(v4);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : r.z));
}
case OP_HFF4: {
vec4 r = bn_hash_float_to_float4(v4);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : (comp == 2 ? r.z : r.w)));
}
case OP_WRAP_MUL: return (uint32_t)bn_wrap_mul(ia, ib);
case OP_WRAP_ADD: return (uint32_t)bn_wrap_add(ia, ib);
case OP_INT_TO_FLOAT_01: return bits_of(bn_int_to_float_01(ia));
case OP_PCG2D: { ivec2 r = bn_hash_pcg2d_i(i2); return (uint32_t)(comp == 0 ? r.x : r.y); }
case OP_PCG3D: { ivec3 r = bn_hash_pcg3d_i(i3);
return (uint32_t)(comp == 0 ? r.x : (comp == 1 ? r.y : r.z)); }
case OP_PCG4D: { ivec4 r = bn_hash_pcg4d_i(i4);
return (uint32_t)(comp == 0 ? r.x : (comp == 1 ? r.y : (comp == 2 ? r.z : r.w))); }
case OP_INT_TO_FLOAT2: { vec2 r = bn_hash_int_to_float2(i2);
return bits_of(comp == 0 ? r.x : r.y); }
case OP_INT_TO_FLOAT3_I3: { vec3 r = bn_hash_int_to_float3(i3);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : r.z)); }
case OP_INT_TO_FLOAT4: { vec4 r = bn_hash_int_to_float4(i4);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : (comp == 2 ? r.z : r.w))); }
case OP_INT_TO_FLOAT3_I2: { vec3 r = bn_hash_int_to_float3(i2);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : r.z)); }
case OP_INT_TO_FLOAT3_I4: { vec3 r = bn_hash_int_to_float3(i4);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : r.z)); }
case OP_PERLIN_NOISE:
if (dim == 1) return bits_of(bp_perlin_noise(x));
if (dim == 2) return bits_of(bp_perlin_noise(v2));
if (dim == 3) return bits_of(bp_perlin_noise(v3));
return bits_of(bp_perlin_noise(v4));
case OP_PERLIN_SIGNED:
if (dim == 1) return bits_of(bp_perlin_signed(x));
if (dim == 2) return bits_of(bp_perlin_signed(v2));
if (dim == 3) return bits_of(bp_perlin_signed(v3));
return bits_of(bp_perlin_signed(v4));
case OP_PERLIN_SELECT:
if (dim == 1) return bits_of(bp_perlin_select_1(x, detail, rough, lac, off, gain, type, nrm));
if (dim == 2) return bits_of(bp_perlin_select_2(v2, detail, rough, lac, off, gain, type, nrm));
if (dim == 3) return bits_of(bp_perlin_select_3(v3, detail, rough, lac, off, gain, type, nrm));
return bits_of(bp_perlin_select_4(v4, detail, rough, lac, off, gain, type, nrm));
case OP_PERLIN_DISTORTED:
if (dim == 1) return bits_of(bp_perlin_fractal_distorted_1(x, detail, rough, lac, off, gain, dist, type, nrm));
if (dim == 2) return bits_of(bp_perlin_fractal_distorted_2(v2, detail, rough, lac, off, gain, dist, type, nrm));
if (dim == 3) return bits_of(bp_perlin_fractal_distorted_3(v3, detail, rough, lac, off, gain, dist, type, nrm));
return bits_of(bp_perlin_fractal_distorted_4(v4, detail, rough, lac, off, gain, dist, type, nrm));
case OP_PERLIN_FLOAT3: {
vec3 r = (dim == 1) ? bp_perlin_float3_fractal_distorted(x, detail, rough, lac, off, gain, dist, type, nrm)
: (dim == 2) ? bp_perlin_float3_fractal_distorted(v2, detail, rough, lac, off, gain, dist, type, nrm)
: (dim == 3) ? bp_perlin_float3_fractal_distorted(v3, detail, rough, lac, off, gain, dist, type, nrm)
: bp_perlin_float3_fractal_distorted(v4, detail, rough, lac, off, gain, dist, type, nrm);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : r.z));
}
case OP_VORONOI_XFX: {
BvParams p = params_from(w);
BvOutput o = (dim == 1) ? bv_fractal_voronoi_x_fx_1(p, x, true)
: (dim == 2) ? bv_fractal_voronoi_x_fx_2(p, v2, true)
: (dim == 3) ? bv_fractal_voronoi_x_fx_3(p, v3, true)
: bv_fractal_voronoi_x_fx_4(p, v4, true);
if (comp == 0) return bits_of(o.distance);
if (comp <= 3) return bits_of(comp == 1 ? o.color.x : (comp == 2 ? o.color.y : o.color.z));
return bits_of(comp == 4 ? o.position.x : (comp == 5 ? o.position.y
: (comp == 6 ? o.position.z : o.position.w)));
}
case OP_VORONOI_EDGE: {
BvParams p = params_from(w);
return bits_of((dim == 1) ? bv_fractal_voronoi_distance_to_edge_1(p, x)
: (dim == 2) ? bv_fractal_voronoi_distance_to_edge_2(p, v2)
: (dim == 3) ? bv_fractal_voronoi_distance_to_edge_3(p, v3)
: bv_fractal_voronoi_distance_to_edge_4(p, v4));
}
case OP_SIMPLEX: return bits_of(ps_simplex_noise_3d(x, y, z));
case OP_SIMPLEX_FBM:
return bits_of(ps_fbm_noise_3d(x, y, z, type, rough, lac));
case OP_RAMP_BLEND: {
// A in a..d, B in the fractal slots, factor in W_DIST -- the record
// is reused rather than widened again.
vec4 A(x, y, z, wv);
vec4 Bc(detail, rough, lac, off);
vec4 r = br_ramp_blend(type, A, dist, Bc);
return bits_of(comp == 0 ? r.x : (comp == 1 ? r.y : (comp == 2 ? r.z : r.w)));
}
case OP_VORONOI_RADIUS: {
BvParams p = params_from(w);
return bits_of((dim == 1) ? bv_voronoi_n_sphere_radius(p, x)
: (dim == 2) ? bv_voronoi_n_sphere_radius(p, v2)
: (dim == 3) ? bv_voronoi_n_sphere_radius(p, v3)
: bv_voronoi_n_sphere_radius(p, v4));
}
// --- erosion_filter.inc, through dcerosion ---------------------------
case OP_ERO_HASH: {
dcerosion::vec2 q = dcerosion::erosion_hash(dcerosion::vec2(ex, ey));
return bits_of(comp == 0 ? q.x : q.y);
}
case OP_ERO_GRAD_NOISED: {
dcerosion::vec3 q = dcerosion::grad_noised_2d(dcerosion::vec2(ex, ey));
return bits_of(comp == 0 ? q.x : (comp == 1 ? q.y : q.z));
}
case OP_ERO_FRACTAL: {
dcerosion::vec3 q = dcerosion::fractal_noise_2d(dcerosion::vec2(ex, ey),
es(w, 0), (int)w[W_E_OCT], es(w, 1), es(w, 2));
return bits_of(comp == 0 ? q.x : (comp == 1 ? q.y : q.z));
}
case OP_ERO_CLAMP01: return bits_of(dcerosion::clamp01(es(w, 0)));
case OP_ERO_POW_INV: return bits_of(dcerosion::pow_inv(es(w, 0), es(w, 1)));
case OP_ERO_EASE_OUT: return bits_of(dcerosion::ease_out(es(w, 0)));
case OP_ERO_SMOOTH_START:
return bits_of(dcerosion::smooth_start(es(w, 0), es(w, 1)));
case OP_ERO_NORMALIZE_SAFE: {
dcerosion::vec2 q = dcerosion::normalize_safe(dcerosion::vec2(ex, ey));
return bits_of(comp == 0 ? q.x : q.y);
}
case OP_ERO_PHACELLE: {
dcerosion::vec4 q = dcerosion::phacelle_noise(dcerosion::vec2(ex, ey),
dcerosion::vec2(enx, eny), es(w, 0), es(w, 1), es(w, 2));
return bits_of(comp == 0 ? q.x : (comp == 1 ? q.y : (comp == 2 ? q.z : q.w)));
}
case OP_ERO_FILTER: {
dcerosion::ErosionResult r = dcerosion::erosion_filter(
dcerosion::vec2(ex, ey), dcerosion::vec3(enx, eny, enz),
es(w, 0), es(w, 1), es(w, 2), es(w, 3),
e4(w, W_E_R0), e4(w, W_E_O0), e2(w, W_E_A0),
es(w, 4), (int)w[W_E_OCT], es(w, 5),
es(w, 6), es(w, 7), es(w, 8), es(w, 9));
return bits_of(comp == 0 ? r.delta.x : (comp == 1 ? r.delta.y
: (comp == 2 ? r.delta.z : (comp == 3 ? r.magnitude : r.ridge))));
}
case OP_ERO_EVAL_FULL: {
dcerosion::vec2 q = dcerosion::erosion_node_eval_full(
dcerosion::vec2(ex, ey), dcerosion::vec3(enx, eny, enz),
es(w, 0), es(w, 1), es(w, 2), (int)w[W_E_OCT], es(w, 3),
es(w, 4), es(w, 5), es(w, 6), es(w, 7), es(w, 8),
es(w, 9), es(w, 10), e2(w, W_E_S0 + 11),
e4(w, W_E_R0), e4(w, W_E_O0), e2(w, W_E_A0), es(w, 13));
return bits_of(comp == 0 ? q.x : q.y);
}
case OP_ERO_EVAL_GRAD: {
dcerosion::vec4 q = dcerosion::erosion_node_eval_grad(
dcerosion::vec2(ex, ey), dcerosion::vec3(enx, eny, enz),
es(w, 0), es(w, 1), es(w, 2), (int)w[W_E_OCT], es(w, 3),
es(w, 4), es(w, 5), es(w, 6), es(w, 7), es(w, 8),
es(w, 9), es(w, 10), e2(w, W_E_S0 + 11),
e4(w, W_E_R0), e4(w, W_E_O0), e2(w, W_E_A0), es(w, 13));
return bits_of(comp == 0 ? q.x : (comp == 1 ? q.y : (comp == 2 ? q.z : q.w)));
}
default:
fprintf(stderr, "FATAL: cpu dispatcher has no arm for op %u\n", op);
exit(1);
}
}
// --- the GLSL half. O(families), not O(cases). ------------------------------
static const char *DISPATCH_GLSL = R"DISPATCH(
layout(local_size_x = 64) in;
layout(std430, binding = 0) buffer Out { uint v[]; } out_v;
layout(std430, binding = 1) readonly buffer In { uint v[]; } in_v;
uint gW(uint b, uint k) { return in_v.v[b + k]; }
float gF(uint b, uint k) { return uintBitsToFloat(in_v.v[b + k]); }
// The erosion block's scalar k, so this side reads it the way es(w, k) does.
float eS(uint b, uint k) { return gF(b, W_E_S0 + k); }
void main() {
uint idx = gl_GlobalInvocationID.x;
if (idx >= CASE_COUNT) { return; }
uint b = idx * STRIDE_W;
uint op = gW(b, 0u);
uint dim = gW(b, 1u);
int type = int(gW(b, 2u));
bool nrm = (gW(b, 3u) & 1u) != 0u;
uint comp = gW(b, 4u);
uint ua = gW(b, 5u), ub = gW(b, 6u), uc = gW(b, 7u), ud = gW(b, 8u);
int ia = int(ua), ib = int(ub), ic = int(uc), id = int(ud);
float x = gF(b, 5u), y = gF(b, 6u), z = gF(b, 7u), w = gF(b, 8u);
float detail = gF(b, 9u), rough = gF(b, 10u), lac = gF(b, 11u);
float off = gF(b, 12u), gain = gF(b, 13u), dist = gF(b, 14u);
vec2 v2 = vec2(x, y); vec3 v3 = vec3(x, y, z); vec4 v4 = vec4(x, y, z, w);
ivec2 i2 = ivec2(ia, ib); ivec3 i3 = ivec3(ia, ib, ic); ivec4 i4 = ivec4(ia, ib, ic, id);
// The erosion block. Slot names come from the C++ enum via SLOT_NAMES.
vec2 ep = vec2(gF(b, W_E_PX), gF(b, W_E_PY));
vec3 en = vec3(gF(b, W_E_NX), gF(b, W_E_NY), gF(b, W_E_NZ));
vec4 e_round = vec4(gF(b, W_E_R0), gF(b, W_E_R0 + 1u), gF(b, W_E_R0 + 2u), gF(b, W_E_R0 + 3u));
vec4 e_onset = vec4(gF(b, W_E_O0), gF(b, W_E_O0 + 1u), gF(b, W_E_O0 + 2u), gF(b, W_E_O0 + 3u));
vec2 e_slope = vec2(gF(b, W_E_A0), gF(b, W_E_A0 + 1u));
int e_oct = int(gW(b, W_E_OCT));
uint r = 0u;
if (op == OP_HASH1) { r = bn_hash(ua); }
else if (op == OP_HASH2) { r = bn_hash(ua, ub); }
else if (op == OP_HASH3) { r = bn_hash(ua, ub, uc); }
else if (op == OP_HASH4) { r = bn_hash(ua, ub, uc, ud); }
else if (op == OP_UINT_TO_FLOAT_01) { r = floatBitsToUint(bn_uint_to_float_01(ua)); }
else if (op == OP_HASH_TO_FLOAT1) { r = floatBitsToUint(bn_hash_to_float(ua)); }
else if (op == OP_HASH_TO_FLOAT2) { r = floatBitsToUint(bn_hash_to_float(ua, ub)); }
else if (op == OP_HASH_TO_FLOAT3) { r = floatBitsToUint(bn_hash_to_float(ua, ub, uc)); }
else if (op == OP_HASH_TO_FLOAT4) { r = floatBitsToUint(bn_hash_to_float(ua, ub, uc, ud)); }
else if (op == OP_HASH_FLOAT) {
if (dim == 1u) { r = bn_hash_float(x); }
else if (dim == 2u) { r = bn_hash_float(v2); }
else if (dim == 3u) { r = bn_hash_float(v3); }
else { r = bn_hash_float(v4); }
}
else if (op == OP_HFF) {
float f;
if (dim == 1u) { f = bn_hash_float_to_float(x); }
else if (dim == 2u) { f = bn_hash_float_to_float(v2); }
else if (dim == 3u) { f = bn_hash_float_to_float(v3); }
else { f = bn_hash_float_to_float(v4); }
r = floatBitsToUint(f);
}
else if (op == OP_HFF2) {
vec2 q;
if (dim == 2u) { q = bn_hash_float_to_float2(v2); }
else if (dim == 3u) { q = bn_hash_float_to_float2(v3); }
else { q = bn_hash_float_to_float2(v4); }
r = floatBitsToUint(comp == 0u ? q.x : q.y);
}
else if (op == OP_HFF3) {
vec3 q;
if (dim == 1u) { q = bn_hash_float_to_float3(x); }
else if (dim == 2u) { q = bn_hash_float_to_float3(v2); }
else if (dim == 3u) { q = bn_hash_float_to_float3(v3); }
else { q = bn_hash_float_to_float3(v4); }
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z));
}
else if (op == OP_HFF4) {
vec4 q = bn_hash_float_to_float4(v4);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : (comp == 2u ? q.z : q.w)));
}
else if (op == OP_WRAP_MUL) { r = uint(bn_wrap_mul(ia, ib)); }
else if (op == OP_WRAP_ADD) { r = uint(bn_wrap_add(ia, ib)); }
else if (op == OP_INT_TO_FLOAT_01) { r = floatBitsToUint(bn_int_to_float_01(ia)); }
else if (op == OP_PCG2D) { ivec2 q = bn_hash_pcg2d_i(i2); r = uint(comp == 0u ? q.x : q.y); }
else if (op == OP_PCG3D) { ivec3 q = bn_hash_pcg3d_i(i3);
r = uint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z)); }
else if (op == OP_PCG4D) { ivec4 q = bn_hash_pcg4d_i(i4);
r = uint(comp == 0u ? q.x : (comp == 1u ? q.y : (comp == 2u ? q.z : q.w))); }
else if (op == OP_INT_TO_FLOAT2) { vec2 q = bn_hash_int_to_float2(i2);
r = floatBitsToUint(comp == 0u ? q.x : q.y); }
else if (op == OP_INT_TO_FLOAT3_I3) { vec3 q = bn_hash_int_to_float3(i3);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z)); }
else if (op == OP_INT_TO_FLOAT4) { vec4 q = bn_hash_int_to_float4(i4);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : (comp == 2u ? q.z : q.w))); }
else if (op == OP_INT_TO_FLOAT3_I2) { vec3 q = bn_hash_int_to_float3(i2);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z)); }
else if (op == OP_INT_TO_FLOAT3_I4) { vec3 q = bn_hash_int_to_float3(i4);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z)); }
else if (op == OP_PERLIN_NOISE) {
float f;
if (dim == 1u) { f = bp_perlin_noise(x); }
else if (dim == 2u) { f = bp_perlin_noise(v2); }
else if (dim == 3u) { f = bp_perlin_noise(v3); }
else { f = bp_perlin_noise(v4); }
r = floatBitsToUint(f);
}
else if (op == OP_PERLIN_SIGNED) {
float f;
if (dim == 1u) { f = bp_perlin_signed(x); }
else if (dim == 2u) { f = bp_perlin_signed(v2); }
else if (dim == 3u) { f = bp_perlin_signed(v3); }
else { f = bp_perlin_signed(v4); }
r = floatBitsToUint(f);
}
else if (op == OP_PERLIN_SELECT) {
float f;
if (dim == 1u) { f = bp_perlin_select_1(x, detail, rough, lac, off, gain, type, nrm); }
else if (dim == 2u) { f = bp_perlin_select_2(v2, detail, rough, lac, off, gain, type, nrm); }
else if (dim == 3u) { f = bp_perlin_select_3(v3, detail, rough, lac, off, gain, type, nrm); }
else { f = bp_perlin_select_4(v4, detail, rough, lac, off, gain, type, nrm); }
r = floatBitsToUint(f);
}
else if (op == OP_PERLIN_DISTORTED) {
float f;
if (dim == 1u) { f = bp_perlin_fractal_distorted_1(x, detail, rough, lac, off, gain, dist, type, nrm); }
else if (dim == 2u) { f = bp_perlin_fractal_distorted_2(v2, detail, rough, lac, off, gain, dist, type, nrm); }
else if (dim == 3u) { f = bp_perlin_fractal_distorted_3(v3, detail, rough, lac, off, gain, dist, type, nrm); }
else { f = bp_perlin_fractal_distorted_4(v4, detail, rough, lac, off, gain, dist, type, nrm); }
r = floatBitsToUint(f);
}
else if (op == OP_PERLIN_FLOAT3) {
vec3 q;
if (dim == 1u) { q = bp_perlin_float3_fractal_distorted(x, detail, rough, lac, off, gain, dist, type, nrm); }
else if (dim == 2u) { q = bp_perlin_float3_fractal_distorted(v2, detail, rough, lac, off, gain, dist, type, nrm); }
else if (dim == 3u) { q = bp_perlin_float3_fractal_distorted(v3, detail, rough, lac, off, gain, dist, type, nrm); }
else { q = bp_perlin_float3_fractal_distorted(v4, detail, rough, lac, off, gain, dist, type, nrm); }
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z));
}
else if (op == OP_SIMPLEX) {
r = floatBitsToUint(ps_simplex_noise_3d(x, y, z));
}
else if (op == OP_SIMPLEX_FBM) {
r = floatBitsToUint(ps_fbm_noise_3d(x, y, z, type, rough, lac));
}
else if (op == OP_RAMP_BLEND) {
vec4 A = vec4(x, y, z, w);
vec4 Bc = vec4(detail, rough, lac, off);
vec4 q = br_ramp_blend(type, A, dist, Bc);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : (comp == 2u ? q.z : q.w)));
}
else if (op == OP_VORONOI_XFX || op == OP_VORONOI_EDGE || op == OP_VORONOI_RADIUS) {
BvParams p;
p.scale = gF(b, 15u);
p.detail = detail;
p.roughness = rough;
p.lacunarity = lac;
p.smoothness = gF(b, 16u);
p.exponent = gF(b, 17u);
p.randomness = gF(b, 18u);
p.max_distance = gF(b, 19u);
p.normalize = nrm;
p.feature = type;
p.metric = int(gW(b, 20u));
if (op == OP_VORONOI_RADIUS) {
float f;
if (dim == 1u) { f = bv_voronoi_n_sphere_radius(p, x); }
else if (dim == 2u) { f = bv_voronoi_n_sphere_radius(p, v2); }
else if (dim == 3u) { f = bv_voronoi_n_sphere_radius(p, v3); }
else { f = bv_voronoi_n_sphere_radius(p, v4); }
r = floatBitsToUint(f);
} else if (op == OP_VORONOI_EDGE) {
float f;
if (dim == 1u) { f = bv_fractal_voronoi_distance_to_edge_1(p, x); }
else if (dim == 2u) { f = bv_fractal_voronoi_distance_to_edge_2(p, v2); }
else if (dim == 3u) { f = bv_fractal_voronoi_distance_to_edge_3(p, v3); }
else { f = bv_fractal_voronoi_distance_to_edge_4(p, v4); }
r = floatBitsToUint(f);
} else {
BvOutput o;
if (dim == 1u) { o = bv_fractal_voronoi_x_fx_1(p, x, true); }
else if (dim == 2u) { o = bv_fractal_voronoi_x_fx_2(p, v2, true); }
else if (dim == 3u) { o = bv_fractal_voronoi_x_fx_3(p, v3, true); }
else { o = bv_fractal_voronoi_x_fx_4(p, v4, true); }
float f;
if (comp == 0u) { f = o.distance; }
else if (comp == 1u) { f = o.color.x; }
else if (comp == 2u) { f = o.color.y; }
else if (comp == 3u) { f = o.color.z; }
else if (comp == 4u) { f = o.position.x; }
else if (comp == 5u) { f = o.position.y; }
else if (comp == 6u) { f = o.position.z; }
else { f = o.position.w; }
r = floatBitsToUint(f);
}
}
else if (op == OP_ERO_HASH) {
vec2 q = erosion_hash(ep);
r = floatBitsToUint(comp == 0u ? q.x : q.y);
}
else if (op == OP_ERO_GRAD_NOISED) {
vec3 q = grad_noised_2d(ep);
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z));
}
else if (op == OP_ERO_FRACTAL) {
vec3 q = fractal_noise_2d(ep, eS(b, 0u), e_oct, eS(b, 1u), eS(b, 2u));
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : q.z));
}
else if (op == OP_ERO_CLAMP01) { r = floatBitsToUint(clamp01(eS(b, 0u))); }
else if (op == OP_ERO_POW_INV) { r = floatBitsToUint(pow_inv(eS(b, 0u), eS(b, 1u))); }
else if (op == OP_ERO_EASE_OUT) { r = floatBitsToUint(ease_out(eS(b, 0u))); }
else if (op == OP_ERO_SMOOTH_START) {
r = floatBitsToUint(smooth_start(eS(b, 0u), eS(b, 1u)));
}
else if (op == OP_ERO_NORMALIZE_SAFE) {
vec2 q = normalize_safe(ep);
r = floatBitsToUint(comp == 0u ? q.x : q.y);
}
else if (op == OP_ERO_PHACELLE) {
vec4 q = phacelle_noise(ep, vec2(en.x, en.y), eS(b, 0u), eS(b, 1u), eS(b, 2u));
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : (comp == 2u ? q.z : q.w)));
}
else if (op == OP_ERO_FILTER) {
ErosionResult q = erosion_filter(ep, en,
eS(b, 0u), eS(b, 1u), eS(b, 2u), eS(b, 3u),
e_round, e_onset, e_slope,
eS(b, 4u), e_oct, eS(b, 5u),
eS(b, 6u), eS(b, 7u), eS(b, 8u), eS(b, 9u));
float f;
if (comp == 0u) { f = q.delta.x; }
else if (comp == 1u) { f = q.delta.y; }
else if (comp == 2u) { f = q.delta.z; }
else if (comp == 3u) { f = q.magnitude; }
else { f = q.ridge; }
r = floatBitsToUint(f);
}
else if (op == OP_ERO_EVAL_FULL) {
vec2 q = erosion_node_eval_full(ep, en,
eS(b, 0u), eS(b, 1u), eS(b, 2u), e_oct, eS(b, 3u),
eS(b, 4u), eS(b, 5u), eS(b, 6u), eS(b, 7u), eS(b, 8u),
eS(b, 9u), eS(b, 10u), vec2(eS(b, 11u), eS(b, 12u)),
e_round, e_onset, e_slope, eS(b, 13u));
r = floatBitsToUint(comp == 0u ? q.x : q.y);
}
else if (op == OP_ERO_EVAL_GRAD) {
vec4 q = erosion_node_eval_grad(ep, en,
eS(b, 0u), eS(b, 1u), eS(b, 2u), e_oct, eS(b, 3u),
eS(b, 4u), eS(b, 5u), eS(b, 6u), eS(b, 7u), eS(b, 8u),
eS(b, 9u), eS(b, 10u), vec2(eS(b, 11u), eS(b, 12u)),
e_round, e_onset, e_slope, eS(b, 13u));
r = floatBitsToUint(comp == 0u ? q.x : (comp == 1u ? q.y : (comp == 2u ? q.z : q.w)));
}
out_v.v[idx] = r;
}
)DISPATCH";
// ============================================================================
// THE CASES
//
// Inputs are chosen where these functions actually break: zero, the all-ones
// word a shift can eat, single-bit words that expose a rotate written as a
// shift, the two's-complement spelling of a negative lattice coordinate (an
// ORDINARY input -- Blender indexes its lattice with signed ints cast to
// uint), real float bit patterns, and for perlin a row past the 100000 wrap
// where the truncated modulo and the half-unit nudge past a million do
// something.
// ============================================================================
struct B {
uint32_t w[STRIDE];
explicit B(uint32_t op) {
memset(w, 0, sizeof(w));
w[W_OP] = op;
w[W_DIM] = 3;
}
B &dim(uint32_t d) { w[W_DIM] = d; return *this; }
B &type(int t) { w[W_TYPE] = (uint32_t)t; return *this; }
B &norm(bool n) { w[W_FLAGS] = n ? 1u : 0u; return *this; }
B &comp(uint32_t c) { w[W_COMP] = c; return *this; }
B &u(uint32_t a, uint32_t b = 0, uint32_t c = 0, uint32_t d = 0) {
w[W_A] = a; w[W_B] = b; w[W_C] = c; w[W_D] = d; return *this;
}
B &f(float a, float b = 0, float c = 0, float d = 0) {
w[W_A] = bits_of(a); w[W_B] = bits_of(b);
w[W_C] = bits_of(c); w[W_D] = bits_of(d); return *this;
}
B &vor(float scale, float smooth, float expo, float rand01, float maxd, int metric) {
w[W_SCALE] = bits_of(scale); w[W_SMOOTH] = bits_of(smooth);
w[W_EXPONENT] = bits_of(expo); w[W_RANDOMNESS] = bits_of(rand01);
w[W_MAXDIST] = bits_of(maxd); w[W_METRIC] = (uint32_t)metric; return *this;
}
// --- the erosion block ---------------------------------------------------
B &ep(float px, float py) {
w[W_E_PX] = bits_of(px); w[W_E_PY] = bits_of(py); return *this;
}
B &en(float nx, float ny, float nz) {
w[W_E_NX] = bits_of(nx); w[W_E_NY] = bits_of(ny); w[W_E_NZ] = bits_of(nz);
return *this;
}
B &es(int k, float v) { w[W_E_S0 + k] = bits_of(v); return *this; }
B &eoct(int o) { w[W_E_OCT] = (uint32_t)o; return *this; }
B &eround(float a_, float b_, float c_, float d_) {
w[W_E_R0] = bits_of(a_); w[W_E_R0 + 1] = bits_of(b_);
w[W_E_R0 + 2] = bits_of(c_); w[W_E_R0 + 3] = bits_of(d_); return *this;
}
B &eonset(float a_, float b_, float c_, float d_) {
w[W_E_O0] = bits_of(a_); w[W_E_O0 + 1] = bits_of(b_);
w[W_E_O0 + 2] = bits_of(c_); w[W_E_O0 + 3] = bits_of(d_); return *this;
}
B &eslope(float a_, float b_) {
w[W_E_A0] = bits_of(a_); w[W_E_A0 + 1] = bits_of(b_); return *this;
}
B &frac(float detail, float rough, float lac, float off, float gain, float dist) {
w[W_DETAIL] = bits_of(detail); w[W_ROUGH] = bits_of(rough);
w[W_LAC] = bits_of(lac); w[W_OFF] = bits_of(off);
w[W_GAIN] = bits_of(gain); w[W_DIST] = bits_of(dist); return *this;
}
};
static std::vector<Rec> g_cases;
static void add(const std::string &name, bool exact, const B &b, double tol_abs = 0.0) {
Rec r;
memcpy(r.w, b.w, sizeof(r.w));
r.name = name;
r.exact = exact;
r.tol_abs = tol_abs;
r.cpu = eval_cpu(r.w);
g_cases.push_back(r);
}
static void build_cases() {
const uint32_t KU[] = {
0u, 1u, 2u, 7u, 16u, 255u,
0x0000FFFFu, 0xFFFF0000u, 0x80000000u, 0xFFFFFFFFu,
0xdeadbeefu, 0xcafebabeu, 0x12345678u,
(uint32_t)-1, (uint32_t)-7, (uint32_t)-1000000
};
const int NU = (int)(sizeof(KU) / sizeof(KU[0]));
const float KF[] = {
0.0f, -0.0f, 1.0f, -1.0f, 0.5f, -0.25f, 2.5f, 7.0f,
3.14159265f, -7.25f, 1e7f, -1e7f, 1e-7f, 123456.789f, -0.001f, 65536.0f
};
const int NF = (int)(sizeof(KF) / sizeof(KF[0]));
// --- EXACT: the hash layer.
for (int i = 0; i < NU; i++) {
const uint32_t a = KU[i], b = KU[(i + 3) % NU], c = KU[(i + 7) % NU], d = KU[(i + 11) % NU];
add("bn_hash1", true, B(OP_HASH1).u(a));
add("bn_hash2", true, B(OP_HASH2).u(a, b));
add("bn_hash3", true, B(OP_HASH3).u(a, b, c));
add("bn_hash4", true, B(OP_HASH4).u(a, b, c, d));
add("bn_uint_to_float_01", true, B(OP_UINT_TO_FLOAT_01).u(a));
add("bn_hash_to_float1", true, B(OP_HASH_TO_FLOAT1).u(a));
add("bn_hash_to_float2", true, B(OP_HASH_TO_FLOAT2).u(a, b));
add("bn_hash_to_float3", true, B(OP_HASH_TO_FLOAT3).u(a, b, c));
add("bn_hash_to_float4", true, B(OP_HASH_TO_FLOAT4).u(a, b, c, d));
// PCG. Signed throughout, so the shift is ARITHMETIC and the wrap goes
// through uint on both sides -- neither language defines signed overflow.
add("bn_wrap_mul", true, B(OP_WRAP_MUL).u(a, b));
add("bn_wrap_add", true, B(OP_WRAP_ADD).u(a, b));
add("bn_int_to_float_01", true, B(OP_INT_TO_FLOAT_01).u(a));
for (uint32_t k = 0; k < 2; k++) {
add("bn_hash_pcg2d_i", true, B(OP_PCG2D).u(a, b).comp(k));
add("bn_hash_int_to_float2", true, B(OP_INT_TO_FLOAT2).u(a, b).comp(k));
}
for (uint32_t k = 0; k < 3; k++) {
add("bn_hash_pcg3d_i", true, B(OP_PCG3D).u(a, b, c).comp(k));
add("bn_hash_int_to_float3(i3)", true, B(OP_INT_TO_FLOAT3_I3).u(a, b, c).comp(k));
add("bn_hash_int_to_float3(i2)", true, B(OP_INT_TO_FLOAT3_I2).u(a, b).comp(k));
add("bn_hash_int_to_float3(i4)", true, B(OP_INT_TO_FLOAT3_I4).u(a, b, c, d).comp(k));
}
for (uint32_t k = 0; k < 4; k++) {
add("bn_hash_pcg4d_i", true, B(OP_PCG4D).u(a, b, c, d).comp(k));
add("bn_hash_int_to_float4", true, B(OP_INT_TO_FLOAT4).u(a, b, c, d).comp(k));
}
}
// The float-keyed families. The bit cast is the interesting part: C++ needs
// memcpy where GLSL has floatBitsToUint.
for (int i = 0; i < NF; i++) {
const float x = KF[i], y = KF[(i + 3) % NF], z = KF[(i + 5) % NF], w = KF[(i + 9) % NF];
for (uint32_t d = 1; d <= 4; d++) {
add("bn_hash_float", true, B(OP_HASH_FLOAT).dim(d).f(x, y, z, w));
add("bn_hff", true, B(OP_HFF).dim(d).f(x, y, z, w));
for (uint32_t k = 0; k < 3; k++) {
add("bn_hff3", true, B(OP_HFF3).dim(d).comp(k).f(x, y, z, w));
}
if (d >= 2) {
for (uint32_t k = 0; k < 2; k++) {
add("bn_hff2", true, B(OP_HFF2).dim(d).comp(k).f(x, y, z, w));
}
}
}
for (uint32_t k = 0; k < 4; k++) {
add("bn_hff4", true, B(OP_HFF4).comp(k).f(x, y, z, w));
}
}
// --- NEAR: perlin.
const float PC[][4] = {
{ 0.0f, 0.0f, 0.0f, 0.0f },
{ 0.5f, 1.5f, 2.5f, 3.5f },
{ -7.25f, 0.5f, 13.75f, -2.5f },
{ 123.456f, -78.9f, 0.001f, 44.0f },
// past the 100000 wrap and past the million where the nudge kicks in
{ 123456.75f, -234567.5f, 1000001.5f, -1000000.5f },
};
const int NP = (int)(sizeof(PC) / sizeof(PC[0]));
// detail deliberately FRACTIONAL in most rows: rmd != 0 takes the
// extra-octave branch that four of the five fractal types have.
const float FP[][5] = {
{ 3.0f, 0.5f, 2.0f, 0.0f, 1.0f },
{ 2.5f, 0.5f, 2.0f, 0.7f, 1.0f },
{ 5.25f, 0.75f, 2.5f, 0.3f, 1.5f },
{ 0.0f, 0.5f, 2.0f, 1.0f, 2.0f },
{ 7.5f, 0.25f, 1.75f, 0.5f, 0.5f },
};
const int NFP = (int)(sizeof(FP) / sizeof(FP[0]));
for (int i = 0; i < NP; i++) {
const float x = PC[i][0], y = PC[i][1], z = PC[i][2], w = PC[i][3];
for (uint32_t d = 1; d <= 4; d++) {
add("bp_perlin_noise", false, B(OP_PERLIN_NOISE).dim(d).f(x, y, z, w));
add("bp_perlin_signed", false, B(OP_PERLIN_SIGNED).dim(d).f(x, y, z, w));
}
for (int f = 0; f < NFP; f++) {
const float dt = FP[f][0], rg = FP[f][1], lc = FP[f][2];
const float of = FP[f][3], gn = FP[f][4];
// Every fractal type and every dimension: twenty functions stamped
// from one macro body, and a wrong suffix would make one dimension
// quietly evaluate another. The per-dimension scale constants --
// 0.2500, 0.6616, 0.9820, 0.8344 -- make that show up as a scale
// error rather than as noise.
for (int t = 0; t < 5; t++) {
for (int n = 0; n < 2; n++) {
const bool nb = (n != 0);
for (uint32_t d = 1; d <= 4; d++) {
char nm[64];
snprintf(nm, sizeof(nm), "bp_perlin_select_%u.t%d", d, t);
add(nm, false, B(OP_PERLIN_SELECT).dim(d).type(t).norm(nb)
.f(x, y, z, w).frac(dt, rg, lc, of, gn, 0.0f));
}
}
}
for (uint32_t d = 1; d <= 4; d++) {
char nm[64];
snprintf(nm, sizeof(nm), "bp_perlin_fractal_distorted_%u", d);
add(nm, false, B(OP_PERLIN_DISTORTED).dim(d).type(1).norm(true)
.f(x, y, z, w).frac(dt, rg, lc, of, gn, 0.75f));
for (uint32_t k = 0; k < 3; k++) {
snprintf(nm, sizeof(nm), "bp_perlin_float3(d%u)", d);
add(nm, false, B(OP_PERLIN_FLOAT3).dim(d).type(1).norm(true).comp(k)
.f(x, y, z, w).frac(dt, rg, lc, of, gn, 0.75f));
}
}
}
}
}
static void build_voronoi_cases() {
// Every dimension x feature x metric, with normalize both ways. They are
// not variations on one function: the feature picks which of three
// entry points runs, the metric changes the distance inside the cell
// search, and the dimension changes both the lattice and which of the
// four stamped copies is called.
const float VC[][4] = {
{ 0.0f, 0.0f, 0.0f, 0.0f },
{ 1.25f, -2.5f, 0.75f, 3.5f },
{ -13.75f, 7.5f, -0.25f, -6.5f },
};
const int NVC = (int)(sizeof(VC) / sizeof(VC[0]));
for (int i = 0; i < NVC; i++) {
const float x = VC[i][0], y = VC[i][1], z = VC[i][2], w = VC[i][3];
const vec2 c2(x, y); const vec3 c3(x, y, z); const vec4 c4(x, y, z, w);
for (uint32_t dim = 1; dim <= 4; dim++) {
for (int feat = 0; feat <= 4; feat++) {
for (int metric = 0; metric <= 3; metric++) {
for (int n = 0; n < 2; n++) {
const bool nb = (n != 0);
const float scale = 1.0f + 0.25f * (float)i;
const float detail = (i == 0) ? 0.0f : 2.5f; // fractional
const float rough = 0.5f;
const float lac = 2.0f;
const float smooth = (feat == 2) ? 0.25f : 0.0f;
const float expo = (metric == 3 && dim != 1) ? 2.5f : 0.0f;
const float rand01 = 0.85f;
const float r = 0.5f + 0.5f * rand01;
const float f2 = (feat == 1) ? 2.0f : 1.0f;
// max_distance: 1.0 for the edge and radius features,
// where the CPU leaves the struct default standing, and
// the per-dimension cell distance for the rest.
float maxd = 1.0f;
if (feat != 3 && feat != 4) {
BvParams tmp;
tmp.scale = scale; tmp.detail = detail; tmp.roughness = rough;
tmp.lacunarity = lac; tmp.smoothness = smooth; tmp.exponent = expo;
tmp.randomness = rand01; tmp.max_distance = 1.0f;
tmp.normalize = nb; tmp.feature = feat; tmp.metric = metric;
if (dim == 1) {
maxd = r * f2;
} else if (dim == 2) {
maxd = bv_voronoi_distance(vec2(0.0f), vec2(r, r), tmp) * f2;
} else if (dim == 3) {
maxd = bv_voronoi_distance(vec3(0.0f), vec3(r, r, r), tmp) * f2;
} else {
maxd = bv_voronoi_distance(vec4(0.0f), vec4(r, r, r, r), tmp) * f2;
}
}
char nm[80];
if (feat == 4) {
snprintf(nm, sizeof(nm), "bv_n_sphere_radius.d%u.m%d", dim, metric);
add(nm, false, B(OP_VORONOI_RADIUS).dim(dim).type(feat).norm(nb)
.f(x, y, z, w).frac(detail, rough, lac, 0.0f, 0.0f, 0.0f)
.vor(scale, smooth, expo, rand01, maxd, metric));
} else if (feat == 3) {
snprintf(nm, sizeof(nm), "bv_edge.d%u.m%d", dim, metric);
add(nm, false, B(OP_VORONOI_EDGE).dim(dim).type(feat).norm(nb)
.f(x, y, z, w).frac(detail, rough, lac, 0.0f, 0.0f, 0.0f)
.vor(scale, smooth, expo, rand01, maxd, metric));
} else {
// distance, all three colour components and all four
// of position -- the W component is a separate output
// on the node and comes from position.w.
for (uint32_t comp = 0; comp < 8; comp++) {
snprintf(nm, sizeof(nm), "bv_x_fx.d%u.f%d.m%d", dim, feat, metric);
add(nm, false, B(OP_VORONOI_XFX).dim(dim).type(feat).norm(nb)
.comp(comp).f(x, y, z, w)
.frac(detail, rough, lac, 0.0f, 0.0f, 0.0f)
.vor(scale, smooth, expo, rand01, maxd, metric));
}
}
}
}
}
}
}
}
static void build_ramp_cases() {
// All nineteen blend modes on all four channels. Several of them are not a
// lerp of anything, four go through HSV, and the alpha channel is left
// alone by every mode except BLEND -- an asymmetry that is Blender's and
// that a "reasonable" implementation would smooth away.
const float CP[][9] = {
// A r,g,b,a B r,g,b,a fac
{ 0.2f, 0.4f, 0.6f, 1.0f, 0.7f, 0.1f, 0.3f, 1.0f, 0.5f },
{ 0.0f, 0.0f, 0.0f, 1.0f, 0.5f, 0.5f, 0.5f, 0.25f, 0.0f }, // black A, fac 0
{ 0.5f, 0.5f, 0.5f, 0.5f, 0.0f, 0.0f, 0.0f, 1.0f, 1.0f }, // grey A, black B, fac 1
{ 0.9f, 0.05f, 0.35f, 0.8f, 0.15f, 0.95f, 0.6f, 0.2f, 0.75f },
};
const int NCP = (int)(sizeof(CP) / sizeof(CP[0]));
for (int i = 0; i < NCP; i++) {
for (int mode = 0; mode <= 18; mode++) {
for (uint32_t comp = 0; comp < 4; comp++) {
char nm[48];
snprintf(nm, sizeof(nm), "br_ramp_blend.m%d", mode);
add(nm, false, B(OP_RAMP_BLEND).type(mode).comp(comp)
.f(CP[i][0], CP[i][1], CP[i][2], CP[i][3])
.frac(CP[i][4], CP[i][5], CP[i][6], CP[i][7], 0.0f, CP[i][8]));
}
}
}
}
static void build_simplex_cases() {
// The Ashima simplex, now that it is a shared source rather than a
// file-local C++ copy paired with a hand-written GLSL one. It is the
// Simplex Noise node's whole arithmetic and Displace's neighbour.
const float SC[][3] = {
{ 0.0f, 0.0f, 0.0f },
{ 0.5f, 1.5f, 2.5f },
{ -7.25f, 13.75f, -0.5f },
{ 123.456f, -78.9f, 0.001f },
};
// ⚠ Coordinates stop at 250, which with eight octaves is an effective
// magnitude of ~3e4. Beyond roughly 1e6 the two machines part company by
// ~1.3e-03: Ashima's permutation (multiply by 34, add 10, mod 289) runs
// out of float precision there and each rounds its own way. That is a
// property of the algorithm in float32, not a divergence between the CPU
// and the GPU texts -- they are the same text -- so it is recorded rather
// than tolerated.
const int NSC = (int)(sizeof(SC) / sizeof(SC[0]));
for (int i = 0; i < NSC; i++) {
const float x = SC[i][0], y = SC[i][1], z = SC[i][2];
add("ps_simplex_noise_3d", false, B(OP_SIMPLEX).f(x, y, z, 0.0f));
// The octave loop is back. It was withdrawn while the CPU side was the
// broken scalar port -- which returned ~±1e5 instead of [-1, 1], so an
// fbm sum was an accumulation of nonsense whose last bits no tolerance
// should have been asked to cover. With the vector form in place the
// function is bounded again and the sum means something.
// ⚠ AN EXPLICIT ABSOLUTE TOLERANCE, and this is the one place in the
// file with one. An fbm sums N octaves of magnitude ~0.9 into a result
// that can be ~0.01 -- the terms largely cancel -- so scaling the
// difference by the RESULT measures the cancellation rather than the
// agreement, and Ashima's permutation loses integer exactness in
// float32 as the effective coordinate grows besides. Measured worst
// over this set: 1.6e-05 absolute at eight octaves. 5e-05 is that with
// headroom, and it is a number taken from a measurement rather than
// one chosen until the gate went green.
for (int oct = 1; oct <= 8; oct++) {
for (int k = 0; k < 2; k++) {
const float rough = k ? 0.75f : 0.5f;
const float lac = k ? 1.75f : 2.0f;
add("ps_fbm_noise_3d", false, B(OP_SIMPLEX_FBM).type(oct)
.f(x, y, z, 0.0f).frac(0.0f, rough, lac, 0.0f, 0.0f, 0.0f), 5e-5);
}
}
}
}
// ============================================================================
// erosion_filter.inc, bottom-up.
//
// The layering is the one the Blender port was built in and for the same
// reason: a divergence in erosion_hash is unreadable once it has been through
// eight octaves, a fade and a gully mask, and each layer here is checked before
// the one that consumes it.
//
// The parameter bundles are the SHIPPED ones -- runevision's reference values,
// read off get_erosion_library_glsl()'s globals and ErosionNoiseNode's socket
// defaults, not invented for the test. A gate on parameters nothing uses proves
// the arithmetic and not the surface.
//
// ---------------------------------------------------------------------------
// ⚠ HOW TO MUTATION-TEST THIS FAMILY, because the obvious way proves nothing.
//
// Every other gate in this repo is mutated by editing the thing under test.
// Here the thing under test is ONE TEXT compiled twice, so editing the .inc
// changes BOTH sides and the gate stays green -- "a transcription used as the
// reference can only ever agree with itself", one level up. The mutation that
// works is to edit the .inc and NOT regenerate the .gen.h: the CPU picks the
// change up through the shim, the GPU keeps the shipped text, and they part.
// That is also exactly the stale-generated-header scenario the generator's own
// docstring promises cannot happen silently.
//
// Measured 2026-09-06, every family caught by its own case:
//
// erosion_hash 2^24 -> 2^24-1 (both lines) 9 red erosion_hash
// quintic fade +10.0 -> +11.0 105 red grad_noised_2d
// fade derivative 30.0 -> 29.0 62 red grad_noised, fractal
// clamp(t, 0.0, 1.0) -> (t, 0.0, 0.9) 521 red clamp01 + downstream
// ease_out 1 - v*v -> 1 - v*v*v 541 red ease_out + downstream
// pow_inv 1 - pow() -> 1 - 0.5*pow() 408 red pow_inv + wrappers
// smooth_start 0.5*s -> 0.4*s 357 red filter, _full, _grad
// normalize_safe 1e-20 -> 1e-1 2 red normalize_safe
// phacelle TAU 6.2831853 -> 6.2831 92 red phacelle_noise
//
// ⛔ AND ONE THAT PROVED NOTHING. TAU 6.28318530717958647 -> 6.28318530717000000
// left the gate GREEN, and it is not a coverage gap: both spellings are the
// same float32 (6.2831855), so the mutation never happened. It is recorded
// because a sweep that had stopped at "9 of 10 caught" would have gone looking
// for a hole in the phacelle cases that is not there.
// ---------------------------------------------------------------------------
// ============================================================================
static void build_erosion_cases() {
// Points chosen where this filter actually breaks: the origin (a lattice
// point, where erosion_hash's floor lands exactly), the negative quadrant
// (GLSL fract is x - floor(x), which for negative x is neither fmod nor
// truncation), a diagonal, and coordinates large enough that the hash's
// 16777216.0 = 2^24 matters.
const float P[][2] = {
{ 0.0f, 0.0f }, { 0.5f, 0.5f }, { 1.0f, 0.0f },
{ 0.37f, 0.81f }, { -0.7f, 1.3f }, { -3.3f, -4.1f },
{ 12.25f, -7.75f }, { 123.5f, 456.25f }, { 1e-4f, -1e-4f },
};
const int NP = (int)(sizeof(P) / sizeof(P[0]));
// EXACT for the hash: it is integer arithmetic to a float conversion, the
// same class as the blender hash layer, and a near-miss is a finding.
for (int i = 0; i < NP; i++) {
for (uint32_t c = 0; c < 2; c++) {
add("erosion_hash", true, B(OP_ERO_HASH).comp(c).ep(P[i][0], P[i][1]));
}
}
// grad_noised_2d is where the f-suffix hazard is sharpest: the quintic fade
// `t*t*t*(t*(t*6.0 - 15.0) + 10.0)` and its derivative `30.0*t*t*(t*(t-2.0)
// + 1.0)` are chained expressions containing bare literals, so on the C++
// side the whole chain evaluates in DOUBLE and narrows once at the end,
// where the shader rounds at every step in float.
for (int i = 0; i < NP; i++) {
for (uint32_t c = 0; c < 3; c++) {
add("grad_noised_2d", false, B(OP_ERO_GRAD_NOISED).comp(c).ep(P[i][0], P[i][1]));
}
}
// The scalar helpers. clamp01 and ease_out are one expression each; pow_inv
// carries a pow(), and smooth_start has the branch whose t >= smoothing side
// is the guard against a division by zero.
const float T[] = { -0.5f, 0.0f, 0.25f, 0.5f, 1.0f, 1.5f, 7.0f };
const int NT = (int)(sizeof(T) / sizeof(T[0]));
for (int i = 0; i < NT; i++) {
add("clamp01", false, B(OP_ERO_CLAMP01).es(0, T[i]));
add("ease_out", false, B(OP_ERO_EASE_OUT).es(0, T[i]));
for (int k = 0; k < 3; k++) {
const float pw[3] = { 0.5f, 1.0f, 2.5f };
add("pow_inv", false, B(OP_ERO_POW_INV).es(0, T[i]).es(1, pw[k]));
const float sm[3] = { 0.0f, 0.1f, 1.0f };
add("smooth_start", false, B(OP_ERO_SMOOTH_START).es(0, T[i]).es(1, sm[k]));
}
}
// normalize_safe's whole reason for existing is the l > 1e-20 branch, so
// the zero vector and a vector just under the threshold are both here.
const float NV[][2] = { { 0.0f, 0.0f }, { 3.0f, 4.0f }, { -1e-21f, 0.0f }, { 1e-8f, -1e-8f } };
for (int i = 0; i < 4; i++) {
for (uint32_t c = 0; c < 2; c++) {
add("normalize_safe", false, B(OP_ERO_NORMALIZE_SAFE).comp(c).ep(NV[i][0], NV[i][1]));
}
}
// fractal_noise_2d: the octave loop over grad_noised_2d, which is where a
// per-step rounding difference would accumulate if there were one.
for (int i = 0; i < NP; i++) {
for (int oct = 1; oct <= 7; oct += 2) {
for (uint32_t c = 0; c < 3; c++) {
add("fractal_noise_2d", false, B(OP_ERO_FRACTAL).comp(c)
.ep(P[i][0], P[i][1]).eoct(oct)
.es(0, 3.0f).es(1, 2.0f).es(2, 0.5f));
}
}
}
// phacelle_noise: sixteen samples around a cell, TAU written to nineteen
// digits -- more than float32 holds -- and 0.01111. The normalisation
// parameter is swept because 0 and 1 take different branches.
for (int i = 0; i < NP; i++) {
for (int k = 0; k < 2; k++) {
for (uint32_t c = 0; c < 4; c++) {
add("phacelle_noise", false, B(OP_ERO_PHACELLE).comp(c)
.ep(P[i][0], P[i][1]).en(0.6f, 0.8f, 0.0f)
.es(0, 1.0f).es(1, k ? 0.37f : 0.0f).es(2, k ? 1.0f : 0.0f));
}
}
}
// The filter and both wrappers, on the shipped bundles.
//
// rounding = (crease 0.1, ridge 0.0, input feature 0.1, per octave 2.0)
// onset = (gully 1.25, octave 1.25, ridge 2.8, ridge octave 1.5)
// assumed = (0.7, 1.0)
//
// max_octaves 999 is what every shipping caller passes and makes each
// octave weight exactly 1.0; 4.5 exercises the fractional taper, which is
// the branch that weights the last partial octave.
const float GN[][3] = {
{ 0.0f, 0.0f, 0.0f }, // flat: the assumed-slope path, since the
// real slope is zero and normalize_safe
// returns the zero vector
{ 0.3f, 0.1f, -0.2f },
{ -0.8f, 0.45f, 0.6f },
{ 0.95f, -1.2f, 0.7f },
};
const float MO[] = { 999.0f, 4.5f };
for (int i = 0; i < NP; i++) {
for (int g = 0; g < 4; g++) {
for (int m = 0; m < 2; m++) {
for (uint32_t c = 0; c < 5; c++) {
add("erosion_filter", false, B(OP_ERO_FILTER).comp(c)
.ep(P[i][0], P[i][1]).en(GN[g][0], GN[g][1], GN[g][2])
.es(0, 0.25f) // fade_target
.es(1, 0.16f) // strength
.es(2, 0.5f) // gully_weight
.es(3, 1.5f) // detail
.es(4, 5.0f) // scale
.es(5, 2.0f) // lacunarity
.es(6, 0.5f) // gain
.es(7, 1.0f) // cell_scale
.es(8, 0.5f) // normalization
.es(9, MO[m]) // max_octaves
.eoct(5)
.eround(0.1f, 0.0f, 0.1f, 2.0f)
.eonset(1.25f, 1.25f, 2.8f, 1.5f)
.eslope(0.7f, 1.0f), 5e-5);
}
for (uint32_t c = 0; c < 2; c++) {
add("erosion_node_eval_full", false, B(OP_ERO_EVAL_FULL).comp(c)
.ep(P[i][0], P[i][1]).en(GN[g][0], GN[g][1], GN[g][2])
.es(0, 5.0f) // height_scale
.es(1, 0.16f) // erosion_strength
.es(2, 1.0f) // cell_size
.es(3, 0.5f) // gain
.es(4, 0.6f) // peak_height_divisor (pn_fade_div)
.es(5, 0.15f) // gully_band_scale (pn_erosion_scale)
.es(6, 0.25f) // input_height_amplitude
.es(7, 2.0f) // octave_lacunarity
.es(8, 0.5f) // gully_depth
.es(9, 1.5f) // fine_gully_slope_restriction
.es(10, 0.5f) // gully_normalize_threshold
.es(11, -0.65f) // height_offset
.es(12, 0.0f) // height_offset_follows_altitude
.es(13, MO[m]) // max_octaves
.eoct(5)
.eround(0.1f, 0.0f, 0.1f, 2.0f)
.eonset(1.25f, 1.25f, 2.8f, 1.5f)
.eslope(0.7f, 1.0f), 5e-5);
}
// _grad exists because reconstructing the gradient by central
// differences cost six more full evaluations per edge -- 491 us
// to 3795 us per 32^3 grid. It is a different function from
// _full and has to be checked as one.
for (uint32_t c = 0; c < 4; c++) {
add("erosion_node_eval_grad", false, B(OP_ERO_EVAL_GRAD).comp(c)
.ep(P[i][0], P[i][1]).en(GN[g][0], GN[g][1], GN[g][2])
.es(0, 5.0f).es(1, 0.16f).es(2, 1.0f).es(3, 0.5f)
.es(4, 0.6f).es(5, 0.15f).es(6, 0.25f).es(7, 2.0f)
.es(8, 0.5f).es(9, 1.5f).es(10, 0.5f)
.es(11, -0.65f).es(12, 0.0f).es(13, MO[m])
.eoct(5)
.eround(0.1f, 0.0f, 0.1f, 2.0f)
.eonset(1.25f, 1.25f, 2.8f, 1.5f)
.eslope(0.7f, 1.0f), 5e-5);
}
}
}
}
}
static long peak_rss_kb() {
std::ifstream f("/proc/self/status");
std::string line;
while (std::getline(f, line)) {
if (line.compare(0, 6, "VmHWM:") == 0) {
return atol(line.c_str() + 6);
}
}
return -1;
}
int main(int argc, char **argv) {
cap_own_memory();
std::vector<std::pair<std::string, std::string> > subjects;
if (argc > 2) {
for (int i = 1; i + 1 < argc; i += 2) {
subjects.push_back(std::make_pair(std::string(argv[i]), std::string(argv[i + 1])));
}
} else {
for (int i = 0; i < SUBJECT_N; i++) {
subjects.push_back(std::make_pair(std::string(SUBJECT_FILES[i]),
std::string(SUBJECT_MARKS[i])));
}
}
std::string lib, subject_desc;
for (size_t i = 0; i < subjects.size(); i++) {
const std::string &file = subjects[i].first;
const bool raw_inc = file.size() > 4 && file.compare(file.size() - 4, 4, ".inc") == 0;
std::string text = read_file(file);
lib += raw_inc ? text : extract_raw_string(text, subjects[i].second);
lib += "\n";
subject_desc += (i ? "\n " : "") + file;
}
if (lib.find("bn_hash_bit_final") == std::string::npos ||
lib.find("bp_perlin_signed") == std::string::npos ||
lib.find("bv_voronoi_f1") == std::string::npos ||
lib.find("br_ramp_blend") == std::string::npos ||
lib.find("ps_simplex_noise_3d") == std::string::npos) {
fprintf(stderr, "FATAL: subject text is missing bn_hash_bit_final, bp_perlin_signed"
" or bv_voronoi_f1 -- wrong subject?\n");
return 1;
}
build_cases();
build_voronoi_cases();
build_ramp_cases();
build_simplex_cases();
build_erosion_cases();
if (g_cases.empty()) { fprintf(stderr, "FATAL: no cases\n"); return 1; }
int n_exact = 0;
for (size_t i = 0; i < g_cases.size(); i++) {
if (g_cases[i].exact) n_exact++;
}
// The opcode names come from the C++ enum, so the two dispatchers cannot
// drift by editing one of them.
std::string defs = "#version 450\n";
{
char b[128];
snprintf(b, sizeof(b), "#define STRIDE_W %uu\n#define CASE_COUNT %uu\n",
STRIDE, (unsigned)g_cases.size());
defs += b;
for (size_t i = 0; i < sizeof(OP_NAMES) / sizeof(OP_NAMES[0]); i++) {
snprintf(b, sizeof(b), "#define %s %uu\n", OP_NAMES[i].name, (unsigned)OP_NAMES[i].op);
defs += b;
}
// And the erosion block's slots, for the same reason.
for (size_t i = 0; i < sizeof(SLOT_NAMES) / sizeof(SLOT_NAMES[0]); i++) {
snprintf(b, sizeof(b), "#define %s %uu\n", SLOT_NAMES[i].name,
(unsigned)SLOT_NAMES[i].slot);
defs += b;
}
}
const std::string shader = defs + lib + DISPATCH_GLSL;
Ctx c{};
init_vulkan(c);
printf("device : %s\n", c.devname);
printf("subject: %s (%zu bytes of GLSL, extracted not retyped)\n",
subject_desc.c_str(), lib.size());
printf("cases : %d EXACT (bit patterns, no epsilon), %d NEAR (rel %.0e)\n",
n_exact, (int)g_cases.size() - n_exact, TOL_REL);
printf("shader : %zu bytes, and its size does NOT depend on the case count\n\n",
shader.size());
std::vector<uint32_t> spv = compile_glsl(shader);
// inputs
std::vector<uint32_t> inwords;
inwords.reserve(g_cases.size() * STRIDE);
for (size_t i = 0; i < g_cases.size(); i++) {
for (uint32_t k = 0; k < STRIDE; k++) {
inwords.push_back(g_cases[i].w[k]);
}
}
VkBuffer buf, inbuf;
VkDeviceMemory mem, inmem;
void *mapped = nullptr;
void *inmapped = nullptr;
const VkDeviceSize out_bytes = sizeof(uint32_t) * g_cases.size();
const VkDeviceSize in_bytes = sizeof(uint32_t) * inwords.size();
{
VkBufferCreateInfo bi{VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO};
bi.size = out_bytes;
bi.usage = VK_BUFFER_USAGE_STORAGE_BUFFER_BIT;
VK(vkCreateBuffer(c.dev, &bi, nullptr, &buf));
VkMemoryRequirements mr;
vkGetBufferMemoryRequirements(c.dev, buf, &mr);
VkMemoryAllocateInfo ai{VK_STRUCTURE_TYPE_MEMORY_ALLOCATE_INFO};
ai.allocationSize = mr.size;
ai.memoryTypeIndex = find_mem(c, mr.memoryTypeBits,
VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT);
VK(vkAllocateMemory(c.dev, &ai, nullptr, &mem));
VK(vkBindBufferMemory(c.dev, buf, mem, 0));
VK(vkMapMemory(c.dev, mem, 0, VK_WHOLE_SIZE, 0, &mapped));
memset(mapped, 0, out_bytes);
bi.size = in_bytes;
VK(vkCreateBuffer(c.dev, &bi, nullptr, &inbuf));
vkGetBufferMemoryRequirements(c.dev, inbuf, &mr);
ai.allocationSize = mr.size;
ai.memoryTypeIndex = find_mem(c, mr.memoryTypeBits,
VK_MEMORY_PROPERTY_HOST_VISIBLE_BIT | VK_MEMORY_PROPERTY_HOST_COHERENT_BIT);
VK(vkAllocateMemory(c.dev, &ai, nullptr, &inmem));
VK(vkBindBufferMemory(c.dev, inbuf, inmem, 0));
VK(vkMapMemory(c.dev, inmem, 0, VK_WHOLE_SIZE, 0, &inmapped));
memcpy(inmapped, inwords.data(), in_bytes);
}
VkDescriptorSetLayoutBinding bind[2]{};
for (int i = 0; i < 2; i++) {
bind[i].binding = (uint32_t)i;
bind[i].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
bind[i].descriptorCount = 1;
bind[i].stageFlags = VK_SHADER_STAGE_COMPUTE_BIT;
}
VkDescriptorSetLayoutCreateInfo dli{VK_STRUCTURE_TYPE_DESCRIPTOR_SET_LAYOUT_CREATE_INFO};
dli.bindingCount = 2; dli.pBindings = bind;
VkDescriptorSetLayout dsl;
VK(vkCreateDescriptorSetLayout(c.dev, &dli, nullptr, &dsl));
VkPipelineLayoutCreateInfo pli{VK_STRUCTURE_TYPE_PIPELINE_LAYOUT_CREATE_INFO};
pli.setLayoutCount = 1; pli.pSetLayouts = &dsl;
VkPipelineLayout playout;
VK(vkCreatePipelineLayout(c.dev, &pli, nullptr, &playout));
VkShaderModuleCreateInfo smi{VK_STRUCTURE_TYPE_SHADER_MODULE_CREATE_INFO};
smi.codeSize = spv.size() * 4; smi.pCode = spv.data();
VkShaderModule mod;
VK(vkCreateShaderModule(c.dev, &smi, nullptr, &mod));
VkPipelineShaderStageCreateInfo st{VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO};
st.stage = VK_SHADER_STAGE_COMPUTE_BIT; st.module = mod; st.pName = "main";
VkComputePipelineCreateInfo cpi{VK_STRUCTURE_TYPE_COMPUTE_PIPELINE_CREATE_INFO};
cpi.stage = st; cpi.layout = playout;
VkPipeline pipe;
VK(vkCreateComputePipelines(c.dev, VK_NULL_HANDLE, 1, &cpi, nullptr, &pipe));
VkDescriptorPoolSize ps{VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 2};
VkDescriptorPoolCreateInfo dpi{VK_STRUCTURE_TYPE_DESCRIPTOR_POOL_CREATE_INFO};
dpi.maxSets = 1; dpi.poolSizeCount = 1; dpi.pPoolSizes = &ps;
VkDescriptorPool dpool;
VK(vkCreateDescriptorPool(c.dev, &dpi, nullptr, &dpool));
VkDescriptorSetAllocateInfo dai{VK_STRUCTURE_TYPE_DESCRIPTOR_SET_ALLOCATE_INFO};
dai.descriptorPool = dpool; dai.descriptorSetCount = 1; dai.pSetLayouts = &dsl;
VkDescriptorSet dset;
VK(vkAllocateDescriptorSets(c.dev, &dai, &dset));
VkDescriptorBufferInfo dbi{buf, 0, VK_WHOLE_SIZE};
VkDescriptorBufferInfo dbi_in{inbuf, 0, VK_WHOLE_SIZE};
VkWriteDescriptorSet w[2]{};
for (int i = 0; i < 2; i++) {
w[i].sType = VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET;
w[i].dstSet = dset; w[i].dstBinding = (uint32_t)i; w[i].descriptorCount = 1;
w[i].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
w[i].pBufferInfo = (i == 0) ? &dbi : &dbi_in;
}
vkUpdateDescriptorSets(c.dev, 2, w, 0, nullptr);
VkCommandBufferAllocateInfo cbi{VK_STRUCTURE_TYPE_COMMAND_BUFFER_ALLOCATE_INFO};
cbi.commandPool = c.cmdpool; cbi.level = VK_COMMAND_BUFFER_LEVEL_PRIMARY; cbi.commandBufferCount = 1;
VkCommandBuffer cmd;
VK(vkAllocateCommandBuffers(c.dev, &cbi, &cmd));
VkCommandBufferBeginInfo bbi{VK_STRUCTURE_TYPE_COMMAND_BUFFER_BEGIN_INFO};
bbi.flags = VK_COMMAND_BUFFER_USAGE_ONE_TIME_SUBMIT_BIT;
VK(vkBeginCommandBuffer(cmd, &bbi));
vkCmdBindPipeline(cmd, VK_PIPELINE_BIND_POINT_COMPUTE, pipe);
vkCmdBindDescriptorSets(cmd, VK_PIPELINE_BIND_POINT_COMPUTE, playout, 0, 1, &dset, 0, nullptr);
vkCmdDispatch(cmd, ((uint32_t)g_cases.size() + 63u) / 64u, 1, 1);
VK(vkEndCommandBuffer(cmd));
VkSubmitInfo si{VK_STRUCTURE_TYPE_SUBMIT_INFO};
si.commandBufferCount = 1; si.pCommandBuffers = &cmd;
VK(vkQueueSubmit(c.queue, 1, &si, VK_NULL_HANDLE));
VK(vkQueueWaitIdle(c.queue));
const uint32_t *got = (const uint32_t *)mapped;
int bad = 0, bad_exact = 0, n_bitident = 0;
double worst_rel = 0.0;
std::string worst_name;
for (size_t i = 0; i < g_cases.size(); i++) {
const Rec &k = g_cases[i];
bool ok;
if (k.exact) {
ok = (got[i] == k.cpu);
if (!ok) bad_exact++;
} else {
const double a = float_of(k.cpu), b = float_of(got[i]);
const bool bitident = (got[i] == k.cpu);
if (bitident) n_bitident++;
// IDENTICAL BITS ARE AGREEMENT, and that clause has to come first.
//
// Without it two identical INFINITIES read as a disagreement:
// inf - inf is nan, and every comparison against nan is false, so
// the arm below rejects a pair that is the same number. Found
// 2026-09-06 by smooth_start(-0.5, 0.0), which divides by the
// smoothing and produces +inf on both sides -- the gate printed
// `cpu=inf gpu=inf FAIL`, which is a sentence about the comparator
// and not about the filter.
const double diff = fabs(a - b);
const double scale = fabs(a) > 1.0 ? fabs(a) : 1.0;
const double rel = diff / scale;
// A non-finite rel must not become the reported worst: it would
// print as the headline number for a pair that agrees.
if (std::isfinite(rel) && rel > worst_rel) {
worst_rel = rel; worst_name = k.name;
}
const double abs_tol = (k.tol_abs > 0.0) ? k.tol_abs : TOL_ABS;
ok = bitident || (diff <= abs_tol) || (rel <= TOL_REL);
}
if (!ok) {
bad++;
if (bad <= 20) {
printf(" FAIL %-30s dim=%u type=%u comp=%u cpu=%.9g gpu=%.9g\n",
k.name.c_str(), k.w[W_DIM], k.w[W_TYPE], k.w[W_COMP],
float_of(k.cpu), float_of(got[i]));
}
}
}
const int n_near = (int)g_cases.size() - n_exact;
printf("\n%zu/%zu cases agree (%d exact, %d of those failing)\n",
g_cases.size() - (size_t)bad, g_cases.size(), n_exact, bad_exact);
if (n_near > 0) {
printf("NEAR: %d/%d bit-identical anyway; worst relative error %.3g (%s)\n",
n_bitident, n_near, worst_rel, worst_name.c_str());
}
printf("peak RSS: %.1f MiB\n", peak_rss_kb() / 1024.0);
return bad ? 1 : 0;
}