// Standalone research probe. Synthetic operator composition inspired by VoxelForge;
// NOT extracted production code. Compile with MSVC /O2 /arch:AVX2 /std:c++20.
#include <immintrin.h>
#include <algorithm>
#include <array>
#include <chrono>
#include <cstdint>
#include <cstdio>
#include <thread>
#include <vector>
constexpr int N=256, B=8, W=6;
using V=__m256i;
static std::vector<int16_t> scalarField(size_t(N)*N*N), vectorField(size_t(N)*N*N);
alignas(32) static int32_t fadeTable[256];
static inline uint32_t hashS(uint32_t x,uint32_t y,uint32_t z,uint32_t seed){
 uint32_t h=seed^(x*0x9e3779b1u)^(y*0x85ebca77u)^(z*0xc2b2ae3du);
 h^=h>>16; h*=0x7feb352du; h^=h>>15; h*=0x846ca68bu; return h^(h>>16);
}
static inline V hashV(V x,V y,V z,V seed){
 V h=_mm256_xor_si256(seed,_mm256_xor_si256(_mm256_mullo_epi32(x,_mm256_set1_epi32(0x9e3779b1u)),
 _mm256_xor_si256(_mm256_mullo_epi32(y,_mm256_set1_epi32(0x85ebca77u)),_mm256_mullo_epi32(z,_mm256_set1_epi32(0xc2b2ae3du)))));
 h=_mm256_xor_si256(h,_mm256_srli_epi32(h,16)); h=_mm256_mullo_epi32(h,_mm256_set1_epi32(0x7feb352du));
 h=_mm256_xor_si256(h,_mm256_srli_epi32(h,15)); h=_mm256_mullo_epi32(h,_mm256_set1_epi32(0x846ca68bu));
 return _mm256_xor_si256(h,_mm256_srli_epi32(h,16));
}
static inline int lerpS(int a,int b,int t){return a+((b-a)*t>>8);}
static inline V lerpV(V a,V b,V t){return _mm256_add_epi32(a,_mm256_srai_epi32(_mm256_mullo_epi32(_mm256_sub_epi32(b,a),t),8));}
static inline int gradS(uint32_t h,int x,int y,int z){
 h&=15;int u=(h&8)?y:x;int v=(h&12)==0?y:((h&13)==12?x:z);
 return ((h&1)?-u:u)+((h&2)?-v:v);
}
static inline V selectV(V mask,V a,V b){return _mm256_blendv_epi8(b,a,mask);}
static inline V gradV(V h,V x,V y,V z){
 V zero=_mm256_setzero_si256(); h=_mm256_and_si256(h,_mm256_set1_epi32(15));
 V u=selectV(_mm256_cmpeq_epi32(_mm256_and_si256(h,_mm256_set1_epi32(8)),zero),x,y);
 V isY=_mm256_cmpeq_epi32(_mm256_and_si256(h,_mm256_set1_epi32(12)),zero);
 V isX=_mm256_cmpeq_epi32(_mm256_and_si256(h,_mm256_set1_epi32(13)),_mm256_set1_epi32(12));
 V v=selectV(isY,y,selectV(isX,x,z));
 V negU=_mm256_cmpeq_epi32(_mm256_and_si256(h,_mm256_set1_epi32(1)),_mm256_set1_epi32(1));
 V negV=_mm256_cmpeq_epi32(_mm256_and_si256(h,_mm256_set1_epi32(2)),_mm256_set1_epi32(2));
 return _mm256_add_epi32(_mm256_sub_epi32(_mm256_xor_si256(u,negU),negU),_mm256_sub_epi32(_mm256_xor_si256(v,negV),negV));
}
static inline int noiseS(int x,int y,int z,int shift,int seed){
 int ix=x>>shift,iy=y>>shift,iz=z>>shift,fx=(x&((1<<shift)-1))<<(8-shift),fy=(y&((1<<shift)-1))<<(8-shift),fz=(z&((1<<shift)-1))<<(8-shift);
 int ux=fadeTable[fx],uy=fadeTable[fy],uz=fadeTable[fz];
 int yz[2]; for(int dz=0;dz<2;++dz){int yy[2];for(int dy=0;dy<2;++dy){
 int a=gradS(hashS(ix,iy+dy,iz+dz,seed),fx,fy-dy*256,fz-dz*256);
 int b=gradS(hashS(ix+1,iy+dy,iz+dz,seed),fx-256,fy-dy*256,fz-dz*256);
 yy[dy]=lerpS(a,b,ux);}
 yz[dz]=lerpS(yy[0],yy[1],uy);} return lerpS(yz[0],yz[1],uz)/4;
}
static inline V noiseV(V x,V y,V z,int shift,int seed){
 V mask=_mm256_set1_epi32((1<<shift)-1), one=_mm256_set1_epi32(1);
 V ix=_mm256_srai_epi32(x,shift),iy=_mm256_srai_epi32(y,shift),iz=_mm256_srai_epi32(z,shift);
 V fx=_mm256_slli_epi32(_mm256_and_si256(x,mask),8-shift),fy=_mm256_slli_epi32(_mm256_and_si256(y,mask),8-shift),fz=_mm256_slli_epi32(_mm256_and_si256(z,mask),8-shift);
 V ux=_mm256_i32gather_epi32(fadeTable,fx,4),uy=_mm256_i32gather_epi32(fadeTable,fy,4),uz=_mm256_i32gather_epi32(fadeTable,fz,4);
 V yz[2]; for(int dz=0;dz<2;++dz){V yy[2];for(int dy=0;dy<2;++dy){
 V y1=_mm256_add_epi32(iy,_mm256_set1_epi32(dy)),z1=_mm256_add_epi32(iz,_mm256_set1_epi32(dz));
 V yyoff=_mm256_sub_epi32(fy,_mm256_set1_epi32(dy*256)),zzoff=_mm256_sub_epi32(fz,_mm256_set1_epi32(dz*256));
 V a=gradV(hashV(ix,y1,z1,_mm256_set1_epi32(seed)),fx,yyoff,zzoff);
 V b=gradV(hashV(_mm256_add_epi32(ix,one),y1,z1,_mm256_set1_epi32(seed)),_mm256_sub_epi32(fx,_mm256_set1_epi32(256)),yyoff,zzoff);
 yy[dy]=lerpV(a,b,ux);}yz[dz]=lerpV(yy[0],yy[1],uy);}
 V raw=lerpV(yz[0],yz[1],uz);V bias=_mm256_and_si256(_mm256_srai_epi32(raw,31),_mm256_set1_epi32(3));
 return _mm256_srai_epi32(_mm256_add_epi32(raw,bias),2);
}
static inline int smoothMaxS(int a,int b){int h=std::max(0,64-std::abs(a-b));return std::max(a,b)+(h*h>>8);}
static inline V smoothMaxV(V a,V b){V d=_mm256_sub_epi32(a,b);V h=_mm256_max_epi32(_mm256_setzero_si256(),_mm256_sub_epi32(_mm256_set1_epi32(64),_mm256_abs_epi32(d)));
 return _mm256_add_epi32(_mm256_max_epi32(a,b),_mm256_srai_epi32(_mm256_mullo_epi32(h,h),8));}
