alice-sdf 1.9.2

//! SIMD-accelerated SDF evaluation (8 points at once)
//!
//! This module provides 8-wide SIMD evaluation using AVX2/AVX-512/NEON,
//! achieving up to 8x speedup over scalar evaluation.
//!
//! Author: Moroya Sakamoto

use super::compiler::CompiledSdf;
use super::opcode::OpCode;
use super::simd::{Quatx8, Vec3x8};
use glam::Vec3;
use wide::{f32x8, CmpGt, CmpLt};

use crate::modifiers::perlin_noise_3d;
use crate::operations::{
    sdf_columns_intersection, sdf_columns_subtraction, sdf_columns_union,
    sdf_exp_smooth_intersection, sdf_exp_smooth_subtraction, sdf_exp_smooth_union,
};
use crate::primitives::*;

/// Maximum stack depth for value stack (SIMD)
const MAX_VALUE_STACK: usize = 64;
/// Maximum stack depth for coordinate transforms (SIMD)
const MAX_COORD_STACK: usize = 32;

/// Coordinate frame for SIMD evaluation
#[derive(Clone, Copy)]
struct CoordFrameSimd {
    point: Vec3x8,
    scale_correction: f32x8,
    opcode: OpCode,
    params: [f32; 4],
    aux_offset: u32,
    aux_len: u32,
    /// Per-lane distance divisor applied at PopTransform (LatticeDeform Jacobian)
    lane_correction: f32x8,
}

impl Default for CoordFrameSimd {
    fn default() -> Self {
        Self {
            point: Vec3x8::zero(),
            scale_correction: f32x8::ONE,
            opcode: OpCode::End,
            params: [0.0; 4],
            aux_offset: 0,
            aux_len: 0,
            lane_correction: f32x8::ONE,
        }
    }
}

