// Standalone throughput probe, not VoxelForge's generator or an Unreal integration.
// Build from a VS 2022 x64 Native Tools prompt:
//   cl /nologo /O2 /std:c++20 /EHsc /arch:AVX2 brick_kernel.cpp /Fe:brick_kernel.exe
#include <algorithm>
#include <array>
#include <chrono>
#include <cstdlib>
#include <cstdint>
#include <cstdio>
#include <thread>
#include <vector>

constexpr int N = 256, B = 8, NB = N / B;
#ifndef VF_WORKERS
#define VF_WORKERS 8
#endif
struct Sphere { int x, y, z, r; };
static std::array<Sphere, 48> spheres;
static std::vector<int16_t> field(size_t(N) * N * N);
static std::vector<uint8_t> reference_sign(size_t(N) * N * N);
static volatile uint64_t sink = 0;

static uint32_t hash3(int x, int y, int z, uint32_t seed) {
    uint32_t h = seed ^ (uint32_t(x) * 0x9e3779b1u) ^
        (uint32_t(y) * 0x85ebca77u) ^ (uint32_t(z) * 0xc2b2ae3du);
    h ^= h >> 16; h *= 0x7feb352du; h ^= h >> 15;
    h *= 0x846ca68bu; return h ^ (h >> 16);
}
// Integer lattice noise: fixed shifts and products; no FP/CRT/FMA.
static int noise(int x, int y, int z) {
    int sum = 0;
    for (int octave = 0; octave != 3; ++octave) {
        int s = 16 >> octave, ix = x / s, iy = y / s, iz = z / s;
        int fx = (x % s) * 256 / s, fy = (y % s) * 256 / s, fz = (z % s) * 256 / s;
        uint32_t seed = 1337u + uint32_t(octave) * 17u;
        int yz[2];
        for (int dz = 0; dz != 2; ++dz) {
            int yy[2];
            for (int dy = 0; dy != 2; ++dy) {
                int a = int(hash3(ix, iy + dy, iz + dz, seed) & 255u);
                int b = int(hash3(ix + 1, iy + dy, iz + dz, seed) & 255u);
                yy[dy] = a + ((b - a) * fx >> 8);
            }
            yz[dz] = yy[0] + ((yy[1] - yy[0]) * fy >> 8);
        }
        int v = yz[0] + ((yz[1] - yz[0]) * fz >> 8);
        sum += (v - 128) >> octave;
    }
    return sum; // in [-224, 222], conservatively use +/-224.
}
static int density(int x, int y, int z) {
    int d = -10000;
    for (const Sphere& s : spheres) {
        int dx=x-s.x, dy=y-s.y, dz=z-s.z;
        int q=s.r*s.r-dx*dx-dy*dy-dz*dz;
        d=std::max(d,q);
    }
    return d + noise(x,y,z);
}
static size_t at(int x,int y,int z) { return (size_t(z)*N+y)*N+x; }
template<class F> static double run(F&& f) {
    auto a=std::chrono::steady_clock::now(); f();
    return std::chrono::duration<double,std::milli>(std::chrono::steady_clock::now()-a).count();
}
static void parallel_z(auto fn) {
    std::array<std::thread,VF_WORKERS> t;
    for (int i=0;i<VF_WORKERS;++i) t[i]=std::thread([=]{
        fn((i*NB/VF_WORKERS)*B, ((i+1)*NB/VF_WORKERS)*B);
    });
    for (auto& x:t) x.join();
}
int main() {
    std::printf("workers=%d\n", VF_WORKERS);
    for(size_t i=0;i<spheres.size();++i) {
        uint32_t h=hash3(int(i),17,71,12345);
        spheres[i]={int((h&255u)),int((h>>8)&255u),int((h>>16)&255u), 12+int((h>>24)%16u)};
    }
    for(int trial=0;trial<3;++trial) {
        double full=run([]{parallel_z([](int lo,int hi){
            for(int z=lo;z<hi;++z) for(int y=0;y<N;++y) for(int x=0;x<N;++x)
                field[at(x,y,z)]=int16_t(density(x,y,z));
        });});
        uint64_t checksum=0; for(size_t i=0;i<field.size();i+=4093) checksum+=uint16_t(field[i]); sink=checksum;
        std::printf("full trial=%d ms=%.3f checksum=%llu\n",trial,full,(unsigned long long)checksum);
    }
    for(size_t i=0;i<field.size();++i) reference_sign[i]=field[i]>0;
    // Conservative sphere/brick bounds with a +/-224 noise interval.
    std::array<uint8_t, size_t(NB)*NB*NB> kind{};
    size_t solid=0,air=0,mixed=0;
    double classify=run([&]{
        for(int bz=0;bz<NB;++bz) for(int by=0;by<NB;++by) for(int bx=0;bx<NB;++bx) {
            int x0=bx*B,y0=by*B,z0=bz*B,x1=x0+B-1,y1=y0+B-1,z1=z0+B-1;
            bool possiblyAir=false, definitelyAir=false;
            for(const Sphere& s:spheres) {
                int dx=std::max({x0-s.x,0,s.x-x1});
                int dy=std::max({y0-s.y,0,s.y-y1});
                int dz=std::max({z0-s.z,0,s.z-z1});
                int near2=dx*dx+dy*dy+dz*dz;
                int fx=std::max(std::abs(x0-s.x),std::abs(x1-s.x));
                int fy=std::max(std::abs(y0-s.y),std::abs(y1-s.y));
                int fz=std::max(std::abs(z0-s.z),std::abs(z1-s.z));
                int far2=fx*fx+fy*fy+fz*fz;
                possiblyAir |= (s.r*s.r-near2+224>=0);
                definitelyAir |= (s.r*s.r-far2-224>0);
            }
            uint8_t k=definitelyAir?2:(possiblyAir?1:0);
            kind[(size_t(bz)*NB+by)*NB+bx]=k;
            if(k==0)++solid; else if(k==2)++air; else ++mixed;
        }
    });
    std::printf("classify ms=%.3f solid=%zu air=%zu mixed=%zu total=%zu\n",classify,solid,air,mixed,kind.size());
    for(int trial=0;trial<3;++trial) {
        double sparse=run([&]{parallel_z([&](int lo,int hi){
            for(int bz=lo/B;bz<hi/B;++bz) for(int by=0;by<NB;++by) for(int bx=0;bx<NB;++bx) {
                uint8_t k=kind[(size_t(bz)*NB+by)*NB+bx];
                for(int z=bz*B;z<(bz+1)*B;++z) for(int y=by*B;y<(by+1)*B;++y)
                    for(int x=bx*B;x<(bx+1)*B;++x)
                        field[at(x,y,z)]= k==1?int16_t(density(x,y,z)):int16_t(k==2?1:-1);
            }
        });});
        uint64_t checksum=0; for(size_t i=0;i<field.size();i+=4093) checksum+=uint16_t(field[i]); sink=checksum;
        std::printf("sparse trial=%d ms=%.3f checksum=%llu\n",trial,sparse,(unsigned long long)checksum);
    }
    size_t mismatches=0;
    for(size_t i=0;i<field.size();++i) mismatches+=(reference_sign[i]!=(field[i]>0));
    std::printf("sign mismatches=%zu of %zu\n",mismatches,field.size());
}