static inline V evalV(V x,V y,V z){
 V c=_mm256_set1_epi32(1024), warp=_mm256_srai_epi32(noiseV(_mm256_add_epi32(x,c),_mm256_add_epi32(y,c),_mm256_add_epi32(z,c),4,17),4);
 V nx=_mm256_add_epi32(_mm256_add_epi32(x,c),warp),ny=_mm256_sub_epi32(_mm256_add_epi32(y,c),warp),nz=_mm256_add_epi32(z,c);
 V fb=_mm256_add_epi32(noiseV(nx,ny,nz,4,100),_mm256_add_epi32(_mm256_srai_epi32(noiseV(nx,ny,nz,3,101),1),_mm256_srai_epi32(noiseV(nx,ny,nz,2,102),2)));
 V ridge=_mm256_sub_epi32(_mm256_set1_epi32(128),_mm256_abs_epi32(noiseV(nx,ny,nz,3,202)));
 // Scalar signed division truncates toward zero; AVX2 shift rounds negative down.
 // Correct sign before shift to make the two exact.
 auto divPow2=[](V a,int shift){V bias=_mm256_and_si256(_mm256_srai_epi32(a,31),_mm256_set1_epi32((1<<shift)-1));return _mm256_srai_epi32(_mm256_add_epi32(a,bias),shift);};
 // Exact truncation toward zero matches scalar division by powers of two.
 V floor=_mm256_add_epi32(_mm256_sub_epi32(z,_mm256_set1_epi32(84)),divPow2(fb,2));
 V ceil=_mm256_add_epi32(_mm256_sub_epi32(_mm256_set1_epi32(170),z),divPow2(ridge,2));
 V stratum=_mm256_min_epi32(floor,ceil),dy=_mm256_sub_epi32(y,_mm256_set1_epi32(128)),dz=_mm256_sub_epi32(z,_mm256_set1_epi32(128));
 V dy2=_mm256_mullo_epi32(dy,dy),dz2=_mm256_mullo_epi32(dz,dz),dx=_mm256_sub_epi32(x,_mm256_set1_epi32(128));
 V passage=divPow2(_mm256_sub_epi32(_mm256_set1_epi32(18*18),_mm256_add_epi32(dy2,dz2)),3);
 V room=divPow2(_mm256_sub_epi32(_mm256_set1_epi32(35*35),_mm256_add_epi32(_mm256_mullo_epi32(dx,dx),_mm256_add_epi32(dy2,dz2))),4);
 return smoothMaxV(stratum,smoothMaxV(passage,room));
}
static inline int evalS(int x,int y,int z){
 int warp=noiseS(x+1024,y+1024,z+1024,4,17)>>4; int nx=x+1024+warp,ny=y+1024-warp,nz=z+1024;
 int fb=noiseS(nx,ny,nz,4,100)+(noiseS(nx,ny,nz,3,101)>>1)+(noiseS(nx,ny,nz,2,102)>>2);
 int ridge=128-std::abs(noiseS(nx,ny,nz,3,202));
 int stratum=std::min(z-84+fb/4,170-z+ridge/4);
 int dy=y-128,dz=z-128,dx=x-128;
 int passage=(18*18-dy*dy-dz*dz)/8,room=(35*35-dx*dx-dy*dy-dz*dz)/16;
 return smoothMaxS(stratum,smoothMaxS(passage,room));
}
static inline size_t idx(int x,int y,int z){return(size_t(z)*N+y)*N+x;}
template<class F> static double ms(F&& f){auto t=std::chrono::steady_clock::now();f();return std::chrono::duration<double,std::milli>(std::chrono::steady_clock::now()-t).count();}
template<class F> static void parallel(F fn){std::array<std::thread,W> threads;for(int k=0;k<W;++k)threads[k]=std::thread([=]{fn(k*N/W,(k+1)*N/W);});for(auto& t:threads)t.join();}
static void fillScalar(){parallel([](int a,int b){for(int z=a;z<b;++z)for(int y=0;y<N;++y)for(int x=0;x<N;++x)scalarField[idx(x,y,z)]=(int16_t)evalS(x,y,z);});}
static void fillVector(){parallel([](int a,int b){for(int z=a;z<b;++z)for(int y=0;y<N;++y)for(int x=0;x<N;x+=8){
 V xx=_mm256_add_epi32(_mm256_set1_epi32(x),_mm256_setr_epi32(0,1,2,3,4,5,6,7)); V d=evalV(xx,_mm256_set1_epi32(y),_mm256_set1_epi32(z));
 alignas(32) int32_t values[8];_mm256_store_si256((V*)values,d);for(int j=0;j<8;++j)vectorField[idx(x+j,y,z)]=(int16_t)values[j];}});}