/// Evaluate compiled SDF at 8 points simultaneously
///
/// This is the core SIMD evaluation function. It processes 8 points
/// in parallel using AVX2/AVX-512 instructions.
///
/// # Arguments
/// * `sdf` - Compiled SDF bytecode
/// * `points` - 8 points packed as Vec3x8
///
/// # Returns
/// 8 distance values as f32x8
#[inline]
pub fn eval_compiled_simd(sdf: &CompiledSdf, points: Vec3x8) -> f32x8 {
    let mut value_stack: [f32x8; MAX_VALUE_STACK] = [f32x8::ZERO; MAX_VALUE_STACK];
    let mut vsp: usize = 0;

    let mut coord_stack: [CoordFrameSimd; MAX_COORD_STACK] =
        [CoordFrameSimd::default(); MAX_COORD_STACK];
    let mut csp: usize = 0;

    let mut p = points;
    let mut scale_correction = f32x8::ONE;

    for inst in &sdf.instructions {
        match inst.opcode {
            // === Primitives ===
            OpCode::Sphere => {
                let r = f32x8::splat(inst.params[0]);
                let d = p.length() - r;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Box3d => {
                let hx = f32x8::splat(inst.params[0]);
                let hy = f32x8::splat(inst.params[1]);
                let hz = f32x8::splat(inst.params[2]);

                // q = abs(p) - half_extents
                let qx = p.x.abs() - hx;
                let qy = p.y.abs() - hy;
                let qz = p.z.abs() - hz;

                // length(max(q, 0)) + min(max(q.x, q.y, q.z), 0)
                let qx_pos = qx.max(f32x8::ZERO);
                let qy_pos = qy.max(f32x8::ZERO);
                let qz_pos = qz.max(f32x8::ZERO);
                let outside = (qx_pos * qx_pos + qy_pos * qy_pos + qz_pos * qz_pos).sqrt();
                let inside = qx.max(qy).max(qz).min(f32x8::ZERO);
                let d = outside + inside;

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Cylinder => {
                let r = f32x8::splat(inst.params[0]);
                let h = f32x8::splat(inst.params[1]);

                // d.x = length(p.xz) - r
                // d.y = abs(p.y) - h
                let dx = (p.x * p.x + p.z * p.z).sqrt() - r;
                let dy = p.y.abs() - h;

                let dx_pos = dx.max(f32x8::ZERO);
                let dy_pos = dy.max(f32x8::ZERO);
                let outside = (dx_pos * dx_pos + dy_pos * dy_pos).sqrt();
                let inside = dx.max(dy).min(f32x8::ZERO);
                let d = outside + inside;

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Torus => {
                let major = f32x8::splat(inst.params[0]);
                let minor = f32x8::splat(inst.params[1]);

                // q = vec2(length(p.xz) - major, p.y)
                let qx = (p.x * p.x + p.z * p.z).sqrt() - major;
                let qy = p.y;
                let d = (qx * qx + qy * qy).sqrt() - minor;

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Plane => {
                let nx = f32x8::splat(inst.params[0]);
                let ny = f32x8::splat(inst.params[1]);
                let nz = f32x8::splat(inst.params[2]);
                let dist = f32x8::splat(inst.params[3]);

                // dot(p, n) - distance: same law as `sdf_plane` (tree / scalar / BVH)
                let d = p.x * nx + p.y * ny + p.z * nz - dist;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Capsule => {
                let ax = f32x8::splat(inst.params[0]);
                let ay = f32x8::splat(inst.params[1]);
                let az = f32x8::splat(inst.params[2]);
                let bx = f32x8::splat(inst.params[3]);
                let by = f32x8::splat(inst.params[4]);
                let bz = f32x8::splat(inst.params[5]);
                let radius = f32x8::splat(inst.get_capsule_radius());

                // pa = p - a, ba = b - a
                let pax = p.x - ax;
                let pay = p.y - ay;
                let paz = p.z - az;
                let bax = bx - ax;
                let bay = by - ay;
                let ba_z = bz - az;

                // h = clamp(dot(pa, ba) / dot(ba, ba), 0, 1)
                let dot_pa_ba = pax * bax + pay * bay + paz * ba_z;
                let dot_ba_ba = bax * bax + bay * bay + ba_z * ba_z;
                // Branchless zero guard for degenerate capsule (a == b)
                let safe_dot = dot_ba_ba.max(f32x8::splat(1e-10));
                let h = (dot_pa_ba / safe_dot).max(f32x8::ZERO).min(f32x8::ONE);

                // length(pa - ba * h) - radius
                let dx = pax - bax * h;
                let dy = pay - bay * h;
                let dz = paz - ba_z * h;
                let d = (dx * dx + dy * dy + dz * dz).sqrt() - radius;

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Cone => {
                let radius = f32x8::splat(inst.params[0]);
                let h = f32x8::splat(inst.params[1]);

                let q_x = (p.x * p.x + p.z * p.z).sqrt();
                let q_y = p.y;

                let k2x = -radius;
                let k2y = h + h;

                // ca_r = q_y < 0 ? radius : 0
                let neg_mask = q_y.cmp_lt(f32x8::ZERO);
                let ca_r = neg_mask.blend(radius, f32x8::ZERO);

                let ca_x = q_x - q_x.min(ca_r);
                let ca_y = q_y.abs() - h;

                let diff_x = -q_x;
                let diff_y = h - q_y;
                let k2_dot = k2x * k2x + k2y * k2y;
                let t = ((diff_x * k2x + diff_y * k2y) / k2_dot)
                    .max(f32x8::ZERO)
                    .min(f32x8::ONE);

                let cb_x = q_x + k2x * t;
                let cb_y = q_y - h + k2y * t;

                // s = (cb_x < 0 && ca_y < 0) ? -1 : 1
                let both_neg = cb_x.cmp_lt(f32x8::ZERO) & ca_y.cmp_lt(f32x8::ZERO);
                let s = both_neg.blend(f32x8::splat(-1.0), f32x8::ONE);

                let d2 = (ca_x * ca_x + ca_y * ca_y).min(cb_x * cb_x + cb_y * cb_y);
                let d = s * d2.sqrt();

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Ellipsoid => {
                // Division Exorcism: 3 divisions instead of 6 (inv_r² = inv_r * inv_r)
                let rx = inst.params[0].max(1e-10);
                let ry = inst.params[1].max(1e-10);
                let rz = inst.params[2].max(1e-10);
                let inv_rx = f32x8::splat(1.0 / rx);
                let inv_ry = f32x8::splat(1.0 / ry);
                let inv_rz = f32x8::splat(1.0 / rz);
                let inv_rx2 = inv_rx * inv_rx;
                let inv_ry2 = inv_ry * inv_ry;
                let inv_rz2 = inv_rz * inv_rz;

                // k0 = length(p * inv_radii)
                let px_r = p.x * inv_rx;
                let py_r = p.y * inv_ry;
                let pz_r = p.z * inv_rz;
                let k0 = (px_r * px_r + py_r * py_r + pz_r * pz_r).sqrt();

                // k1 = length(p * inv_radii²)
                let px_rr = p.x * inv_rx2;
                let py_rr = p.y * inv_ry2;
                let pz_rr = p.z * inv_rz2;
                let k1 = (px_rr * px_rr + py_rr * py_rr + pz_rr * pz_rr).sqrt();

                let d = k0 * (k0 - f32x8::ONE) / k1.max(f32x8::splat(1e-10));
                // sdf_ellipsoid law: at the centre (k1 ≈ 0) return -min(radii)
                let centre = f32x8::splat(-rx.min(ry).min(rz));
                let d = k1.cmp_lt(f32x8::splat(1e-10)).blend(centre, d);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::RoundedCone => {
                let r1 = f32x8::splat(inst.params[0]);
                let r2 = f32x8::splat(inst.params[1]);
                let half_height = f32x8::splat(inst.params[2]);

                let h = half_height + half_height;
                let q_x = (p.x * p.x + p.z * p.z).sqrt();
                let q_y = p.y + half_height;

                let b_val = (r1 - r2) / h;
                let a_val = (f32x8::ONE - b_val * b_val).sqrt();
                let k = q_x * (-b_val) + q_y * a_val;

                // Case 1: k < 0 → bottom sphere
                let d_bottom = (q_x * q_x + q_y * q_y).sqrt() - r1;
                // Case 2: k > a*h → top sphere
                let dy_top = q_y - h;
                let d_top = (q_x * q_x + dy_top * dy_top).sqrt() - r2;
                // Case 3: mantle
                let d_mantle = q_x * a_val + q_y * b_val - r1;

                let mask_bottom = k.cmp_lt(f32x8::ZERO);
                let mask_top = k.cmp_gt(a_val * h);
                let d = mask_bottom.blend(d_bottom, mask_top.blend(d_top, d_mantle));

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Pyramid => {
                let half_height = f32x8::splat(inst.params[0]);
                let h = half_height + half_height;
                let m2 = h * h + f32x8::splat(0.25);
                let half = f32x8::splat(0.5);

                let py = p.y + half_height;
                let abs_px = p.x.abs();
                let abs_pz = p.z.abs();

                // if pz > px { swap(px, pz) }
                let swap_mask = abs_pz.cmp_gt(abs_px);
                let px = swap_mask.blend(abs_pz, abs_px) - half;
                let pz = swap_mask.blend(abs_px, abs_pz) - half;

                let qx = pz;
                let qy = h * py - half * px;
                let qz = h * px + half * py;

                let s = (-qx).max(f32x8::ZERO);
                let t = ((qy - half * pz) / (m2 + f32x8::splat(0.25)))
                    .max(f32x8::ZERO)
                    .min(f32x8::ONE);

                let a = m2 * (qx + s) * (qx + s) + qy * qy;
                let half_t = half * t;
                let b = m2 * (qx + half_t) * (qx + half_t) + (qy - m2 * t) * (qy - m2 * t);

                // d2 = if qy.min(-qx * m2 - qy * 0.5) > 0 { 0 } else { min(a, b) }
                let inner = (-qx * m2 - qy * half).min(qy);
                let zero_mask = inner.cmp_gt(f32x8::ZERO);
                let d2 = zero_mask.blend(f32x8::ZERO, a.min(b));

                // sign = signum(max(qz, -py))
                let sign_input = qz.max(-py);
                let pos = sign_input.cmp_gt(f32x8::ZERO);
                let neg = sign_input.cmp_lt(f32x8::ZERO);
                let sign_val = pos.blend(f32x8::ONE, neg.blend(f32x8::splat(-1.0), f32x8::ZERO));

                let d = ((d2 + qz * qz) / m2).sqrt() * sign_val;

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Octahedron => {
                let s = f32x8::splat(inst.params[0]);
                let three = f32x8::splat(3.0);
                let inv_sqrt3 = f32x8::splat(0.57735027);
                let half = f32x8::splat(0.5);

                let ax = p.x.abs();
                let ay = p.y.abs();
                let az = p.z.abs();
                let m = ax + ay + az - s;

                // d_flat: the "else" branch = m * inv_sqrt3
                let d_flat = m * inv_sqrt3;

                // Exclusive masks for the 3 permutation cases
                let mask1 = (three * ax).cmp_lt(m);
                let mask2 = (three * ay).cmp_lt(m) & !mask1;
                let mask3 = (three * az).cmp_lt(m) & !mask1 & !mask2;
                let mask_any = mask1 | mask2 | mask3;

                // Select q permutation based on which case
                let qx = mask1.blend(ax, mask2.blend(ay, mask3.blend(az, f32x8::ZERO)));
                let qy = mask1.blend(ay, mask2.blend(az, mask3.blend(ax, f32x8::ZERO)));
                let qz = mask1.blend(az, mask2.blend(ax, mask3.blend(ay, f32x8::ZERO)));

                let k = (half * (qz - qy + s)).max(f32x8::ZERO).min(s);
                let vx = qx;
                let vy = qy - s + k;
                let vz = qz - k;
                let d_edge = (vx * vx + vy * vy + vz * vz).sqrt();

                let d = mask_any.blend(d_edge, d_flat);

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::HexPrism => {
                let hex_radius = f32x8::splat(inst.params[0]);
                let half_height = f32x8::splat(inst.params[1]);

                let kx = f32x8::splat(-0.8660254_f32);
                let ky = f32x8::splat(0.5_f32);
                let kz_c = f32x8::splat(0.57735027_f32);
                let two = f32x8::splat(2.0);

                let mut px = p.x.abs();
                let mut py = p.y.abs();
                let pz = p.z.abs();

                // Reflect across hex symmetry
                let dot_kxy = kx * px + ky * py;
                let reflect = two * dot_kxy.min(f32x8::ZERO);
                px -= reflect * kx;
                py -= reflect * ky;

                // Clamp and compute XY distance
                let clamped_x = px.max(-kz_c * hex_radius).min(kz_c * hex_radius);
                let dx = px - clamped_x;
                let dy = py - hex_radius;
                let d_xy_len = (dx * dx + dy * dy).sqrt();

                // signum(dy)
                let pos = dy.cmp_gt(f32x8::ZERO);
                let neg = dy.cmp_lt(f32x8::ZERO);
                let dy_sign = pos.blend(f32x8::ONE, neg.blend(f32x8::splat(-1.0), f32x8::ZERO));
                let d_xy = d_xy_len * dy_sign;

                let d_z = pz - half_height;

                // max(d_xy, d_z).min(0) + sqrt(max(d_xy,0)^2 + max(d_z,0)^2)
                let d_xy_pos = d_xy.max(f32x8::ZERO);
                let d_z_pos = d_z.max(f32x8::ZERO);
                let d = d_xy.max(d_z).min(f32x8::ZERO)
                    + (d_xy_pos * d_xy_pos + d_z_pos * d_z_pos).sqrt();

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Link => {
                let half_length = f32x8::splat(inst.params[0]);
                let r1 = f32x8::splat(inst.params[1]);
                let r2 = f32x8::splat(inst.params[2]);

                let qx = p.x;
                let qy = (p.y.abs() - half_length).max(f32x8::ZERO);
                let qz = p.z;

                let xy_len = (qx * qx + qy * qy).sqrt() - r1;
                let d = (xy_len * xy_len + qz * qz).sqrt() - r2;

                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            // === Extended Primitives (native SIMD) ===
            OpCode::RoundedBox => {
                // q = abs(p) - half_extents; max(q, 0).length() + min(max(q.x, q.y, q.z), 0) - r
                let hx = f32x8::splat(inst.params[0]);
                let hy = f32x8::splat(inst.params[1]);
                let hz = f32x8::splat(inst.params[2]);
                let rr = f32x8::splat(inst.params[3]);
                let qx = p.x.abs() - hx;
                let qy = p.y.abs() - hy;
                let qz = p.z.abs() - hz;
                let qx_pos = qx.max(f32x8::ZERO);
                let qy_pos = qy.max(f32x8::ZERO);
                let qz_pos = qz.max(f32x8::ZERO);
                let outer = (qx_pos * qx_pos + qy_pos * qy_pos + qz_pos * qz_pos).sqrt();
                let inner = qx.max(qy).max(qz).min(f32x8::ZERO);
                let d = outer + inner - rr;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::CappedCone => {
                // 2D profile in (length(p.xz), p.y) space
                let h = f32x8::splat(inst.params[0]);
                let r1 = f32x8::splat(inst.params[1]);
                let r2 = f32x8::splat(inst.params[2]);
                let qx = (p.x * p.x + p.z * p.z).sqrt();
                let qy = p.y;
                // k2 = (r2 - r1, 2*h)
                let k2x = r2 - r1;
                let k2y = h * f32x8::splat(2.0);
                let k2_dot = k2x * k2x + k2y * k2y;
                // ca = (qx - min(qx, if qy<0 {r1} else {r2}), abs(qy) - h)
                let neg_mask = qy.cmp_lt(f32x8::ZERO);
                let min_r = neg_mask.blend(r1, r2);
                let ca_x = qx - qx.min(min_r);
                let ca_y = qy.abs() - h;
                // t = clamp(dot(k1-q, k2) / dot(k2,k2), 0, 1); k1=(r2,h)
                let d_to_k1_x = r2 - qx;
                let d_to_k1_y = h - qy;
                let num = d_to_k1_x * k2x + d_to_k1_y * k2y;
                let safe_k2_dot = k2_dot.max(f32x8::splat(0.0001));
                let t = (num / safe_k2_dot).max(f32x8::ZERO).min(f32x8::ONE);
                // cb = q - k1 + k2*t
                let cb_x = qx - r2 + k2x * t;
                let cb_y = qy - h + k2y * t;
                let ca_d2 = ca_x * ca_x + ca_y * ca_y;
                let cb_d2 = cb_x * cb_x + cb_y * cb_y;
                // s = -1 if cb.x<0 && ca.y<0, else 1
                let both_neg = cb_x.cmp_lt(f32x8::ZERO) & ca_y.cmp_lt(f32x8::ZERO);
                let s = both_neg.blend(f32x8::splat(-1.0), f32x8::ONE);
                let d = s * ca_d2.min(cb_d2).sqrt();
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::CappedTorus => {
                let major_r = f32x8::splat(inst.params[0]);
                let minor_r = f32x8::splat(inst.params[1]);
                let sc_sin = f32x8::splat(inst.params[2].sin());
                let sc_cos = f32x8::splat(inst.params[2].cos());
                let px = p.x.abs();
                // k = sc.cos*px > sc.sin*py ? sc.sin*px + sc.cos*py : sqrt(px² + py²)
                let dot_val = sc_sin * px + sc_cos * p.y;
                let len_val = (px * px + p.y * p.y).sqrt();
                let mask = (sc_cos * px).cmp_gt(sc_sin * p.y);
                let k = mask.blend(dot_val, len_val);
                // sqrt(px² + py² + pz² + R² - 2*R*k) - r
                let inner = px * px + p.y * p.y + p.z * p.z + major_r * major_r
                    - f32x8::splat(2.0) * major_r * k;
                let d = inner.max(f32x8::ZERO).sqrt() - minor_r;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::RoundedCylinder => {
                let radius = f32x8::splat(inst.params[0]);
                let round_r = f32x8::splat(inst.params[1]);
                let half_h = f32x8::splat(inst.params[2]);
                let dx = (p.x * p.x + p.z * p.z).sqrt() - f32x8::splat(2.0) * radius + round_r;
                let dy = p.y.abs() - half_h;
                let dx_pos = dx.max(f32x8::ZERO);
                let dy_pos = dy.max(f32x8::ZERO);
                let d = dx.max(dy).min(f32x8::ZERO) + (dx_pos * dx_pos + dy_pos * dy_pos).sqrt()
                    - round_r;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::TriangularPrism => {
                let width = f32x8::splat(inst.params[0]);
                let half_depth = f32x8::splat(inst.params[1]);
                let qx = p.x.abs();
                let qy = p.y; // not abs for y
                let qz = p.z.abs();
                // 0.866025 = sqrt(3)/2
                let sqrt3_half = f32x8::splat(0.866025);
                let half = f32x8::splat(0.5);
                let d =
                    (qz - half_depth).max((qx * sqrt3_half + qy * half).max(-qy) - width * half);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::CutSphere => {
                let radius = f32x8::splat(inst.params[0]);
                let ch = f32x8::splat(inst.params[1]); // cut_height
                let w = (radius * radius - ch * ch).max(f32x8::ZERO).sqrt();
                let qx = (p.x * p.x + p.z * p.z).sqrt();
                let qy = p.y;
                let q_len = (qx * qx + qy * qy).sqrt();
                // Three regions via branchless blend
                let s1 = (ch - radius) * qx * qx + w * w * (ch + radius - f32x8::splat(2.0) * qy);
                let s2 = ch * qx - w * qy;
                let s = s1.max(s2);
                // d_sphere = length(q) - r
                let d_sphere = q_len - radius;
                // d_plane = h - q.y
                let d_plane = ch - qy;
                // d_edge = length(q - (w, h))
                let ex = qx - w;
                let ey = qy - ch;
                let d_edge = (ex * ex + ey * ey).sqrt();
                // if s < 0 -> d_sphere; elif qx < w -> d_plane; else -> d_edge
                let mask_s_neg = s.cmp_lt(f32x8::ZERO);
                let mask_qx_lt_w = qx.cmp_lt(w);
                let d = mask_s_neg.blend(d_sphere, mask_qx_lt_w.blend(d_plane, d_edge));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::CutHollowSphere => {
                let radius = f32x8::splat(inst.params[0]);
                let ch = f32x8::splat(inst.params[1]);
                let thickness = f32x8::splat(inst.params[2]);
                let w = (radius * radius - ch * ch).max(f32x8::ZERO).sqrt();
                let qx = (p.x * p.x + p.z * p.z).sqrt();
                let qy = p.y;
                // if h*qx < w*qy -> length(q - (w,h)) - t; else -> abs(length(q) - r) - t
                let mask = (ch * qx).cmp_lt(w * qy);
                let ex = qx - w;
                let ey = qy - ch;
                let d_cap = (ex * ex + ey * ey).sqrt() - thickness;
                let q_len = (qx * qx + qy * qy).sqrt();
                let d_shell = (q_len - radius).abs() - thickness;
                let d = mask.blend(d_cap, d_shell);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::DeathStar => {
                let ra = f32x8::splat(inst.params[0]);
                let rb = f32x8::splat(inst.params[1]);
                let dd = f32x8::splat(inst.params[2]);
                let two = f32x8::splat(2.0);
                // a = (ra² - rb² + d²) / (2d)
                #[allow(clippy::suspicious_operation_groupings)]
                let a = (ra * ra - rb * rb + dd * dd) / (two * dd);
                let b = (ra * ra - a * a).max(f32x8::ZERO).sqrt();
                let p2x = p.x;
                let p2y = (p.y * p.y + p.z * p.z).sqrt();
                // Condition: p2.x*b - p2.y*a > d*max(b - p2.y, 0)
                let lhs = p2x * b - p2y * a;
                let rhs = dd * (b - p2y).max(f32x8::ZERO);
                let mask = lhs.cmp_gt(rhs);
                // d_edge = length(p2 - (a, b))
                let ex = p2x - a;
                let ey = p2y - b;
                let d_edge = (ex * ex + ey * ey).sqrt();
                // d_main = max(length(p2) - ra, -(length(p2 - (d,0)) - rb))
                let p2_len = (p2x * p2x + p2y * p2y).sqrt();
                let dx = p2x - dd;
                let d_ra = p2_len - ra;
                let d_rb = -((dx * dx + p2y * p2y).sqrt() - rb);
                let d_main = d_ra.max(d_rb);
                let d = mask.blend(d_edge, d_main);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::SolidAngle => {
                let c_sin = f32x8::splat(inst.params[0].sin());
                let c_cos = f32x8::splat(inst.params[0].cos());
                let radius = f32x8::splat(inst.params[1]);
                let qx = (p.x * p.x + p.z * p.z).sqrt();
                let qy = p.y;
                let q_len = (qx * qx + qy * qy).sqrt();
                let l = q_len - radius;
                // dot(q, c) clamped to [0, radius]
                let q_dot_c = (qx * c_sin + qy * c_cos).max(f32x8::ZERO).min(radius);
                // m = length(q - c * clamp(dot(q,c), 0, r))
                let proj_x = qx - c_sin * q_dot_c;
                let proj_y = qy - c_cos * q_dot_c;
                let m = (proj_x * proj_x + proj_y * proj_y).sqrt();
                // sign = c.y*q.x - c.x*q.y < 0 ? -1 : 1
                let sign_val = c_cos * qx - c_sin * qy;
                let neg_mask = sign_val.cmp_lt(f32x8::ZERO);
                let sign = neg_mask.blend(f32x8::splat(-1.0), f32x8::ONE);
                let d = l.max(m * sign);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Rhombus => {
                let la = f32x8::splat(inst.params[0]);
                let lb = f32x8::splat(inst.params[1]);
                let half_h = f32x8::splat(inst.params[2]);
                let rr = f32x8::splat(inst.params[3]);
                let ax = p.x.abs();
                let ay = p.y.abs();
                let az = p.z.abs();
                // ndot(b, b - 2*(px,pz)) = la*(la-2*px) - lb*(lb-2*pz)
                //                        = la² - 2*la*px - lb² + 2*lb*pz
                let b_dot_b = la * la + lb * lb;
                let ndot_val =
                    la * (la - f32x8::splat(2.0) * ax) - lb * (lb - f32x8::splat(2.0) * az);
                let f = (ndot_val / b_dot_b).max(f32x8::splat(-1.0)).min(f32x8::ONE);
                // q_xz = length((px,pz) - 0.5*b*(1-f, 1+f))
                let half = f32x8::splat(0.5);
                let proj_x = ax - half * la * (f32x8::ONE - f);
                let proj_z = az - half * lb * (f32x8::ONE + f);
                let qxz_len = (proj_x * proj_x + proj_z * proj_z).sqrt();
                // sign(px*lb + pz*la - la*lb)
                let sign_input = ax * lb + az * la - la * lb;
                let pos = sign_input.cmp_gt(f32x8::ZERO);
                let neg = sign_input.cmp_lt(f32x8::ZERO);
                let sign = pos.blend(f32x8::ONE, neg.blend(f32x8::splat(-1.0), f32x8::ZERO));
                let dx = qxz_len * sign - rr;
                let dy = ay - half_h;
                let dx_pos = dx.max(f32x8::ZERO);
                let dy_pos = dy.max(f32x8::ZERO);
                let d = dx.max(dy).min(f32x8::ZERO) + (dx_pos * dx_pos + dy_pos * dy_pos).sqrt();
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Horseshoe => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_horseshoe(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                        inst.params[4],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Vesica => {
                let d = eval_per_lane(&p, |pt| sdf_vesica(pt, inst.params[0], inst.params[1]));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::InfiniteCylinder => {
                // Simple SIMD: length(p.xz) - r
                let r = f32x8::splat(inst.params[0]);
                let d = (p.x * p.x + p.z * p.z).sqrt() - r;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::InfiniteCone => {
                let d = eval_per_lane(&p, |pt| sdf_infinite_cone(pt, inst.params[0]));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Gyroid => {
                // ★ Native SIMD: sin(x)*cos(y) + sin(y)*cos(z) + sin(z)*cos(x)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let spx = p.x * scale;
                let spy = p.y * scale;
                let spz = p.z * scale;
                let sx = sin_approx(spx);
                let cx = cos_approx(spx);
                let sy = sin_approx(spy);
                let cy = cos_approx(spy);
                let sz = sin_approx(spz);
                let cz = cos_approx(spz);
                let d = (sx * cy + sy * cz + sz * cx).abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Heart => {
                let d = eval_per_lane(&p, |pt| sdf_heart(pt, inst.params[0]));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Tube => {
                // SIMD: hollow cylinder
                let outer_r = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let half_h = f32x8::splat(inst.params[2]);
                let xz_len = (p.x * p.x + p.z * p.z).sqrt();
                let dx = (xz_len - outer_r).abs() - thickness;
                let dy = p.y.abs() - half_h;
                let dx_pos = dx.max(f32x8::ZERO);
                let dy_pos = dy.max(f32x8::ZERO);
                let d = (dx_pos * dx_pos + dy_pos * dy_pos).sqrt() + dx.max(dy).min(f32x8::ZERO);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Barrel => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_barrel(pt, inst.params[0], inst.params[1], inst.params[2])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Diamond => {
                let d = eval_per_lane(&p, |pt| sdf_diamond(pt, inst.params[0], inst.params[1]));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::ChamferedCube => {
                let d = eval_per_lane(&p, |pt| {
                    let he = Vec3::new(inst.params[0], inst.params[1], inst.params[2]);
                    sdf_chamfered_cube(pt, he, inst.params[3])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::SchwarzP => {
                // ★ Native SIMD: cos(x) + cos(y) + cos(z)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let d =
                    (cos_approx(p.x * scale) + cos_approx(p.y * scale) + cos_approx(p.z * scale))
                        .abs()
                        * inv_scale
                        - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Superellipsoid => {
                let d = eval_per_lane(&p, |pt| {
                    let he = Vec3::new(inst.params[0], inst.params[1], inst.params[2]);
                    sdf_superellipsoid(pt, he, inst.params[3], inst.params[4])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::RoundedX => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_rounded_x(pt, inst.params[0], inst.params[1], inst.params[2])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Pie => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_pie(pt, inst.params[0], inst.params[1], inst.params[2])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Trapezoid => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_trapezoid(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Parallelogram => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_parallelogram(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Tunnel => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_tunnel(pt, inst.params[0], inst.params[1], inst.params[2])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::UnevenCapsule => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_uneven_capsule(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Egg => {
                let d = eval_per_lane(&p, |pt| sdf_egg(pt, inst.params[0], inst.params[1]));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::ArcShape => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_arc_shape(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Moon => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_moon(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::CrossShape => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_cross_shape(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::BlobbyCross => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_blobby_cross(pt, inst.params[0], inst.params[1])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::ParabolaSegment => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_parabola_segment(pt, inst.params[0], inst.params[1], inst.params[2])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::RegularPolygon => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_regular_polygon(pt, inst.params[0], inst.params[1], inst.params[2])
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::StarPolygon => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_star_polygon(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Stairs => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_stairs(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Helix => {
                let d = eval_per_lane(&p, |pt| {
                    sdf_helix(
                        pt,
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    )
                });
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Tetrahedron => {
                // ★ Native SIMD: 4 dot products + max (no abs — tetrahedron normals are signed)
                let radius = f32x8::splat(inst.params[0]);
                let s = f32x8::splat(0.577_350_26_f32); // 1/sqrt(3)
                let ns = f32x8::splat(-0.577_350_26_f32);
                // n0=(s,s,s) n1=(-s,-s,s) n2=(-s,s,-s) n3=(s,-s,-s)
                let d0 = p.x * s + p.y * s + p.z * s;
                let d1 = p.x * ns + p.y * ns + p.z * s;
                let d2 = p.x * ns + p.y * s + p.z * ns;
                let d3 = p.x * s + p.y * ns + p.z * ns;
                let d = d0.max(d1).max(d2).max(d3) - radius;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Dodecahedron => {
                // ★ Native SIMD: 6 abs-dot products (GDF_DODECAHEDRON normals)
                let radius = f32x8::splat(inst.params[0]);
                let a = f32x8::splat(0.850_650_8_f32); // ICO_B
                let b = f32x8::splat(0.525_731_1_f32); // ICO_A
                                                       // n0=(0,a,b) n1=(0,a,-b) n2=(a,b,0) n3=(-a,b,0) n4=(b,0,a) n5=(b,0,-a)
                let d0 = (p.y * a + p.z * b).abs();
                let d1 = (p.y * a - p.z * b).abs();
                let d2 = (p.x * a + p.y * b).abs();
                let d3 = (p.x * a - p.y * b).abs();
                let d4 = (p.x * b + p.z * a).abs();
                let d5 = (p.x * b - p.z * a).abs();
                let d = d0.max(d1).max(d2).max(d3).max(d4).max(d5) - radius;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Icosahedron => {
                // ★ Native SIMD: 10 abs-dot products (octahedron[4] + icosahedron[6])
                let radius = f32x8::splat(inst.params[0]);
                let s = f32x8::splat(0.577_350_26_f32); // 1/sqrt(3)
                let ia = f32x8::splat(0.525_731_1_f32); // ICO_A
                let ib = f32x8::splat(0.850_650_8_f32); // ICO_B
                                                        // Octahedron normals (4): abs(dot) with (±s,±s,±s) variants
                let d0 = (p.x * s + p.y * s + p.z * s).abs();
                let d1 = (-p.x * s + p.y * s + p.z * s).abs();
                let d2 = (p.x * s - p.y * s + p.z * s).abs();
                let d3 = (p.x * s + p.y * s - p.z * s).abs();
                // Icosahedron normals (6): abs(dot) with (0,±a,±b) permutations
                let d4 = (p.y * ia + p.z * ib).abs();
                let d5 = (p.y * ia - p.z * ib).abs();
                let d6 = (p.x * ia + p.y * ib).abs();
                let d7 = (p.x * ia - p.y * ib).abs();
                let d8 = (p.x * ib + p.z * ia).abs();
                let d9 = (p.x * ib - p.z * ia).abs();
                let d = d0
                    .max(d1)
                    .max(d2)
                    .max(d3)
                    .max(d4)
                    .max(d5)
                    .max(d6)
                    .max(d7)
                    .max(d8)
                    .max(d9)
                    - radius;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::TruncatedOctahedron => {
                // ★ Native SIMD: 7 abs-dot products (cube[3] + octahedron[4])
                let radius = f32x8::splat(inst.params[0]);
                let s = f32x8::splat(0.577_350_26_f32);
                // Cube normals (3): abs of each axis
                let d0 = p.x.abs();
                let d1 = p.y.abs();
                let d2 = p.z.abs();
                // Octahedron normals (4)
                let d3 = (p.x * s + p.y * s + p.z * s).abs();
                let d4 = (-p.x * s + p.y * s + p.z * s).abs();
                let d5 = (p.x * s - p.y * s + p.z * s).abs();
                let d6 = (p.x * s + p.y * s - p.z * s).abs();
                let d = d0.max(d1).max(d2).max(d3).max(d4).max(d5).max(d6) - radius;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::TruncatedIcosahedron => {
                // ★ Native SIMD: 16 abs-dot products (oct[4] + ico[6] + dodec[6])
                let radius = f32x8::splat(inst.params[0]);
                let s = f32x8::splat(0.577_350_26_f32);
                let ia = f32x8::splat(0.525_731_1_f32);
                let ib = f32x8::splat(0.850_650_8_f32);
                // Octahedron (4)
                let d0 = (p.x * s + p.y * s + p.z * s).abs();
                let d1 = (-p.x * s + p.y * s + p.z * s).abs();
                let d2 = (p.x * s - p.y * s + p.z * s).abs();
                let d3 = (p.x * s + p.y * s - p.z * s).abs();
                // Icosahedron (6)
                let d4 = (p.y * ia + p.z * ib).abs();
                let d5 = (p.y * ia - p.z * ib).abs();
                let d6 = (p.x * ia + p.y * ib).abs();
                let d7 = (p.x * ia - p.y * ib).abs();
                let d8 = (p.x * ib + p.z * ia).abs();
                let d9 = (p.x * ib - p.z * ia).abs();
                // Dodecahedron (6)
                let d10 = (p.y * ib + p.z * ia).abs();
                let d11 = (p.y * ib - p.z * ia).abs();
                let d12 = (p.x * ib + p.y * ia).abs();
                let d13 = (p.x * ib - p.y * ia).abs();
                let d14 = (p.x * ia + p.z * ib).abs();
                let d15 = (p.x * ia - p.z * ib).abs();
                let d = d0
                    .max(d1)
                    .max(d2)
                    .max(d3)
                    .max(d4)
                    .max(d5)
                    .max(d6)
                    .max(d7)
                    .max(d8)
                    .max(d9)
                    .max(d10)
                    .max(d11)
                    .max(d12)
                    .max(d13)
                    .max(d14)
                    .max(d15)
                    - radius;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::BoxFrame => {
                // ★ Native SIMD: abs, max, min, sqrt — all native ops
                let bx = f32x8::splat(inst.params[0]);
                let by = f32x8::splat(inst.params[1]);
                let bz = f32x8::splat(inst.params[2]);
                let e = f32x8::splat(inst.params[3]);
                let zero = f32x8::ZERO;
                // p = abs(p) - half_extents
                let px = p.x.abs() - bx;
                let py = p.y.abs() - by;
                let pz = p.z.abs() - bz;
                // q = abs(p + e) - e
                let qx = (px + e).abs() - e;
                let qy = (py + e).abs() - e;
                let qz = (pz + e).abs() - e;
                // d1 = length(max(vec3(px,qy,qz),0)) + min(max(px,max(qy,qz)),0)
                let v1x = px.max(zero);
                let v1y = qy.max(zero);
                let v1z = qz.max(zero);
                let d1 = (v1x * v1x + v1y * v1y + v1z * v1z).sqrt() + px.max(qy.max(qz)).min(zero);
                // d2
                let v2x = qx.max(zero);
                let v2y = py.max(zero);
                let v2z = qz.max(zero);
                let d2 = (v2x * v2x + v2y * v2y + v2z * v2z).sqrt() + qx.max(py.max(qz)).min(zero);
                // d3
                let v3x = qx.max(zero);
                let v3y = qy.max(zero);
                let v3z = pz.max(zero);
                let d3 = (v3x * v3x + v3y * v3y + v3z * v3z).sqrt() + qx.max(qy.max(pz)).min(zero);
                let d = d1.min(d2).min(d3);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::DiamondSurface => {
                // ★ Native SIMD: sin(x)*sin(y)*sin(z) + sin(x)*cos(y)*cos(z) + cos(x)*sin(y)*cos(z) + cos(x)*cos(y)*sin(z)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let spx = p.x * scale;
                let spy = p.y * scale;
                let spz = p.z * scale;
                let sx = sin_approx(spx);
                let cx = cos_approx(spx);
                let sy = sin_approx(spy);
                let cy = cos_approx(spy);
                let sz = sin_approx(spz);
                let cz = cos_approx(spz);
                let d = sx * sy * sz + sx * cy * cz + cx * sy * cz + cx * cy * sz;
                let d = d.abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Neovius => {
                // ★ Native SIMD: 3*(cos(x)+cos(y)+cos(z)) + 4*cos(x)*cos(y)*cos(z)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let three = f32x8::splat(3.0);
                let four = f32x8::splat(4.0);
                let cx = cos_approx(p.x * scale);
                let cy = cos_approx(p.y * scale);
                let cz = cos_approx(p.z * scale);
                let d = three * (cx + cy + cz) + four * cx * cy * cz;
                let d = d.abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Lidinoid => {
                // ★ Native SIMD + double-angle identities (0 extra trig calls)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let two = f32x8::splat(2.0);
                let one = f32x8::ONE;
                let half = f32x8::splat(0.5);
                let spx = p.x * scale;
                let spy = p.y * scale;
                let spz = p.z * scale;
                let sx = sin_approx(spx);
                let cx = cos_approx(spx);
                let sy = sin_approx(spy);
                let cy = cos_approx(spy);
                let sz = sin_approx(spz);
                let cz = cos_approx(spz);
                // sin(2x) = 2*sx*cx, cos(2x) = 2*cx*cx - 1
                let s2x = two * sx * cx;
                let s2y = two * sy * cy;
                let s2z = two * sz * cz;
                let c2x = two * cx * cx - one;
                let c2y = two * cy * cy - one;
                let c2z = two * cz * cz - one;
                let term1 = half * (s2x * cy * sz + sx * s2y * cz + cx * sy * s2z);
                let term2 = half * (c2x * c2y + c2y * c2z + c2z * c2x);
                let d = (term1 - term2 + f32x8::splat(0.15)).abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::IWP => {
                // ★ Native SIMD + double-angle: cos(2x) = 2*cos²(x)-1
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let two = f32x8::splat(2.0);
                let one = f32x8::ONE;
                let cx = cos_approx(p.x * scale);
                let cy = cos_approx(p.y * scale);
                let cz = cos_approx(p.z * scale);
                let c2x = two * cx * cx - one;
                let c2y = two * cy * cy - one;
                let c2z = two * cz * cz - one;
                let d = two * (cx * cy + cy * cz + cz * cx) - (c2x + c2y + c2z);
                let d = d.abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::FRD => {
                // ★ Native SIMD + double-angle: cos(2x) = 2*cos²(x)-1
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let two = f32x8::splat(2.0);
                let one = f32x8::ONE;
                let spx = p.x * scale;
                let spy = p.y * scale;
                let spz = p.z * scale;
                let sx = sin_approx(spx);
                let cx = cos_approx(spx);
                let sy = sin_approx(spy);
                let cy = cos_approx(spy);
                let sz = sin_approx(spz);
                let cz = cos_approx(spz);
                let c2x = two * cx * cx - one;
                let c2y = two * cy * cy - one;
                let c2z = two * cz * cz - one;
                let d = c2x * sy * cz + cx * c2y * sz + sx * cy * c2z;
                let d = d.abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::FischerKochS => {
                // ★ Native SIMD + double-angle (same as FRD with -0.4 offset)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let two = f32x8::splat(2.0);
                let one = f32x8::ONE;
                let spx = p.x * scale;
                let spy = p.y * scale;
                let spz = p.z * scale;
                let sx = sin_approx(spx);
                let cx = cos_approx(spx);
                let sy = sin_approx(spy);
                let cy = cos_approx(spy);
                let sz = sin_approx(spz);
                let cz = cos_approx(spz);
                let c2x = two * cx * cx - one;
                let c2y = two * cy * cy - one;
                let c2z = two * cz * cz - one;
                let d = c2x * sy * cz + cx * c2y * sz + sx * cy * c2z - f32x8::splat(0.4);
                let d = d.abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::PMY => {
                // ★ Native SIMD + double-angle: sin(2x) = 2*sin(x)*cos(x)
                let scale = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let inv_scale = f32x8::splat(1.0 / inst.params[0]);
                let two = f32x8::splat(2.0);
                let spx = p.x * scale;
                let spy = p.y * scale;
                let spz = p.z * scale;
                let sx = sin_approx(spx);
                let cx = cos_approx(spx);
                let sy = sin_approx(spy);
                let cy = cos_approx(spy);
                let sz = sin_approx(spz);
                let cz = cos_approx(spz);
                let s2x = two * sx * cx;
                let s2y = two * sy * cy;
                let s2z = two * sz * cz;
                let d = two * cx * cy * cz + s2x * sy + sx * s2z + s2y * sz;
                let d = d.abs() * inv_scale - thickness;
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            // === Binary Operations ===
            OpCode::Union => {
                vsp -= 1;
                let b = value_stack[vsp];
                value_stack[vsp - 1] = value_stack[vsp - 1].min(b);
            }

            OpCode::Intersection => {
                vsp -= 1;
                let b = value_stack[vsp];
                value_stack[vsp - 1] = value_stack[vsp - 1].max(b);
            }

            OpCode::Subtraction => {
                vsp -= 1;
                let b = value_stack[vsp];
                value_stack[vsp - 1] = value_stack[vsp - 1].max(-b);
            }

            OpCode::SmoothUnion => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let k = f32x8::splat(inst.params[0]);
                let rk = f32x8::splat(inst.params[1]); // Division Exorcism: precomputed 1/k
                value_stack[vsp - 1] = smooth_min_simd_rk(a, b, k, rk);
            }

            OpCode::SmoothIntersection => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let k = f32x8::splat(inst.params[0]);
                let rk = f32x8::splat(inst.params[1]);
                value_stack[vsp - 1] = -smooth_min_simd_rk(-a, -b, k, rk);
            }

            OpCode::SmoothSubtraction => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let k = f32x8::splat(inst.params[0]);
                let rk = f32x8::splat(inst.params[1]);
                value_stack[vsp - 1] = -smooth_min_simd_rk(-a, b, k, rk);
            }

            OpCode::ChamferUnion => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = f32x8::splat(inst.params[0]);
                value_stack[vsp - 1] = chamfer_min_simd(a, b, r);
            }

            OpCode::ChamferIntersection => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = f32x8::splat(inst.params[0]);
                value_stack[vsp - 1] = -chamfer_min_simd(-a, -b, r);
            }

            OpCode::ChamferSubtraction => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = f32x8::splat(inst.params[0]);
                value_stack[vsp - 1] = -chamfer_min_simd(-a, b, r);
            }

            OpCode::StairsUnion => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                value_stack[vsp - 1] = stairs_min_simd(a, b, inst.params[0], inst.params[1]);
            }

            OpCode::StairsIntersection => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                value_stack[vsp - 1] = -stairs_min_simd(-a, -b, inst.params[0], inst.params[1]);
            }

            OpCode::StairsSubtraction => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                value_stack[vsp - 1] = -stairs_min_simd(-a, b, inst.params[0], inst.params[1]);
            }

            OpCode::XOR => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                value_stack[vsp - 1] = a.min(b).max(-a.max(b));
            }

            OpCode::Morph => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let t = f32x8::splat(inst.params[0]);
                let one_minus_t = f32x8::splat(1.0 - inst.params[0]);
                value_stack[vsp - 1] = a * one_minus_t + b * t;
            }

            OpCode::ColumnsUnion => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = inst.params[0];
                let n = inst.params[1];
                let d = eval_per_lane_binary(a, b, |av, bv| sdf_columns_union(av, bv, r, n));
                value_stack[vsp - 1] = d;
            }

            OpCode::ColumnsIntersection => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = inst.params[0];
                let n = inst.params[1];
                let d = eval_per_lane_binary(a, b, |av, bv| sdf_columns_intersection(av, bv, r, n));
                value_stack[vsp - 1] = d;
            }

            OpCode::ColumnsSubtraction => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = inst.params[0];
                let n = inst.params[1];
                let d = eval_per_lane_binary(a, b, |av, bv| sdf_columns_subtraction(av, bv, r, n));
                value_stack[vsp - 1] = d;
            }

            OpCode::Pipe => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = f32x8::splat(inst.params[0]);
                value_stack[vsp - 1] = (a * a + b * b).sqrt() - r;
            }

            OpCode::Engrave => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let r = f32x8::splat(inst.params[0]);
                let s = f32x8::splat(std::f32::consts::FRAC_1_SQRT_2);
                let abs_b = b.abs();
                value_stack[vsp - 1] = a.max((a + r - abs_b) * s);
            }

            OpCode::Groove => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let ra = f32x8::splat(inst.params[0]);
                let rb = f32x8::splat(inst.params[1]);
                let abs_b = b.abs();
                value_stack[vsp - 1] = a.max((a + ra).min(rb - abs_b));
            }

            OpCode::Tongue => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let ra = f32x8::splat(inst.params[0]);
                let rb = f32x8::splat(inst.params[1]);
                let abs_b = b.abs();
                value_stack[vsp - 1] = a.min((a - ra).max(abs_b - rb));
            }

            // === Transforms ===
            OpCode::Translate => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Translate,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let offset = Vec3x8 {
                    x: f32x8::splat(inst.params[0]),
                    y: f32x8::splat(inst.params[1]),
                    z: f32x8::splat(inst.params[2]),
                };
                p = p - offset;
            }

            OpCode::Rotate => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Rotate,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let q = glam::Quat::from_xyzw(
                    inst.params[0],
                    inst.params[1],
                    inst.params[2],
                    inst.params[3],
                );
                let qx8 = Quatx8::splat(q.inverse());
                p = qx8.mul_vec3(p);
            }

            OpCode::Scale => {
                let inv_factor = inst.params[0]; // precomputed 1.0/factor
                let factor = inst.params[1]; // original factor
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Scale,
                    params: [factor, 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // Multiply by precomputed inverse (no division)
                p = p * f32x8::splat(inv_factor);
                scale_correction *= f32x8::splat(factor);
            }

            OpCode::ScaleNonUniform => {
                let inv_sx = inst.params[0]; // precomputed 1.0/sx
                let inv_sy = inst.params[1]; // precomputed 1.0/sy
                let inv_sz = inst.params[2]; // precomputed 1.0/sz
                let min_factor = inst.params[3]; // precomputed min(sx,sy,sz)

                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::ScaleNonUniform,
                    params: [min_factor, 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // Multiply by precomputed inverses (no division)
                p = Vec3x8 {
                    x: p.x * f32x8::splat(inv_sx),
                    y: p.y * f32x8::splat(inv_sy),
                    z: p.z * f32x8::splat(inv_sz),
                };
                scale_correction *= f32x8::splat(min_factor);
            }

            // === Modifiers ===
            OpCode::Twist => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Twist,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let k = f32x8::splat(inst.params[0]);
                let angle = k * p.y;
                let c = cos_approx(angle);
                let s = sin_approx(angle);
                let new_x = c * p.x - s * p.z;
                let new_z = s * p.x + c * p.z;
                p = Vec3x8 {
                    x: new_x,
                    y: p.y,
                    z: new_z,
                };
            }

            OpCode::Bend => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Bend,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let k = f32x8::splat(inst.params[0]);
                let angle = k * p.x;
                let c = cos_approx(angle);
                let s = sin_approx(angle);
                let new_x = c * p.x - s * p.y;
                let new_y = s * p.x + c * p.y;
                p = Vec3x8 {
                    x: new_x,
                    y: new_y,
                    z: p.z,
                };
            }

            OpCode::RepeatInfinite => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::RepeatInfinite,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let sx = f32x8::splat(inst.params[0]);
                let sy = f32x8::splat(inst.params[1]);
                let sz = f32x8::splat(inst.params[2]);
                // Division Exorcism: precomputed reciprocal spacing
                let rsx = f32x8::splat(inst.params[3]);
                let rsy = f32x8::splat(inst.params[4]);
                let rsz = f32x8::splat(inst.params[5]);
                let half = f32x8::splat(0.5);

                // modulo operation: p - spacing * floor(p * recip_spacing + 0.5)
                p = Vec3x8 {
                    x: p.x - sx * (p.x * rsx + half).floor(),
                    y: p.y - sy * (p.y * rsy + half).floor(),
                    z: p.z - sz * (p.z * rsz + half).floor(),
                };
            }

            OpCode::RepeatFinite => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::RepeatFinite,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // modifier_repeat_finite law: limit = count * 0.5
                let cx = f32x8::splat(inst.params[0] * 0.5);
                let cy = f32x8::splat(inst.params[1] * 0.5);
                let cz = f32x8::splat(inst.params[2] * 0.5);
                let sx = f32x8::splat(inst.params[3]);
                let sy = f32x8::splat(inst.params[4]);
                let sz = f32x8::splat(inst.params[5]);

                // Division Exorcism: precompute scalar reciprocals, splat to SIMD
                let rsx = f32x8::splat(1.0 / inst.params[3]);
                let rsy = f32x8::splat(1.0 / inst.params[4]);
                let rsz = f32x8::splat(1.0 / inst.params[5]);

                // clamp(round(p * recip_spacing), -count, count) * spacing
                let ix = ((p.x * rsx).round()).max(-cx).min(cx);
                let iy = ((p.y * rsy).round()).max(-cy).min(cy);
                let iz = ((p.z * rsz).round()).max(-cz).min(cz);

                p = Vec3x8 {
                    x: p.x - ix * sx,
                    y: p.y - iy * sy,
                    z: p.z - iz * sz,
                };
            }

            OpCode::Round => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Round,
                    params: [inst.params[0], 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
            }

            OpCode::Onion => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Onion,
                    params: [inst.params[0], 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
            }

            OpCode::Elongate => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Elongate,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let ax = f32x8::splat(inst.params[0]);
                let ay = f32x8::splat(inst.params[1]);
                let az = f32x8::splat(inst.params[2]);

                // p - clamp(p, -amount, amount)
                p = Vec3x8 {
                    x: p.x - p.x.max(-ax).min(ax),
                    y: p.y - p.y.max(-ay).min(ay),
                    z: p.z - p.z.max(-az).min(az),
                };
            }

            OpCode::Noise => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Noise,
                    params: [inst.params[0], inst.params[1], inst.params[2], 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
                // Noise doesn't modify point, only post-processes distance
            }

            OpCode::Mirror => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Mirror,
                    params: [inst.params[0], inst.params[1], inst.params[2], 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                p = Vec3x8 {
                    x: if inst.params[0] != 0.0 {
                        p.x.abs()
                    } else {
                        p.x
                    },
                    y: if inst.params[1] != 0.0 {
                        p.y.abs()
                    } else {
                        p.y
                    },
                    z: if inst.params[2] != 0.0 {
                        p.z.abs()
                    } else {
                        p.z
                    },
                };
            }

            OpCode::OctantMirror => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::OctantMirror,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // abs
                let ax = p.x.abs();
                let ay = p.y.abs();
                let az = p.z.abs();

                // sort: x >= y >= z using min/max (branchless per lane)
                let max_xy = ax.max(ay);
                let min_xy = ax.min(ay);
                let x1 = max_xy.max(az);
                let min_xz = max_xy.min(az);
                let y1 = min_xy.max(min_xz);
                let z1 = min_xy.min(min_xz);

                p = Vec3x8 {
                    x: x1,
                    y: y1,
                    z: z1,
                };
            }

            OpCode::Revolution => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Revolution,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let offset = f32x8::splat(inst.params[0]);
                let q = (p.x * p.x + p.z * p.z).sqrt() - offset;
                p = Vec3x8 {
                    x: q,
                    y: p.y,
                    z: f32x8::ZERO,
                };
            }

            OpCode::Extrude => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Extrude,
                    params: [inst.params[0], 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                p = Vec3x8 {
                    x: p.x,
                    y: p.y,
                    z: f32x8::ZERO,
                };
            }

            OpCode::Taper => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Taper,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let factor = f32x8::splat(inst.params[0]);
                let s = f32x8::ONE / (f32x8::ONE - p.y * factor);
                p = Vec3x8 {
                    x: p.x * s,
                    y: p.y,
                    z: p.z * s,
                };
            }

            OpCode::Displacement => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Displacement,
                    // amplitude + per-axis frequency (legacy Displacement = 5,5,5)
                    params: [
                        inst.params[0],
                        inst.params[1],
                        inst.params[2],
                        inst.params[3],
                    ],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
                // Displacement doesn't modify point, only post-processes distance
            }

            OpCode::SweepBezier => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::SweepBezier,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let p0x = f32x8::splat(inst.params[0]);
                let p0z = f32x8::splat(inst.params[1]);
                let p1x = f32x8::splat(inst.params[2]);
                let p1z = f32x8::splat(inst.params[3]);
                let p2x = f32x8::splat(inst.params[4]);
                let p2z = f32x8::splat(inst.params[5]);
                let two = f32x8::splat(2.0);

                // Query point in XZ plane
                let qx = p.x;
                let qz = p.z;

                // Coarse search: 5 samples
                let mut best_t = f32x8::ZERO;
                let mut best_d2 = f32x8::splat(f32::MAX);
                for i in 0..5u32 {
                    let t = f32x8::splat(i as f32 * 0.25);
                    let omt = f32x8::ONE - t;
                    let omt2 = omt * omt;
                    let t2 = t * t;
                    let omt_t_2 = two * omt * t;
                    let bx = p0x * omt2 + p1x * omt_t_2 + p2x * t2;
                    let bz = p0z * omt2 + p1z * omt_t_2 + p2z * t2;
                    let dx = qx - bx;
                    let dz = qz - bz;
                    let d2 = dx * dx + dz * dz;
                    let mask = d2.cmp_lt(best_d2);
                    best_d2 = mask.blend(d2, best_d2);
                    best_t = mask.blend(t, best_t);
                }

                // B''(t) is constant for quadratic
                let bddx = two * (p0x - two * p1x + p2x);
                let bddz = two * (p0z - two * p1z + p2z);

                // Newton iterations
                let mut t = best_t;
                let eps = f32x8::splat(1e-10);
                for _ in 0..5 {
                    let omt = f32x8::ONE - t;
                    let omt2 = omt * omt;
                    let t2 = t * t;
                    let omt_t_2 = two * omt * t;
                    let bx = p0x * omt2 + p1x * omt_t_2 + p2x * t2;
                    let bz = p0z * omt2 + p1z * omt_t_2 + p2z * t2;
                    let tdx = (p1x - p0x) * (two * omt) + (p2x - p1x) * (two * t);
                    let tdz = (p1z - p0z) * (two * omt) + (p2z - p1z) * (two * t);
                    let diffx = bx - qx;
                    let diffz = bz - qz;
                    let num = diffx * tdx + diffz * tdz;
                    #[allow(clippy::suspicious_operation_groupings)]
                    let den = tdx * tdx + tdz * tdz + diffx * bddx + diffz * bddz;
                    let safe_den = den.abs().cmp_lt(eps).blend(f32x8::ONE, den);
                    t = (t - num / safe_den).max(f32x8::ZERO).min(f32x8::ONE);
                }

                // Distance to closest point on bezier
                let omt = f32x8::ONE - t;
                let cx = p0x * (omt * omt) + p1x * (two * omt * t) + p2x * (t * t);
                let cz = p0z * (omt * omt) + p1z * (two * omt * t) + p2z * (t * t);
                let dx = qx - cx;
                let dz = qz - cz;
                let d = (dx * dx + dz * dz).sqrt();

                p = Vec3x8 {
                    x: d,
                    y: p.y,
                    z: f32x8::ZERO,
                };
            }

            // ★ Deep Fried: SIMD PolarRepeat with atan2 approximation
            OpCode::PolarRepeat => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::PolarRepeat,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // Division Exorcism: params[1]=sector, params[2]=recip_sector
                let sector = f32x8::splat(inst.params[1]);
                let inv_sector = f32x8::splat(inst.params[2]);
                let angle = atan2_approx(p.z, p.x);
                let r = (p.x * p.x + p.z * p.z).sqrt();

                // remainder = angle - round(angle / sector) * sector (same as scalar law)
                let remainder = angle - (angle * inv_sector).round() * sector;

                p = Vec3x8 {
                    x: r * cos_approx(remainder),
                    y: p.y,
                    z: r * sin_approx(remainder),
                };
            }

            // === 2D Primitives (extruded) ===
            OpCode::Circle2D => {
                let r = f32x8::splat(inst.params[0]);
                let half_h = f32x8::splat(inst.params[1]);
                let d2d = (p.x * p.x + p.y * p.y).sqrt() - r;
                let dz = p.z.abs() - half_h;
                let wx = d2d.max(f32x8::ZERO);
                let wy = dz.max(f32x8::ZERO);
                let d = (wx * wx + wy * wy).sqrt() + d2d.max(dz).min(f32x8::ZERO);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Rect2D => {
                let hx = f32x8::splat(inst.params[0]);
                let hy = f32x8::splat(inst.params[1]);
                let half_h = f32x8::splat(inst.params[2]);
                let dx = p.x.abs() - hx;
                let dy = p.y.abs() - hy;
                let d2d = (dx.max(f32x8::ZERO) * dx.max(f32x8::ZERO)
                    + dy.max(f32x8::ZERO) * dy.max(f32x8::ZERO))
                .sqrt()
                    + dx.max(dy).min(f32x8::ZERO);
                let dz = p.z.abs() - half_h;
                let wx = d2d.max(f32x8::ZERO);
                let wy = dz.max(f32x8::ZERO);
                let d = (wx * wx + wy * wy).sqrt() + d2d.max(dz).min(f32x8::ZERO);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::RoundedRect2D => {
                let hx = f32x8::splat(inst.params[0]);
                let hy = f32x8::splat(inst.params[1]);
                let round_r = f32x8::splat(inst.params[2]);
                let half_h = f32x8::splat(inst.params[3]);
                let dx = p.x.abs() - hx + round_r;
                let dy = p.y.abs() - hy + round_r;
                let d2d = (dx.max(f32x8::ZERO) * dx.max(f32x8::ZERO)
                    + dy.max(f32x8::ZERO) * dy.max(f32x8::ZERO))
                .sqrt()
                    + dx.max(dy).min(f32x8::ZERO)
                    - round_r;
                let dz = p.z.abs() - half_h;
                let wx = d2d.max(f32x8::ZERO);
                let wy = dz.max(f32x8::ZERO);
                let d = (wx * wx + wy * wy).sqrt() + d2d.max(dz).min(f32x8::ZERO);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Annular2D => {
                let outer_r = f32x8::splat(inst.params[0]);
                let thickness = f32x8::splat(inst.params[1]);
                let half_h = f32x8::splat(inst.params[2]);
                let d2d = ((p.x * p.x + p.y * p.y).sqrt() - outer_r).abs() - thickness;
                let dz = p.z.abs() - half_h;
                let wx = d2d.max(f32x8::ZERO);
                let wy = dz.max(f32x8::ZERO);
                let d = (wx * wx + wy * wy).sqrt() + d2d.max(dz).min(f32x8::ZERO);
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Segment2D => {
                let a = glam::Vec2::new(inst.params[0], inst.params[1]);
                let b = glam::Vec2::new(inst.params[2], inst.params[3]);
                let (thickness, half_h) = (inst.params[4], inst.params[5]);
                let d = eval_per_lane(&p, |q| sdf_segment_2d(q, a, b, thickness, half_h));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            OpCode::Polygon2D => {
                // Vertices live in aux_data as flat [x0, y0, x1, y1, ...]
                let aux_off = inst.aux_offset as usize;
                let flat = &sdf.aux_data[aux_off..aux_off + inst.aux_len as usize];
                let half_h = inst.params[0];
                let d = eval_per_lane(&p, |q| sdf_polygon_2d_flat(q, flat, half_h));
                value_stack[vsp] = d * scale_correction;
                vsp += 1;
            }

            // === Exponential Smooth Operations ===
            // Per-lane scalar via the shared law in `crate::operations` (exact parity
            // with the tree / scalar evaluators). The previous Schraudolph exp + Padé ln
            // approximation was already per-lane scalar and its ln term was off by ~2x.
            OpCode::ExpSmoothUnion => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let k = inst.params[0];
                value_stack[vsp - 1] =
                    eval_per_lane_binary(a, b, |x, y| sdf_exp_smooth_union(x, y, k));
            }

            OpCode::ExpSmoothIntersection => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let k = inst.params[0];
                value_stack[vsp - 1] =
                    eval_per_lane_binary(a, b, |x, y| sdf_exp_smooth_intersection(x, y, k));
            }

            OpCode::ExpSmoothSubtraction => {
                vsp -= 1;
                let b = value_stack[vsp];
                let a = value_stack[vsp - 1];
                let k = inst.params[0];
                value_stack[vsp - 1] =
                    eval_per_lane_binary(a, b, |x, y| sdf_exp_smooth_subtraction(x, y, k));
            }

            // === New Modifiers ===
            OpCode::Shear => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Shear,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let sx = f32x8::splat(inst.params[0]);
                let sy = f32x8::splat(inst.params[1]);
                let sz = f32x8::splat(inst.params[2]);
                p = Vec3x8 {
                    x: p.x,
                    y: p.y - sx * p.x,
                    z: p.z - sy * p.x - sz * p.y,
                };
            }

            OpCode::Animated => {
                // Static SIMD evaluation: just pass through
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::Animated,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
            }

            // === New Transforms (3) ===
            OpCode::ProjectiveTransform => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::ProjectiveTransform,
                    params: [inst.params[0], 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // Per-lane projective transform using aux_data
                let aux_off = inst.aux_offset as usize;
                if inst.aux_len >= 16 {
                    let inv_matrix: &[f32] = &sdf.aux_data[aux_off..aux_off + 16];
                    let mut inv_m = [0.0f32; 16];
                    inv_m.copy_from_slice(inv_matrix);
                    let mut result = Vec3x8::zero();
                    for lane in 0..8 {
                        let pt = glam::Vec3::new(
                            p.x.as_array_ref()[lane],
                            p.y.as_array_ref()[lane],
                            p.z.as_array_ref()[lane],
                        );
                        let (q, _) =
                            crate::transforms::projective::projective_transform(pt, &inv_m);
                        result.x.as_array_mut()[lane] = q.x;
                        result.y.as_array_mut()[lane] = q.y;
                        result.z.as_array_mut()[lane] = q.z;
                    }
                    p = result;
                }
            }

            OpCode::LatticeDeform => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::LatticeDeform,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // Per-lane lattice deform using aux_data
                let aux_off = inst.aux_offset as usize;
                if inst.aux_len >= 9 {
                    let aux = &sdf.aux_data[aux_off..aux_off + inst.aux_len as usize];
                    let nx = aux[0] as u32;
                    let ny = aux[1] as u32;
                    let nz = aux[2] as u32;
                    let bbox_min = glam::Vec3::new(aux[3], aux[4], aux[5]);
                    let bbox_max = glam::Vec3::new(aux[6], aux[7], aux[8]);
                    let cp_data = &aux[9..];
                    let num_cps = cp_data.len() / 3;
                    let control_points: Vec<glam::Vec3> = (0..num_cps)
                        .map(|i| {
                            glam::Vec3::new(cp_data[i * 3], cp_data[i * 3 + 1], cp_data[i * 3 + 2])
                        })
                        .collect();
                    let mut result = Vec3x8::zero();
                    let mut corrections = [1.0f32; 8];
                    #[allow(clippy::needless_range_loop)]
                    for lane in 0..8 {
                        let pt = glam::Vec3::new(
                            p.x.as_array_ref()[lane],
                            p.y.as_array_ref()[lane],
                            p.z.as_array_ref()[lane],
                        );
                        let (q, correction) = crate::transforms::lattice::lattice_deform(
                            pt,
                            &control_points,
                            nx,
                            ny,
                            nz,
                            bbox_min,
                            bbox_max,
                        );
                        result.x.as_array_mut()[lane] = q.x;
                        result.y.as_array_mut()[lane] = q.y;
                        result.z.as_array_mut()[lane] = q.z;
                        corrections[lane] = correction;
                    }
                    p = result;
                    // Tree law: eval(child, q) / correction — applied at PopTransform
                    coord_stack[csp - 1].lane_correction = f32x8::new(corrections);
                }
            }

            OpCode::SdfSkinning => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::SdfSkinning,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                // Per-lane skinning using aux_data
                let aux_off = inst.aux_offset as usize;
                if inst.aux_len >= 1 {
                    let aux = &sdf.aux_data[aux_off..aux_off + inst.aux_len as usize];
                    let bone_count = aux[0] as usize;
                    let mut bones = Vec::with_capacity(bone_count);
                    let mut cursor = 1;
                    for _ in 0..bone_count {
                        let mut inv_bind = [0.0f32; 16];
                        let mut cur_pose = [0.0f32; 16];
                        inv_bind.copy_from_slice(&aux[cursor..cursor + 16]);
                        cursor += 16;
                        cur_pose.copy_from_slice(&aux[cursor..cursor + 16]);
                        cursor += 16;
                        let weight = aux[cursor];
                        cursor += 1;
                        bones.push(crate::transforms::skinning::BoneTransform {
                            inv_bind_pose: inv_bind,
                            current_pose: cur_pose,
                            weight,
                        });
                    }
                    let mut result = Vec3x8::zero();
                    for lane in 0..8 {
                        let pt = glam::Vec3::new(
                            p.x.as_array_ref()[lane],
                            p.y.as_array_ref()[lane],
                            p.z.as_array_ref()[lane],
                        );
                        let (q, _) = crate::transforms::skinning::sdf_skinning(pt, &bones);
                        result.x.as_array_mut()[lane] = q.x;
                        result.y.as_array_mut()[lane] = q.y;
                        result.z.as_array_mut()[lane] = q.z;
                    }
                    p = result;
                }
            }

            // === New Modifiers (4) ===
            OpCode::IcosahedralSymmetry => {
                // Per-lane scalar fallback: icosahedral fold has branches
                // that prevent vectorization, so dispatch per lane.
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::IcosahedralSymmetry,
                    params: [0.0; 4],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let mut result = Vec3x8::zero();
                for lane in 0..8 {
                    let pt = glam::Vec3::new(
                        p.x.as_array_ref()[lane],
                        p.y.as_array_ref()[lane],
                        p.z.as_array_ref()[lane],
                    );
                    let q = crate::modifiers::icosahedral_fold(pt);
                    result.x.as_array_mut()[lane] = q.x;
                    result.y.as_array_mut()[lane] = q.y;
                    result.z.as_array_mut()[lane] = q.z;
                }
                p = result;
            }

            OpCode::IFS => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::IFS,
                    params: [inst.params[0], 0.0, 0.0, 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;

                let iterations = inst.params[0] as u32;
                let aux_off = inst.aux_offset as usize;
                if inst.aux_len >= 1 {
                    // Deserialize transforms from aux_data
                    let aux = &sdf.aux_data[aux_off..aux_off + inst.aux_len as usize];
                    let transform_count = aux[0] as usize;
                    let mut transforms = Vec::with_capacity(transform_count);
                    for i in 0..transform_count {
                        let base = 1 + i * 16;
                        let mut mat = [0.0f32; 16];
                        mat.copy_from_slice(&aux[base..base + 16]);
                        transforms.push(mat);
                    }
                    let mut result = Vec3x8::zero();
                    for lane in 0..8 {
                        let pt = glam::Vec3::new(
                            p.x.as_array_ref()[lane],
                            p.y.as_array_ref()[lane],
                            p.z.as_array_ref()[lane],
                        );
                        let (q, _) =
                            crate::modifiers::ifs_fold_with_scale(pt, &transforms, iterations);
                        result.x.as_array_mut()[lane] = q.x;
                        result.y.as_array_mut()[lane] = q.y;
                        result.z.as_array_mut()[lane] = q.z;
                    }
                    p = result;
                } else {
                    // Fallback: simple abs-fold contraction
                    let mut result = Vec3x8::zero();
                    for lane in 0..8 {
                        let mut px = p.x.as_array_ref()[lane];
                        let mut py = p.y.as_array_ref()[lane];
                        let mut pz = p.z.as_array_ref()[lane];
                        for _ in 0..iterations {
                            px = px.abs() * 0.5;
                            py = py.abs() * 0.5;
                            pz = pz.abs() * 0.5;
                        }
                        result.x.as_array_mut()[lane] = px;
                        result.y.as_array_mut()[lane] = py;
                        result.z.as_array_mut()[lane] = pz;
                    }
                    p = result;
                }
            }

            OpCode::HeightmapDisplacement => {
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::HeightmapDisplacement,
                    params: [inst.params[0], inst.params[1], 0.0, 0.0],
                    aux_offset: inst.aux_offset,
                    aux_len: inst.aux_len,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
            }

            OpCode::SurfaceRoughness => {
                // Per-lane scalar fallback for surface roughness
                // Note: This is a post-process modifier using FBM noise
                coord_stack[csp] = CoordFrameSimd {
                    point: p,
                    scale_correction,
                    opcode: OpCode::SurfaceRoughness,
                    params: [inst.params[0], inst.params[1], inst.params[2], 0.0],
                    aux_offset: 0,
                    aux_len: 0,
                    lane_correction: f32x8::ONE,
                };
                csp += 1;
                // No point transformation - roughness is applied after distance eval
            }

            // === Control ===
            OpCode::PopTransform => {
                csp -= 1;
                let frame = coord_stack[csp];

                match frame.opcode {
                    OpCode::Round => {
                        let radius = f32x8::splat(frame.params[0]);
                        value_stack[vsp - 1] -= radius;
                    }
                    OpCode::Onion => {
                        let thickness = f32x8::splat(frame.params[0]);
                        value_stack[vsp - 1] = value_stack[vsp - 1].abs() - thickness;
                    }
                    OpCode::Noise => {
                        // Process noise for each lane
                        let amplitude = frame.params[0];
                        let frequency = frame.params[1];
                        let seed = frame.params[2] as u32;
                        let mut noise_vals = [0.0f32; 8];
                        #[allow(clippy::needless_range_loop)]
                        for i in 0..8 {
                            let pt = glam::Vec3::new(
                                frame.point.x.as_array_ref()[i] * frequency,
                                frame.point.y.as_array_ref()[i] * frequency,
                                frame.point.z.as_array_ref()[i] * frequency,
                            );
                            noise_vals[i] = perlin_noise_3d(pt.x, pt.y, pt.z, seed) * amplitude;
                        }
                        value_stack[vsp - 1] += f32x8::new(noise_vals);
                    }
                    OpCode::Extrude => {
                        let half_height = f32x8::splat(frame.params[0]);
                        let child_dist = value_stack[vsp - 1];
                        let w_y = frame.point.z.abs() - half_height;
                        let w_x_pos = child_dist.max(f32x8::ZERO);
                        let w_y_pos = w_y.max(f32x8::ZERO);
                        let outside = (w_x_pos * w_x_pos + w_y_pos * w_y_pos).sqrt();
                        let inside = child_dist.max(w_y).min(f32x8::ZERO);
                        value_stack[vsp - 1] = outside + inside;
                    }
                    OpCode::Displacement => {
                        let amplitude = frame.params[0];
                        let frequency =
                            glam::Vec3::new(frame.params[1], frame.params[2], frame.params[3]);
                        let d = value_stack[vsp - 1];
                        let da = d.as_array_ref();
                        let mut out = [0.0f32; 8];
                        #[allow(clippy::needless_range_loop)]
                        for i in 0..8 {
                            let pt = glam::Vec3::new(
                                frame.point.x.as_array_ref()[i],
                                frame.point.y.as_array_ref()[i],
                                frame.point.z.as_array_ref()[i],
                            );
                            out[i] = crate::modifiers::modifier_sine_displacement(
                                da[i], pt, amplitude, frequency,
                            );
                        }
                        value_stack[vsp - 1] = f32x8::new(out);
                    }
                    OpCode::HeightmapDisplacement => {
                        let amplitude = frame.params[0];
                        let scale = frame.params[1];
                        let aux_off = frame.aux_offset as usize;
                        // Same law as the scalar evaluator: d -= heightmap_displacement(...)
                        if frame.aux_len >= 2 {
                            let aux = &sdf.aux_data[aux_off..aux_off + frame.aux_len as usize];
                            let w = aux[0] as u32;
                            let h = aux[1] as u32;
                            let hmap = &aux[2..];
                            let mut displacement_vals = [0.0f32; 8];
                            #[allow(clippy::needless_range_loop)]
                            for i in 0..8 {
                                let pt = glam::Vec3::new(
                                    frame.point.x.as_array_ref()[i],
                                    frame.point.y.as_array_ref()[i],
                                    frame.point.z.as_array_ref()[i],
                                );
                                displacement_vals[i] = crate::modifiers::heightmap_displacement(
                                    pt, hmap, w, h, amplitude, scale,
                                );
                            }
                            value_stack[vsp - 1] -= f32x8::new(displacement_vals);
                        }
                    }
                    OpCode::SurfaceRoughness => {
                        // FBM noise for surface roughness
                        let frequency = frame.params[0];
                        let amplitude = frame.params[1];
                        let octaves = frame.params[2] as u32;
                        // Same law as the scalar evaluator (crate::modifiers::surface_roughness)
                        let d = value_stack[vsp - 1];
                        let da = d.as_array_ref();
                        let mut out = [0.0f32; 8];
                        #[allow(clippy::needless_range_loop)]
                        for i in 0..8 {
                            let pt = glam::Vec3::new(
                                frame.point.x.as_array_ref()[i],
                                frame.point.y.as_array_ref()[i],
                                frame.point.z.as_array_ref()[i],
                            );
                            out[i] = crate::modifiers::surface_roughness(
                                pt, da[i], frequency, amplitude, octaves,
                            );
                        }
                        value_stack[vsp - 1] = f32x8::new(out);
                    }
                    OpCode::ProjectiveTransform => {
                        // Lipschitz correction for projective transform
                        let lipschitz_bound = f32x8::splat(frame.params[0]);
                        value_stack[vsp - 1] *= lipschitz_bound;
                    }
                    OpCode::LatticeDeform => {
                        value_stack[vsp - 1] /= frame.lane_correction;
                    }
                    _ => {}
                }

                p = frame.point;
                scale_correction = frame.scale_correction;
            }

            OpCode::End => {
                break;
            }
        }
    }

    if vsp > 0 {
        value_stack[0]
    } else {
        f32x8::splat(f32::MAX)
    }
}

/// Smooth min — Division Exorcism edition.
/// Takes precomputed `rk = 1/k` to eliminate SIMD division.
#[inline(always)]
fn smooth_min_simd_rk(a: f32x8, b: f32x8, k: f32x8, rk: f32x8) -> f32x8 {
    let h = (f32x8::ONE - (a - b).abs() * rk).max(f32x8::ZERO);
    a.min(b) - h * h * k * f32x8::splat(0.25)
}

/// Chamfer minimum for SIMD: 45-degree beveled blend
/// `min(min(a,b), (a + b) * FRAC_1_SQRT_2 - r)`
#[inline(always)]
fn chamfer_min_simd(a: f32x8, b: f32x8, r: f32x8) -> f32x8 {
    let s = f32x8::splat(std::f32::consts::FRAC_1_SQRT_2);
    a.min(b).min((a + b) * s - r)
}

/// Stairs minimum for SIMD: stepped/terraced blend (Mercury hg_sdf)
/// Processes 8 lanes with scalar params r and n
#[inline(always)]
fn stairs_min_simd(a: f32x8, b: f32x8, r: f32, n: f32) -> f32x8 {
    let s = f32x8::splat(std::f32::consts::FRAC_1_SQRT_2);
    let s2 = f32x8::splat(std::f32::consts::SQRT_2);
    let half = f32x8::splat(0.5);
    // Division Exorcism: precompute scalar reciprocals, then splat
    let rn = r / n; // scalar division (1 cycle amortized over 8 lanes)
    let rn_v = f32x8::splat(rn);
    let step = r * std::f32::consts::SQRT_2 / n;
    let inv_step = 1.0 / step; // scalar reciprocal
    let step_v = f32x8::splat(step);
    let inv_step_v = f32x8::splat(inv_step);
    let hs_v = step_v * half;
    let off_v = f32x8::splat((r - rn) * 0.5 * std::f32::consts::SQRT_2);
    let edge_v = f32x8::splat(0.5 * rn);

    let d = a.min(b);

    // pR45
    let mut px = (a + b) * s;
    let mut py = (b - a) * s;
    // swap
    std::mem::swap(&mut px, &mut py);

    px -= off_v;
    py -= off_v;
    px += half * s2 * rn_v;

    // pMod1: px = glsl_mod(px + hs, step) - hs — Division Exorcism: multiply by inv_step
    let t = px + hs_v;
    px = t - step_v * (t * inv_step_v).floor() - hs_v;

    let d = d.min(py);

    // Second pR45
    let npx = (px + py) * s;
    let npy = (py - px) * s;

    d.min((npx - edge_v).max(npy - edge_v))
}

/// SIMD cosine (`wide` polynomial with range reduction, ~1e-6 abs error).
///
/// Replaced the Bhaskara I approximation (1.6e-3 abs error) in 1.9.1: TPMS
/// primitives and twist/bend amplified that error to >10% vs the tree law.
#[inline(always)]
fn cos_approx(x: f32x8) -> f32x8 {
    x.cos()
}

/// SIMD sine (`wide` polynomial with range reduction, ~1e-6 abs error).
#[inline(always)]
fn sin_approx(x: f32x8) -> f32x8 {
    x.sin()
}

/// SIMD atan2 (`wide` implementation, ~1e-6 rad).
///
/// Replaced the minimax polynomial (3.8e-3 rad) in 1.9.1 for parity with the
/// scalar `modifier_polar_repeat_rk` law.
#[inline(always)]
fn atan2_approx(y: f32x8, x: f32x8) -> f32x8 {
    y.atan2(x)
}

/// Per-lane scalar evaluation helper for binary operations
///
/// Evaluates a scalar binary operation for each of the 8 SIMD lanes independently.
#[inline(always)]
fn eval_per_lane_binary(a: f32x8, b: f32x8, f: impl Fn(f32, f32) -> f32) -> f32x8 {
    let aa = a.as_array_ref();
    let ba = b.as_array_ref();
    let mut results = [0.0f32; 8];
    for i in 0..8 {
        results[i] = f(aa[i], ba[i]);
    }
    f32x8::new(results)
}

/// Per-lane scalar evaluation helper for extended primitives
///
/// Evaluates a scalar SDF function for each of the 8 SIMD lanes independently.
/// This is used for complex primitives where full SIMD implementation would be
/// error-prone. The SIMD benefit still comes from parallelizing tree traversal
/// (transforms, operations) across 8 points.
#[inline(always)]
fn eval_per_lane(p: &Vec3x8, f: impl Fn(Vec3) -> f32) -> f32x8 {
    let px = p.x.as_array_ref();
    let py = p.y.as_array_ref();
    let pz = p.z.as_array_ref();
    let mut results = [0.0f32; 8];
    for i in 0..8 {
        results[i] = f(Vec3::new(px[i], py[i], pz[i]));
    }
    f32x8::new(results)
}

/// Batch evaluate compiled SDF using SIMD (8 points at a time)
///
/// This is the main entry point for SIMD-accelerated batch evaluation.
/// Points are processed 8 at a time, with any remainder handled by
/// scalar evaluation.
pub fn eval_compiled_batch_simd(sdf: &CompiledSdf, points: &[Vec3]) -> Vec<f32> {
    let n = points.len();
    let mut results = vec![0.0f32; n];

    // Process 8 points at a time
    let chunks = n / 8;
    for i in 0..chunks {
        let base = i * 8;
        let p = Vec3x8::from_vecs([
            points[base],
            points[base + 1],
            points[base + 2],
            points[base + 3],
            points[base + 4],
            points[base + 5],
            points[base + 6],
            points[base + 7],
        ]);

        let d = eval_compiled_simd(sdf, p);
        let arr = d.to_array();
        results[base..base + 8].copy_from_slice(&arr);
    }

    // Handle remainder with scalar evaluation
    let remainder = n % 8;
    if remainder > 0 {
        let base = chunks * 8;
        for i in 0..remainder {
            results[base + i] = super::eval::eval_compiled(sdf, points[base + i]);
        }
    }

    results
}

/// Parallel batch evaluate compiled SDF using SIMD
///
/// Combines SIMD (8-wide) with multi-threading for maximum performance.
pub fn eval_compiled_batch_simd_parallel(sdf: &CompiledSdf, points: &[Vec3]) -> Vec<f32> {
    use rayon::prelude::*;

    let n = points.len();
    if n < 64 {
        // Not worth parallelizing for small inputs
        return eval_compiled_batch_simd(sdf, points);
    }

    // Process in chunks of 64 (8 SIMD lanes * 8 iterations per thread)
    let chunk_size = 64;
    let mut results = vec![0.0f32; n];

    results
        .par_chunks_mut(chunk_size)
        .enumerate()
        .for_each(|(chunk_idx, chunk)| {
            let base = chunk_idx * chunk_size;
            let chunk_len = chunk.len();

            // Process 8 points at a time within this chunk
            let simd_iters = chunk_len / 8;
            for i in 0..simd_iters {
                let p_base = base + i * 8;
                let p = Vec3x8::from_vecs([
                    points[p_base],
                    points[p_base + 1],
                    points[p_base + 2],
                    points[p_base + 3],
                    points[p_base + 4],
                    points[p_base + 5],
                    points[p_base + 6],
                    points[p_base + 7],
                ]);

                let d = eval_compiled_simd(sdf, p);
                let arr = d.to_array();
                chunk[i * 8..(i + 1) * 8].copy_from_slice(&arr);
            }

            // Handle remainder within chunk
            let remainder = chunk_len % 8;
            if remainder > 0 {
                let r_base = simd_iters * 8;
                for i in 0..remainder {
                    chunk[r_base + i] = super::eval::eval_compiled(sdf, points[base + r_base + i]);
                }
            }
        });

    results
}

// ============================================================================
// Gradient (Normal) Computation - SIMD
// ============================================================================

/// Compute gradient (normal direction) using finite differences - SIMD 8-wide
///
/// Returns unnormalized gradient vectors (gx, gy, gz) for 8 points simultaneously.
/// The caller can normalize if needed.
///
/// # Performance
/// - 6 SDF evaluations per call (central differences)
/// - ~6x slower than distance-only, but still SIMD-accelerated
/// - For 100K points: ~0.6ms on modern CPU
#[inline]
pub fn eval_gradient_simd(sdf: &CompiledSdf, p: Vec3x8, epsilon: f32) -> (f32x8, f32x8, f32x8) {
    let e = f32x8::splat(epsilon);
    let ne = f32x8::splat(-epsilon);

    // Tetrahedral method: 4 evaluations instead of 6
    let v0 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + e,
            y: p.y + ne,
            z: p.z + ne,
        },
    ); // (+,-,-)
    let v1 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + ne,
            y: p.y + ne,
            z: p.z + e,
        },
    ); // (-,-,+)
    let v2 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + ne,
            y: p.y + e,
            z: p.z + ne,
        },
    ); // (-,+,-)
    let v3 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + e,
            y: p.y + e,
            z: p.z + e,
        },
    ); // (+,+,+)

    let gx = v0 - v1 - v2 + v3;
    let gy = -v0 - v1 + v2 + v3;
    let gz = -v0 + v1 - v2 + v3;

    (gx, gy, gz)
}

/// Compute both distance and gradient in one call - SIMD 8-wide
///
/// More efficient when you need both values, as it shares the center evaluation.
#[inline]
pub fn eval_distance_and_gradient_simd(
    sdf: &CompiledSdf,
    p: Vec3x8,
    epsilon: f32,
) -> (f32x8, f32x8, f32x8, f32x8) {
    // Center distance
    let d_center = eval_compiled_simd(sdf, p);

    // Tetrahedral method: 4 evaluations for gradient
    let e = f32x8::splat(epsilon);
    let ne = f32x8::splat(-epsilon);

    let v0 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + e,
            y: p.y + ne,
            z: p.z + ne,
        },
    ); // (+,-,-)
    let v1 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + ne,
            y: p.y + ne,
            z: p.z + e,
        },
    ); // (-,-,+)
    let v2 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + ne,
            y: p.y + e,
            z: p.z + ne,
        },
    ); // (-,+,-)
    let v3 = eval_compiled_simd(
        sdf,
        Vec3x8 {
            x: p.x + e,
            y: p.y + e,
            z: p.z + e,
        },
    ); // (+,+,+)

    let gx = v0 - v1 - v2 + v3;
    let gy = -v0 - v1 + v2 + v3;
    let gz = -v0 + v1 - v2 + v3;

    (d_center, gx, gy, gz)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::compiled::CompiledSdf;
    use crate::eval::eval;
    use crate::types::SdfNode;

    #[test]
    fn test_simd_sphere() {
        let node = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&node);

        let points = Vec3x8::from_vecs([
            Vec3::ZERO,
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(2.0, 0.0, 0.0),
            Vec3::new(0.0, 1.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            Vec3::new(0.5, 0.5, 0.5),
            Vec3::new(-1.0, 0.0, 0.0),
            Vec3::new(0.0, -1.0, 0.0),
        ]);

        let d = eval_compiled_simd(&compiled, points);
        let arr = d.to_array();

        // Check against scalar evaluation
        assert!((arr[0] - (-1.0)).abs() < 0.001); // origin
        assert!(arr[1].abs() < 0.001); // on surface
        assert!((arr[2] - 1.0).abs() < 0.001); // outside
    }

    #[test]
    fn test_simd_box() {
        let node = SdfNode::box3d(1.0, 1.0, 1.0);
        let compiled = CompiledSdf::compile(&node);

        let points = Vec3x8::from_vecs([
            Vec3::ZERO,
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(2.0, 0.0, 0.0),
            Vec3::new(0.5, 0.5, 0.5),
            Vec3::new(1.5, 1.5, 0.0),
            Vec3::new(-0.5, -0.5, -0.5),
            Vec3::new(0.0, 0.0, 1.0),
            Vec3::new(0.0, 0.0, 2.0),
        ]);

        let d = eval_compiled_simd(&compiled, points);
        let arr = d.to_array();

        // Compare with scalar
        let test_vecs = [
            Vec3::ZERO,
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(2.0, 0.0, 0.0),
        ];
        for (i, v) in test_vecs.iter().enumerate() {
            let scalar = eval(&node, *v);
            assert!(
                (arr[i] - scalar).abs() < 0.001,
                "Mismatch at {}: simd={}, scalar={}",
                i,
                arr[i],
                scalar
            );
        }
    }

    #[test]
    fn test_simd_union() {
        let node = SdfNode::sphere(1.0).union(SdfNode::sphere(1.0).translate(2.0, 0.0, 0.0));
        let compiled = CompiledSdf::compile(&node);

        let points = Vec3x8::from_vecs([
            Vec3::ZERO,
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(2.0, 0.0, 0.0),
            Vec3::new(3.0, 0.0, 0.0),
            Vec3::new(1.0, 1.0, 0.0),
            Vec3::new(-1.0, 0.0, 0.0),
            Vec3::new(0.0, 2.0, 0.0),
            Vec3::new(2.0, 2.0, 0.0),
        ]);

        let d = eval_compiled_simd(&compiled, points);
        let arr = d.to_array();

        // Verify against scalar
        let test_vecs = [
            Vec3::ZERO,
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(2.0, 0.0, 0.0),
        ];
        for (i, v) in test_vecs.iter().enumerate() {
            let scalar = eval(&node, *v);
            assert!(
                (arr[i] - scalar).abs() < 0.01,
                "Union mismatch at {}: simd={}, scalar={}",
                i,
                arr[i],
                scalar
            );
        }
    }

    #[test]
    fn test_simd_batch() {
        let node = SdfNode::sphere(1.0).smooth_union(SdfNode::box3d(0.8, 0.8, 0.8), 0.1);
        let compiled = CompiledSdf::compile(&node);

        let points: Vec<Vec3> = (0..100)
            .map(|i| {
                let t = i as f32 / 100.0;
                Vec3::new(
                    (t * 12.34).sin() * 2.0,
                    (t * 23.45).sin() * 2.0,
                    (t * 34.56).sin() * 2.0,
                )
            })
            .collect();

        let simd_results = eval_compiled_batch_simd(&compiled, &points);
        let scalar_results: Vec<f32> = points.iter().map(|p| eval(&node, *p)).collect();

        for (i, (simd, scalar)) in simd_results.iter().zip(scalar_results.iter()).enumerate() {
            assert!(
                (simd - scalar).abs() < 0.01,
                "Batch mismatch at {}: simd={}, scalar={}",
                i,
                simd,
                scalar
            );
        }
    }

    #[test]
    fn test_simd_complex() {
        let node = SdfNode::sphere(1.0)
            .smooth_union(SdfNode::box3d(0.8, 0.8, 0.8), 0.1)
            .translate(0.5, 0.0, 0.0)
            .scale(1.5);
        let compiled = CompiledSdf::compile(&node);

        let test_vecs = [
            Vec3::ZERO,
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(0.5, 0.5, 0.5),
            Vec3::new(-1.0, -1.0, -1.0),
        ];

        let points = Vec3x8::from_vecs([
            test_vecs[0],
            test_vecs[1],
            test_vecs[2],
            test_vecs[3],
            Vec3::ZERO,
            Vec3::ZERO,
            Vec3::ZERO,
            Vec3::ZERO,
        ]);

        let d = eval_compiled_simd(&compiled, points);
        let arr = d.to_array();

        for (i, v) in test_vecs.iter().enumerate() {
            let scalar = eval(&node, *v);
            assert!(
                (arr[i] - scalar).abs() < 0.05,
                "Complex mismatch at {:?}: simd={}, scalar={}",
                v,
                arr[i],
                scalar
            );
        }
    }
}