static inline int dist2box(int a,int lo,int hi){int d=std::max({lo-a,0,a-hi});return d*d;}
int main(){
 for(int t=0;t<256;++t){int64_t q=t;int64_t q2=q*q,q3=q2*q,q4=q3*q,q5=q4*q;
 fadeTable[t]=int32_t((6*q5-15*q4*256+10*q3*256*256)/(256LL*256*256*256));}
 for(int t=0;t<3;++t){double s=ms(fillScalar),v=ms(fillVector);size_t mismatch=0;for(size_t i=0;i<scalarField.size();++i)mismatch+=scalarField[i]!=vectorField[i];
 std::printf("full trial=%d scalar_ms=%.3f avx2_ms=%.3f speedup=%.3f mismatch=%zu\n",t,s,v,s/v,mismatch);}
 struct Point{int x,y,z;};std::vector<Point> points(200000);uint32_t rng=0x12345678;
 for(auto& p:points){rng^=rng<<13;rng^=rng>>17;rng^=rng<<5;p.x=int(rng&255);
 rng^=rng<<13;rng^=rng>>17;rng^=rng<<5;p.y=int(rng&255);
 rng^=rng<<13;rng^=rng>>17;rng^=rng<<5;p.z=int(rng&255);}
 for(int t=0;t<3;++t){int64_t checksum=0;double q=ms([&]{for(const auto& p:points)checksum+=evalS(p.x,p.y,p.z);});
 std::printf("random_query trial=%d count=%zu ms=%.3f ns_per_query=%.1f checksum=%lld\n",t,points.size(),q,q*1e6/points.size(),(long long)checksum);}
 // Radius 10m at 25cm grid, center on tile boundary. Count exact AABB touched bricks.
 const int cx=128,cy=128,cz=128,r=40;int touched=0;for(int bz=0;bz<N/B;++bz)for(int by=0;by<N/B;++by)for(int bx=0;bx<N/B;++bx){
 int x0=bx*B,x1=x0+B-1,y0=by*B,y1=y0+B-1,z0=bz*B,z1=z0+B-1;
 if(dist2box(cx,x0,x1)+dist2box(cy,y0,y1)+dist2box(cz,z0,z1)<=r*r)++touched;}
 std::printf("crater radius_vox=%d touched_bricks=%d tile_AABB_max=%d samples_AABB=%d\n",r,touched,64,81*81*81);
 // One logical transaction: versioned edit, canonical max composition, copy-on-write buffer.
 std::vector<int16_t> next(vectorField.size());double copy=ms([&]{next=vectorField;});size_t changed=0;double edit=ms([&]{parallel([&](int za,int zb){for(int z=std::max(za,cz-r);z<std::min(zb,cz+r+1);++z)for(int y=cy-r;y<=cy+r;++y)for(int x=cx-r;x<=cx+r;++x){
 int dx=x-cx,dy=y-cy,dz=z-cz,q=r*r-dx*dx-dy*dy-dz*dz;
 if(q>=0){int d=std::max<int>(vectorField[idx(x,y,z)],q/16);next[idx(x,y,z)]=(int16_t)d;}}});});
 for(int z=cz-r;z<=cz+r;++z)for(int y=cy-r;y<=cy+r;++y)for(int x=cx-r;x<=cx+r;++x)changed+=next[idx(x,y,z)]!=vectorField[idx(x,y,z)];
 std::printf("crater copy_ms=%.3f edit_ms=%.3f changed_lattice=%zu copy_bytes=%zu\n",copy,edit,changed,vectorField.size()*sizeof(int16_t));
}