alice-sdf 2.1.0

//! Sphere tracing (raymarching) implementation (Deep Fried v2)
//!
//! # Deep Fried v2 Optimizations
//!
//! - **Integer Division Exorcism**: Row-based iteration eliminates `i % w` / `i / w`.
//! - **Incremental Ray Generation**: Pre-calculate row vectors per pixel.
//! - **Forced Inlining**: Hot loops are forced inline.
//! - **Over-Relaxation**: Enhanced sphere tracing with safety checks.
//! - **Compiled Backend**: `CompiledSdf` VM eliminates tree traversal overhead.
//! - **JIT Backend**: `JitCompiledSdf` native code, zero interpreter overhead.
//! - **SIMD Packet Tracing**: 8 rays marched simultaneously via `Vec3x8`/`JitSimdSdf`.
//! - **Tile-Based Rendering**: Cache-friendly 8-pixel horizontal tiles.
//!
//! Author: Moroya Sakamoto

use crate::eval::{eval, normal};
use crate::interval::eval_lipschitz;
use crate::types::{Hit, Ray};
use crate::SdfNode;
use glam::Vec3;
use rayon::prelude::*;

use crate::compiled::{
    eval_compiled, eval_compiled_normal, eval_compiled_simd, CompiledSdf, Vec3x8,
};
use wide::{f32x8, CmpGe, CmpGt, CmpLt};

/// Raymarch configuration
#[derive(Debug, Clone, Copy)]
pub struct RaymarchConfig {
    /// Maximum number of marching steps. A ray that exhausts the budget is
    /// reported as a miss, so the budget is a distance budget: with steps of
    /// `d / lipschitz`, [`with_bound`](Self::with_bound) scales it by the
    /// bound (a TPMS at L = √3 gets √3 × the steps of a unit-Lipschitz
    /// field). The default (256) leaves ≲ 0.1 % of random rays through a
    /// thin TPMS shell unresolved — rays that graze the shell and creep by
    /// `≈ ε` per step; raise it or use `high_quality` for those.
    pub max_steps: u32,
    /// Distance threshold for surface hit
    pub epsilon: f32,
    /// Minimum step in field units (divided by `lipschitz` like `d`);
    /// prevents infinite loops. Keep `min_step ≤ epsilon` so the floor
    /// cannot carry a sample across the `|d| < epsilon` hit band.
    pub min_step: f32,
    /// Step relaxation factor (omega for over-relaxation, ω ∈ [1, 2))
    pub omega: f32,
    /// Lipschitz constant of the SDF (≥ 1.0). Step = d * ω / L.
    pub lipschitz: f32,
}

impl Default for RaymarchConfig {
    fn default() -> Self {
        Self {
            max_steps: 256,
            epsilon: 0.0001,
            min_step: 0.0001,
            omega: 1.0,
            lipschitz: 1.0,
        }
    }
}

impl RaymarchConfig {
    /// High quality configuration
    pub const fn high_quality() -> Self {
        Self {
            max_steps: 512,
            epsilon: 0.00001,
            min_step: 0.00001,
            omega: 1.0,
            lipschitz: 1.0,
        }
    }

    /// Fast configuration for preview (Deep Fried)
    ///
    /// Uses moderate over-relaxation for faster convergence while maintaining stability.
    pub const fn fast() -> Self {
        Self {
            max_steps: 64,
            epsilon: 0.0005,
            min_step: 0.0001,
            omega: 1.2,
            lipschitz: 1.0,
        }
    }

    /// This configuration with the Lipschitz bound raised to `bound` when
    /// `bound` is finite and larger than the current one.
    ///
    /// Every marcher applies the bound of the field it traces — the tree's
    /// `eval_lipschitz` or the value a `CompiledSdf` / `JitCompiledSdf`
    /// recorded at compile time — so a TPMS surface (L up to 7) traces
    /// correctly with `RaymarchConfig::default()`. An infinite bound (no
    /// provable step size) leaves the configuration unchanged.
    ///
    /// The step budget grows with the bound (`max_steps · bound / lipschitz`,
    /// rounded up): steps are `d / L`, so the same number of steps covers
    /// `1 / L` of the distance. Until 2.1.0 the budget stayed at 128 and a
    /// gyroid with the default configuration lost ~0.5 % of random rays to
    /// budget exhaustion.
    #[inline]
    #[must_use]
    pub fn with_bound(mut self, bound: f32) -> Self {
        if bound.is_finite() && bound > self.lipschitz {
            let scale = bound / self.lipschitz;
            self.max_steps = ((self.max_steps as f32) * scale).ceil() as u32;
            self.lipschitz = bound;
        }
        self
    }

    /// Relaxed sphere tracing configuration.
    ///
    /// Auto-computes the Lipschitz bound from the SDF tree
    /// ([`eval_lipschitz`]) and uses ω = 1.6 over-relaxation with safe step
    /// `d · ω / L`; fewer steps than plain tracing at grazing incidence.
    ///
    /// When the tree has no finite Lipschitz bound (`eval_lipschitz` returns
    /// `INFINITY`: domain repetition, taper, columns, …) no step size is
    /// provably safe, so this falls back to plain sphere tracing (ω = 1,
    /// L = 1) — the same best-effort behaviour as [`raymarch`].
    pub fn relaxed(node: &SdfNode) -> Self {
        let lip = eval_lipschitz(node);
        let (omega, lipschitz) = if lip.is_finite() {
            (1.6, lip.max(1.0))
        } else {
            (1.0, 1.0)
        };
        Self {
            // same distance budget as `default().with_bound(lipschitz)`
            max_steps: (256.0 * lipschitz).ceil() as u32,
            epsilon: 0.0001,
            min_step: 0.0001,
            omega,
            lipschitz,
        }
    }
}

/// Over-relaxed sphere tracing step law (Keinert et al. 2014, "Enhanced
/// Sphere Tracing", §3.1) — the single implementation shared by the tree,
/// compiled and JIT marchers so the three paths cannot drift.
///
/// Plain sphere tracing steps by `d / L`. Over-relaxation steps by `ω·d / L`
/// (ω ∈ [1, 2)) and checks at the next sample whether the two unbounding
/// spheres still overlap: if `|d| + |d_prev| < step_prev` (scaled by `1/L`)
/// they do not, so the relaxed step may have skipped the surface. The only
/// sound response is to **go back** to the last provably safe position
/// `t_prev + |d_prev| / L` (the boundary of the previous unbounding sphere)
/// and continue with ω = 1 for the rest of the ray. Advancing from the
/// overshot position instead (what 1.10.x did) leaves the ray behind the
/// surface with `d < 0`, where every step collapses to `min_step` and the
/// march exhausts `max_steps` — a unit sphere lost 76 % of its rays at ω = 1.6.
struct RelaxedStepper {
    /// Current relaxation factor; drops to 1 after the first detected overshoot.
    omega: f32,
    inv_lip: f32,
    min_step: f32,
    prev_t: f32,
    /// `|d_prev| / L`: radius of the previous unbounding sphere.
    prev_radius: f32,
    /// Length of the previous (possibly relaxed) step; 0 before the first step.
    prev_step: f32,
}

impl RelaxedStepper {
    #[inline(always)]
    fn new(config: &RaymarchConfig) -> Self {
        let inv_lip = 1.0 / config.lipschitz.max(1.0);
        Self {
            omega: config.omega.max(1.0),
            inv_lip,
            // The floor lives in field units like `d`: a floor of `min_step`
            // in ray units would move the sample by up to `L·min_step` in
            // field value and, for `L·min_step > ε`, jump over the `|d| < ε`
            // hit band into the interior (Neovius, L = 7, lost 16 % of its
            // rays that way). With `min_step ≤ ε` the next sample satisfies
            // `d' ≥ d − min_step > 0`, so plain tracing never overshoots.
            min_step: config.min_step * inv_lip,
            prev_t: 0.0,
            prev_radius: 0.0,
            prev_step: 0.0,
        }
    }

    /// Overshoot check for the signed distance `d` sampled after the previous
    /// step. Returns the ray parameter to retreat to when the relaxed step was
    /// not covered — the unbounding spheres of the two samples do not overlap,
    /// or the sign of `d` flipped (the ray is now behind the surface; with an
    /// exact SDF the spheres then touch exactly at the surface, which `<`
    /// alone misses) — and switches the rest of the ray to plain tracing.
    ///
    /// Must run **before** the `|d| < ε` hit test: a relaxed step can land
    /// just past a thin feature with `|d| < ε` on the far side, which is not
    /// a hit (Keinert et al. 2014: `if (!sorFail && radius < ε) break;`).
    #[inline(always)]
    fn overshoot(&mut self, d: f32) -> Option<f32> {
        let radius = d.abs() * self.inv_lip;
        if self.omega > 1.0
            && self.prev_step > 0.0
            && (radius + self.prev_radius < self.prev_step || d < 0.0)
        {
            self.omega = 1.0;
            self.prev_step = 0.0;
            return Some(self.prev_t + self.prev_radius);
        }
        None
    }

    /// Next ray parameter given the signed distance `d` sampled at `t`
    /// (call after [`Self::overshoot`] returned `None`).
    #[inline(always)]
    fn advance(&mut self, t: f32, d: f32) -> f32 {
        let step = (d * self.inv_lip * self.omega).max(self.min_step);
        self.prev_t = t;
        self.prev_radius = d.abs() * self.inv_lip;
        self.prev_step = step;
        t + step
    }
}

/// Raymarch result with detailed information
#[derive(Debug, Clone, Copy)]
pub struct RaymarchResult {
    /// Whether a hit was found
    pub hit: bool,
    /// Distance traveled
    pub distance: f32,
    /// Hit point (if hit)
    pub point: Vec3,
    /// Surface normal (if hit)
    pub normal: Vec3,
    /// Number of steps taken
    pub steps: u32,
}

// ============================================================================
// Interpreter Backend (SdfNode tree traversal)
// ============================================================================

/// Perform sphere tracing along a ray (Interpreter path)
///
/// Uses `eval()` interpreter to traverse the SDF tree. This is the simplest
/// path but slowest. Prefer `raymarch_compiled()` or `raymarch_jit()` for
/// production use.
///
/// Steps by `d / L` with `L = eval_lipschitz(node)` computed per call (a
/// tree walk; a twist / bend child costs an AABB pass) — for many rays use
/// [`raymarch_batch`] / [`render_depth`], which compute it once, or the
/// compiled marchers, which record it at compile time.
#[inline(always)]
pub fn raymarch(node: &SdfNode, origin: Vec3, direction: Vec3, max_distance: f32) -> Option<Hit> {
    let config = RaymarchConfig::default().with_bound(eval_lipschitz(node));
    raymarch_with_config(node, origin, direction, max_distance, &config)
}

/// Perform sphere tracing with custom configuration
///
/// Supports relaxed sphere tracing via `omega` and `lipschitz` fields:
/// step = `ω·d / L`; when the unbounding spheres of two consecutive samples
/// stop overlapping the ray retreats to the last safe point and continues
/// unrelaxed (Keinert et al. 2014).
#[inline(always)]
pub fn raymarch_with_config(
    node: &SdfNode,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
    config: &RaymarchConfig,
) -> Option<Hit> {
    let dir = direction.normalize();
    let mut t = 0.0;
    let mut steps = 0;
    let mut stepper = RelaxedStepper::new(config);

    while t < max_distance && steps < config.max_steps {
        let point = origin + dir * t;
        let d = eval(node, point);

        if let Some(back) = stepper.overshoot(d) {
            t = back;
            steps += 1;
            continue;
        }

        if d.abs() < config.epsilon {
            return Some(Hit {
                distance: t,
                point,
                normal: normal(node, point, config.epsilon),
                steps,
            });
        }

        t = stepper.advance(t, d);
        steps += 1;
    }

    None
}

/// Relaxed sphere tracing with auto-computed Lipschitz constant.
///
/// Uses ω=1.6 over-relaxation with Lipschitz-safe step sizes.
/// Typically converges in ~40% fewer steps than standard tracing for exact SDFs.
#[inline(always)]
pub fn raymarch_relaxed(
    node: &SdfNode,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
) -> Option<Hit> {
    let config = RaymarchConfig::relaxed(node);
    raymarch_with_config(node, origin, direction, max_distance, &config)
}

/// Raymarch returning detailed result
#[inline(always)]
pub fn raymarch_detailed(
    node: &SdfNode,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
    config: &RaymarchConfig,
) -> RaymarchResult {
    let dir = direction.normalize();
    let mut t = 0.0;
    let mut steps = 0;
    let mut stepper = RelaxedStepper::new(config);

    while t < max_distance && steps < config.max_steps {
        let point = origin + dir * t;
        let d = eval(node, point);

        if let Some(back) = stepper.overshoot(d) {
            t = back;
            steps += 1;
            continue;
        }

        if d.abs() < config.epsilon {
            return RaymarchResult {
                hit: true,
                distance: t,
                point,
                normal: normal(node, point, config.epsilon),
                steps,
            };
        }

        t = stepper.advance(t, d);
        steps += 1;
    }

    RaymarchResult {
        hit: false,
        distance: t,
        point: origin + dir * t,
        normal: Vec3::ZERO,
        steps,
    }
}

/// Batch raymarch (single-threaded)
pub fn raymarch_batch(node: &SdfNode, rays: &[Ray], max_distance: f32) -> Vec<Option<Hit>> {
    let config = RaymarchConfig::default().with_bound(eval_lipschitz(node));
    rays.iter()
        .map(|ray| raymarch_with_config(node, ray.origin, ray.direction, max_distance, &config))
        .collect()
}

/// Batch raymarch (parallel)
pub fn raymarch_batch_parallel(
    node: &SdfNode,
    rays: &[Ray],
    max_distance: f32,
) -> Vec<Option<Hit>> {
    let config = RaymarchConfig::default().with_bound(eval_lipschitz(node));
    rays.par_iter()
        .map(|ray| raymarch_with_config(node, ray.origin, ray.direction, max_distance, &config))
        .collect()
}

/// Render a depth buffer (Interpreter path)
#[allow(clippy::too_many_arguments)]
pub fn render_depth(
    node: &SdfNode,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<f32> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);

    let mut buffer = vec![0.0f32; width * height];
    let config = RaymarchConfig::default().with_bound(eval_lipschitz(node));

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            for (x, pixel) in row.iter_mut().enumerate() {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                *pixel = raymarch_with_config(node, camera_pos, ray_dir, max_distance, &config)
                    .map_or(f32::MAX, |hit| hit.distance);
            }
        });

    buffer
}

/// Render normals as RGB values (Interpreter path)
#[allow(clippy::too_many_arguments)]
pub fn render_normals(
    node: &SdfNode,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<[u8; 3]> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);

    let mut buffer = vec![[0u8; 3]; width * height];
    let config = RaymarchConfig::default().with_bound(eval_lipschitz(node));

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            for (x, pixel) in row.iter_mut().enumerate() {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                *pixel = raymarch_with_config(node, camera_pos, ray_dir, max_distance, &config)
                    .map_or([0, 0, 0], |hit| {
                        [
                            (hit.normal.x.mul_add(0.5, 0.5) * 255.0) as u8,
                            (hit.normal.y.mul_add(0.5, 0.5) * 255.0) as u8,
                            (hit.normal.z.mul_add(0.5, 0.5) * 255.0) as u8,
                        ]
                    });
            }
        });

    buffer
}

// ============================================================================
// Compiled Backend (CompiledSdf VM — ~3x faster than interpreter)
// ============================================================================

/// Sphere tracing using CompiledSdf (cache-efficient VM)
///
/// # Deep Fried Optimization
///
/// Eliminates Arc pointer chasing and branch misprediction from tree traversal.
/// The CompiledSdf flat instruction array is much more cache-friendly.
#[inline(always)]
///
/// Steps by `d / L` with the bound the `CompiledSdf` recorded at compile time.
pub fn raymarch_compiled(
    sdf: &CompiledSdf,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
) -> Option<Hit> {
    raymarch_compiled_with_config(
        sdf,
        origin,
        direction,
        max_distance,
        &RaymarchConfig::default(),
    )
}

/// Sphere tracing with config using CompiledSdf
#[inline(always)]
pub fn raymarch_compiled_with_config(
    sdf: &CompiledSdf,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
    config: &RaymarchConfig,
) -> Option<Hit> {
    let dir = direction.normalize();
    let mut t = 0.0;
    let mut steps = 0;
    let config = &config.with_bound(sdf.lipschitz);
    let mut stepper = RelaxedStepper::new(config);

    while t < max_distance && steps < config.max_steps {
        let point = origin + dir * t;
        let d = eval_compiled(sdf, point);

        if let Some(back) = stepper.overshoot(d) {
            t = back;
            steps += 1;
            continue;
        }

        if d.abs() < config.epsilon {
            return Some(Hit {
                distance: t,
                point,
                normal: eval_compiled_normal(sdf, point, config.epsilon),
                steps,
            });
        }

        t = stepper.advance(t, d);
        steps += 1;
    }

    None
}

/// Batch raymarch using CompiledSdf (parallel)
pub fn raymarch_compiled_batch_parallel(
    sdf: &CompiledSdf,
    rays: &[Ray],
    max_distance: f32,
) -> Vec<Option<Hit>> {
    rays.par_iter()
        .map(|ray| raymarch_compiled(sdf, ray.origin, ray.direction, max_distance))
        .collect()
}

/// Render depth buffer using CompiledSdf
#[allow(clippy::too_many_arguments)]
pub fn render_depth_compiled(
    sdf: &CompiledSdf,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<f32> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);

    let mut buffer = vec![0.0f32; width * height];

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            for (x, pixel) in row.iter_mut().enumerate() {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                *pixel = raymarch_compiled(sdf, camera_pos, ray_dir, max_distance)
                    .map_or(f32::MAX, |hit| hit.distance);
            }
        });

    buffer
}

/// Render normals using CompiledSdf
#[allow(clippy::too_many_arguments)]
pub fn render_normals_compiled(
    sdf: &CompiledSdf,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<[u8; 3]> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);

    let mut buffer = vec![[0u8; 3]; width * height];
    let config = RaymarchConfig::default();

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            for (x, pixel) in row.iter_mut().enumerate() {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                *pixel =
                    raymarch_compiled_with_config(sdf, camera_pos, ray_dir, max_distance, &config)
                        .map_or([0, 0, 0], |hit| {
                            [
                                (hit.normal.x.mul_add(0.5, 0.5) * 255.0) as u8,
                                (hit.normal.y.mul_add(0.5, 0.5) * 255.0) as u8,
                                (hit.normal.z.mul_add(0.5, 0.5) * 255.0) as u8,
                            ]
                        });
            }
        });

    buffer
}

// ============================================================================
// SIMD Packet Tracing (CompiledSdf + Vec3x8 — ~6x faster)
// ============================================================================

/// March 8 rays simultaneously using SIMD CompiledSdf evaluation
///
/// # Deep Fried Optimization
///
/// - 8 rays evaluated in parallel via `eval_compiled_simd()`
/// - Masked control flow: finished lanes keep their final distance
/// - No branch misprediction from per-lane hit/miss decisions
///
/// Returns array of 8 `Option<(f32, Vec3)>` (distance, hit_point) — normals
/// should be computed separately for hit lanes only to avoid waste.
pub fn raymarch_simd_8(
    sdf: &CompiledSdf,
    origins: Vec3x8,
    directions: Vec3x8,
    max_distance: f32,
    config: &RaymarchConfig,
) -> [Option<(f32, Vec3, u32)>; 8] {
    let eps = f32x8::splat(config.epsilon);
    let max_d = f32x8::splat(max_distance);
    let zero = f32x8::splat(0.0);
    // step = d / L with the compiled bound (see `RaymarchConfig::with_bound`)
    let inv_lip = f32x8::splat(1.0 / config.with_bound(sdf.lipschitz).lipschitz.max(1.0));

    let mut t = f32x8::splat(0.0);

    // finished: 0.0 = active, 1.0 = done
    let mut finished = f32x8::splat(0.0);
    // hit: 0.0 = miss, 1.0 = hit
    let mut hit_flags = f32x8::splat(0.0);
    let mut step_counts = [0u32; 8];

    for step in 0..config.max_steps {
        // p = origin + direction * t
        let px = origins.x + directions.x * t;
        let py = origins.y + directions.y * t;
        let pz = origins.z + directions.z * t;

        let points = Vec3x8 {
            x: px,
            y: py,
            z: pz,
        };
        let d = eval_compiled_simd(sdf, points);

        // Check hit: |d| < epsilon
        let abs_d = d.abs();
        let is_hit = abs_d.cmp_lt(eps);

        // Check too far: t > max_distance
        let is_far = t.cmp_gt(max_d);

        // Lanes that just became done
        let one = f32x8::splat(1.0);
        let newly_hit = is_hit.blend(one, zero) * (one - finished);
        let newly_far = is_far.blend(one, zero) * (one - finished);

        // Update hit flags (only for newly hit, not far)
        hit_flags += newly_hit;

        // Mark finished lanes
        finished = finished + newly_hit + newly_far;
        // Clamp finished to 1.0
        finished = finished.fast_min(one);

        // Update step counts for active lanes
        let finished_arr: [f32; 8] = finished.into();
        for i in 0..8 {
            if finished_arr[i] < 0.5 {
                step_counts[i] = step + 1;
            }
        }

        // Early exit if all lanes done
        if finished.cmp_ge(one).all() {
            break;
        }

        // Advance t for active lanes only (finished lanes get zero advance)
        let active = one - finished;
        t += d * inv_lip * active;
    }

    // Extract results
    let t_arr: [f32; 8] = t.into();
    let hit_arr: [f32; 8] = hit_flags.into();
    let (ox, oy, oz) = origins.to_array();
    let (dx, dy, dz) = directions.to_array();

    let mut results: [Option<(f32, Vec3, u32)>; 8] = [None; 8];

    for i in 0..8 {
        if hit_arr[i] > 0.5 {
            let p = Vec3::new(
                dx[i].mul_add(t_arr[i], ox[i]),
                dy[i].mul_add(t_arr[i], oy[i]),
                dz[i].mul_add(t_arr[i], oz[i]),
            );
            let _n = eval_compiled_normal(sdf, p, config.epsilon);
            results[i] = Some((t_arr[i], p, step_counts[i]));
        }
    }

    results
}

/// Render depth buffer using SIMD 8-ray packet tracing (Tile-Based)
///
/// # Deep Fried Optimization
///
/// Processes 8 horizontal pixels at once, maximizing SIMD utilization.
/// Each tile evaluates 8 points per step through `eval_compiled_simd()`.
#[allow(clippy::too_many_arguments)]
pub fn render_depth_compiled_simd(
    sdf: &CompiledSdf,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<f32> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);
    let config = RaymarchConfig::default();

    let mut buffer = vec![0.0f32; width * height];

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            // Process in 8-pixel tiles
            let mut x = 0;
            while x + 8 <= width {
                // Build 8 ray directions
                let mut dir_x = [0.0f32; 8];
                let mut dir_y = [0.0f32; 8];
                let mut dir_z = [0.0f32; 8];

                for i in 0..8 {
                    let u = ((x + i) as f32 * inv_width).mul_add(2.0, -1.0);
                    let rd = (row_vec + right_scaled * u).normalize();
                    dir_x[i] = rd.x;
                    dir_y[i] = rd.y;
                    dir_z[i] = rd.z;
                }

                let origins = Vec3x8::splat(camera_pos);
                let directions = Vec3x8::new(dir_x, dir_y, dir_z);

                let results = raymarch_simd_8(sdf, origins, directions, max_distance, &config);

                for i in 0..8 {
                    row[x + i] = match results[i] {
                        Some((dist, _, _)) => dist,
                        None => f32::MAX,
                    };
                }

                x += 8;
            }

            // Scalar fallback for remaining pixels
            while x < width {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                row[x] = raymarch_compiled(sdf, camera_pos, ray_dir, max_distance)
                    .map_or(f32::MAX, |hit| hit.distance);
                x += 1;
            }
        });

    buffer
}

// ============================================================================
// JIT Scalar Backend (Native code — ~5x faster)
// ============================================================================

/// Sphere tracing using JIT-compiled native code
///
/// # Deep Fried Optimization
///
/// Calls the JIT function pointer directly — zero interpreter overhead,
/// zero tree traversal, just a `call` instruction to native machine code.
#[cfg(feature = "jit")]
#[inline(always)]
pub fn raymarch_jit(
    sdf: &crate::compiled::jit::JitCompiledSdf,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
) -> Option<Hit> {
    raymarch_jit_with_config(
        sdf,
        origin,
        direction,
        max_distance,
        &RaymarchConfig::default(),
    )
}

/// Sphere tracing with config using JIT
#[cfg(feature = "jit")]
#[inline(always)]
pub fn raymarch_jit_with_config(
    sdf: &crate::compiled::jit::JitCompiledSdf,
    origin: Vec3,
    direction: Vec3,
    max_distance: f32,
    config: &RaymarchConfig,
) -> Option<Hit> {
    let dir = direction.normalize();
    let mut t = 0.0;
    let mut steps = 0;
    let config = &config.with_bound(sdf.lipschitz());
    let mut stepper = RelaxedStepper::new(config);

    while t < max_distance && steps < config.max_steps {
        let point = origin + dir * t;
        let d = sdf.eval(point);

        if let Some(back) = stepper.overshoot(d) {
            t = back;
            steps += 1;
            continue;
        }

        if d.abs() < config.epsilon {
            return Some(Hit {
                distance: t,
                point,
                normal: calc_normal_jit(sdf, point, config.epsilon),
                steps,
            });
        }

        t = stepper.advance(t, d);
        steps += 1;
    }

    None
}

/// Calculate surface normal using JIT (tetrahedral method)
#[cfg(feature = "jit")]
#[inline(always)]
fn calc_normal_jit(sdf: &crate::compiled::jit::JitCompiledSdf, p: Vec3, eps: f32) -> Vec3 {
    let k0 = Vec3::new(1.0, -1.0, -1.0);
    let k1 = Vec3::new(-1.0, -1.0, 1.0);
    let k2 = Vec3::new(-1.0, 1.0, -1.0);
    let k3 = Vec3::new(1.0, 1.0, 1.0);

    (k0 * sdf.eval(p + k0 * eps)
        + k1 * sdf.eval(p + k1 * eps)
        + k2 * sdf.eval(p + k2 * eps)
        + k3 * sdf.eval(p + k3 * eps))
    .normalize()
}

/// Batch raymarch using JIT (parallel)
#[cfg(feature = "jit")]
pub fn raymarch_jit_batch_parallel(
    sdf: &crate::compiled::jit::JitCompiledSdf,
    rays: &[Ray],
    max_distance: f32,
) -> Vec<Option<Hit>> {
    rays.par_iter()
        .map(|ray| raymarch_jit(sdf, ray.origin, ray.direction, max_distance))
        .collect()
}

/// Render depth buffer using JIT scalar
#[cfg(feature = "jit")]
#[allow(clippy::too_many_arguments)]
pub fn render_depth_jit(
    sdf: &crate::compiled::jit::JitCompiledSdf,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<f32> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);

    let mut buffer = vec![0.0f32; width * height];

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            for (x, pixel) in row.iter_mut().enumerate() {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                *pixel = raymarch_jit(sdf, camera_pos, ray_dir, max_distance)
                    .map_or(f32::MAX, |hit| hit.distance);
            }
        });

    buffer
}

// ============================================================================
// JIT SIMD Packet Tracing (8 rays via JitSimdSdf — ~15-20x faster)
// ============================================================================

/// March 8 rays simultaneously using JIT SIMD native code
///
/// # Deep Fried Optimization
///
/// - Calls `JitSimdSdf::eval_8()` which executes native AVX2/NEON instructions
/// - 8 SDF evaluations per function call
/// - Masked control flow: no branch penalty for mixed hit/miss lanes
/// - Combined with tile-based rendering for cache efficiency
#[cfg(feature = "jit")]
pub fn raymarch_jit_simd_8(
    sdf: &crate::compiled::jit::JitSimdSdf,
    compiled: &CompiledSdf,
    origins: Vec3x8,
    directions: Vec3x8,
    max_distance: f32,
    config: &RaymarchConfig,
) -> [Option<(f32, Vec3, u32)>; 8] {
    let eps = f32x8::splat(config.epsilon);
    let max_d = f32x8::splat(max_distance);
    let zero = f32x8::splat(0.0);
    let one = f32x8::splat(1.0);
    // step = d / L with the compiled bound (see `RaymarchConfig::with_bound`)
    let inv_lip = f32x8::splat(1.0 / config.with_bound(compiled.lipschitz).lipschitz.max(1.0));

    let mut t = f32x8::splat(0.0);
    let mut finished = f32x8::splat(0.0);
    let mut hit_flags = f32x8::splat(0.0);
    let mut step_counts = [0u32; 8];

    for step in 0..config.max_steps {
        // p = origin + direction * t
        let px = origins.x + directions.x * t;
        let py = origins.y + directions.y * t;
        let pz = origins.z + directions.z * t;

        // JIT SIMD eval: 8 SDF evaluations in one native call
        let px_arr: [f32; 8] = px.into();
        let py_arr: [f32; 8] = py.into();
        let pz_arr: [f32; 8] = pz.into();
        // SAFETY: `eval_8` calls a JIT-compiled native function pointer that was
        // generated from a validated CompiledSdf. The input arrays are stack-
        // allocated [f32; 8] with valid data. The JIT code only reads from these
        // arrays and writes to its return value, with no aliasing concerns.
        let d_arr = unsafe { sdf.eval_8(&px_arr, &py_arr, &pz_arr) };
        let d = f32x8::from(d_arr);

        // Check hit: |d| < epsilon
        let abs_d = d.abs();
        let is_hit = abs_d.cmp_lt(eps);
        let is_far = t.cmp_gt(max_d);

        let newly_hit = is_hit.blend(one, zero) * (one - finished);
        let newly_far = is_far.blend(one, zero) * (one - finished);

        hit_flags += newly_hit;
        finished = (finished + newly_hit + newly_far).fast_min(one);

        let finished_arr: [f32; 8] = finished.into();
        for i in 0..8 {
            if finished_arr[i] < 0.5 {
                step_counts[i] = step + 1;
            }
        }

        if finished.cmp_ge(one).all() {
            break;
        }

        let active = one - finished;
        t += d * inv_lip * active;
    }

    // Extract results
    let t_arr: [f32; 8] = t.into();
    let hit_arr: [f32; 8] = hit_flags.into();
    let (ox, oy, oz) = origins.to_array();
    let (dx, dy, dz) = directions.to_array();

    let mut results: [Option<(f32, Vec3, u32)>; 8] = [None; 8];

    for i in 0..8 {
        if hit_arr[i] > 0.5 {
            let p = Vec3::new(
                dx[i].mul_add(t_arr[i], ox[i]),
                dy[i].mul_add(t_arr[i], oy[i]),
                dz[i].mul_add(t_arr[i], oz[i]),
            );
            // Normal from compiled path (no SIMD normal helper needed)
            results[i] = Some((t_arr[i], p, step_counts[i]));
        }
    }

    results
}

/// Render depth buffer using JIT SIMD (Tile-Based, ~15-20x faster)
///
/// # Deep Fried Optimization
///
/// Combines JIT SIMD packet tracing with tile-based rendering:
/// - 8 horizontal pixels form one SIMD packet
/// - `JitSimdSdf::eval_8()` evaluates all 8 in one native call
/// - Row-based parallelization via rayon
/// - Scalar fallback for edge pixels (width % 8 != 0)
#[cfg(feature = "jit")]
#[allow(clippy::too_many_arguments)]
pub fn render_depth_jit_simd(
    sdf: &crate::compiled::jit::JitSimdSdf,
    compiled: &CompiledSdf,
    camera_pos: Vec3,
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
    max_distance: f32,
) -> Vec<f32> {
    let (forward, right, up, half_width, half_height) =
        camera_basis(camera_dir, camera_up, width, height, fov);
    let config = RaymarchConfig::default();

    let mut buffer = vec![0.0f32; width * height];

    buffer
        .par_chunks_mut(width)
        .enumerate()
        .for_each(|(y, row)| {
            let v = (y as f32 / height as f32).mul_add(-2.0, 1.0);
            let row_vec = forward + up * (v * half_height);
            let right_scaled = right * half_width;
            let inv_width = 1.0 / width as f32;

            let mut x = 0;
            while x + 8 <= width {
                let mut dir_x = [0.0f32; 8];
                let mut dir_y = [0.0f32; 8];
                let mut dir_z = [0.0f32; 8];

                for i in 0..8 {
                    let u = ((x + i) as f32 * inv_width).mul_add(2.0, -1.0);
                    let rd = (row_vec + right_scaled * u).normalize();
                    dir_x[i] = rd.x;
                    dir_y[i] = rd.y;
                    dir_z[i] = rd.z;
                }

                let origins = Vec3x8::splat(camera_pos);
                let directions = Vec3x8::new(dir_x, dir_y, dir_z);

                let results =
                    raymarch_jit_simd_8(sdf, compiled, origins, directions, max_distance, &config);

                for i in 0..8 {
                    row[x + i] = match results[i] {
                        Some((dist, _, _)) => dist,
                        None => f32::MAX,
                    };
                }

                x += 8;
            }

            // Scalar fallback for edge pixels
            while x < width {
                let u = (x as f32 * inv_width).mul_add(2.0, -1.0);
                let ray_dir = (row_vec + right_scaled * u).normalize();

                row[x] = raymarch_compiled(compiled, camera_pos, ray_dir, max_distance)
                    .map_or(f32::MAX, |hit| hit.distance);
                x += 1;
            }
        });

    buffer
}

// ============================================================================
// Shared Utilities
// ============================================================================

/// Pre-calculate camera basis vectors (shared by all render functions)
///
/// Returns (forward, right, up, half_width, half_height)
#[inline(always)]
fn camera_basis(
    camera_dir: Vec3,
    camera_up: Vec3,
    width: usize,
    height: usize,
    fov: f32,
) -> (Vec3, Vec3, Vec3, f32, f32) {
    let forward = camera_dir.normalize();
    let right = forward.cross(camera_up).normalize();
    let up = right.cross(forward);

    let aspect = width as f32 / height as f32;
    let half_height = (fov * 0.5).tan();
    let half_width = half_height * aspect;

    (forward, right, up, half_width, half_height)
}

#[cfg(test)]
mod tests {
    use super::*;

    // ============ Interpreter Tests ============

    #[test]
    fn test_raymarch_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch(&sphere, origin, direction, 10.0);
        assert!(hit.is_some());

        let hit = hit.unwrap();
        assert!((hit.distance - 4.0).abs() < 0.01);
        assert!((hit.normal.x - (-1.0)).abs() < 0.01);
    }

    #[test]
    fn test_raymarch_miss() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 5.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch(&sphere, origin, direction, 20.0);
        assert!(hit.is_none());
    }

    #[test]
    fn test_raymarch_config_quality() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let config_fast = RaymarchConfig::fast();
        let config_hq = RaymarchConfig::high_quality();

        let hit_fast = raymarch_with_config(&sphere, origin, direction, 10.0, &config_fast);
        let hit_hq = raymarch_with_config(&sphere, origin, direction, 10.0, &config_hq);

        assert!(hit_fast.is_some());
        assert!(hit_hq.is_some());
    }

    #[test]
    fn test_raymarch_detailed() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let config = RaymarchConfig::default();
        let result = raymarch_detailed(&sphere, origin, direction, 10.0, &config);

        assert!(result.hit);
        assert!(result.steps > 0);
        assert!(result.steps < config.max_steps);
    }

    #[test]
    fn test_render_depth() {
        let sphere = SdfNode::sphere(1.0);
        let camera_pos = Vec3::new(0.0, 0.0, -5.0);
        let camera_dir = Vec3::new(0.0, 0.0, 1.0);
        let camera_up = Vec3::new(0.0, 1.0, 0.0);

        let depth = render_depth(
            &sphere,
            camera_pos,
            camera_dir,
            camera_up,
            8,
            8,
            std::f32::consts::FRAC_PI_4,
            10.0,
        );

        assert_eq!(depth.len(), 64);
        let center_depth = depth[8 * 4 + 4];
        assert!(center_depth < f32::MAX);
        assert!((center_depth - 4.0).abs() < 0.5);
    }

    #[test]
    fn test_render_normals() {
        let sphere = SdfNode::sphere(1.0);
        let camera_pos = Vec3::new(0.0, 0.0, -5.0);
        let camera_dir = Vec3::new(0.0, 0.0, 1.0);
        let camera_up = Vec3::new(0.0, 1.0, 0.0);

        let normals = render_normals(
            &sphere,
            camera_pos,
            camera_dir,
            camera_up,
            8,
            8,
            std::f32::consts::FRAC_PI_4,
            10.0,
        );

        assert_eq!(normals.len(), 64);
        let center_normal = normals[8 * 4 + 4];
        assert!(center_normal[2] < 64);
    }

    #[test]
    fn test_over_relaxation() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let config = RaymarchConfig::fast();
        let hit = raymarch_with_config(&sphere, origin, direction, 10.0, &config);

        assert!(hit.is_some());
        let hit = hit.unwrap();
        assert!((hit.distance - 4.0).abs() < 0.1);
    }

    // ============ Relaxed Sphere Tracing Tests ============

    #[test]
    fn test_relaxed_sphere_hit() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch_relaxed(&sphere, origin, direction, 10.0);
        assert!(hit.is_some());
        let hit = hit.unwrap();
        assert!(
            (hit.distance - 4.0).abs() < 0.01,
            "Relaxed hit distance: {}",
            hit.distance
        );
    }

    #[test]
    fn test_relaxed_sphere_miss() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 5.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch_relaxed(&sphere, origin, direction, 20.0);
        assert!(hit.is_none());
    }

    #[test]
    fn test_relaxed_fewer_steps() {
        let sphere = SdfNode::sphere(1.0);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        // Standard tracing
        let config_std = RaymarchConfig::default();
        let result_std = raymarch_detailed(&sphere, origin, direction, 10.0, &config_std);

        // Relaxed tracing
        let config_rel = RaymarchConfig::relaxed(&sphere);
        let result_rel = raymarch_detailed(&sphere, origin, direction, 10.0, &config_rel);

        assert!(result_std.hit);
        assert!(result_rel.hit);
        assert!((result_rel.distance - result_std.distance).abs() < 1e-3);
        // Head-on the relaxed ray overshoots once and retreats (+1 step); the
        // saving shows at grazing incidence, where plain tracing crawls.
        let plane = SdfNode::plane(Vec3::Y, 0.0);
        let origin = Vec3::new(-3.0, 0.5, 0.0);
        let (s, c) = 6.0f32.to_radians().sin_cos();
        let direction = Vec3::new(c, -s, 0.0);
        let std = raymarch_detailed(&plane, origin, direction, 20.0, &config_std);
        let rel = raymarch_detailed(
            &plane,
            origin,
            direction,
            20.0,
            &RaymarchConfig::relaxed(&plane),
        );
        assert!(std.hit && rel.hit);
        assert!((std.distance - rel.distance).abs() < 1e-3);
        assert!(
            rel.steps < std.steps,
            "Relaxed steps={} should be < standard steps={} at grazing incidence",
            rel.steps,
            std.steps
        );
        // Hit distance should be very similar
        assert!(
            (result_std.distance - result_rel.distance).abs() < 0.01,
            "Standard={}, Relaxed={}",
            result_std.distance,
            result_rel.distance
        );
    }

    #[test]
    fn test_relaxed_complex_scene() {
        let scene = SdfNode::sphere(1.0)
            .translate(1.0, 0.0, 0.0)
            .smooth_union(SdfNode::box3d(0.5, 0.5, 0.5).translate(-1.0, 0.0, 0.0), 0.3);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit_std = raymarch(&scene, origin, direction, 20.0);
        let hit_rel = raymarch_relaxed(&scene, origin, direction, 20.0);

        assert!(hit_std.is_some());
        assert!(hit_rel.is_some());
        assert!((hit_std.unwrap().distance - hit_rel.unwrap().distance).abs() < 0.05);
    }

    #[test]
    fn test_relaxed_config_lipschitz() {
        let sphere = SdfNode::sphere(1.0);
        let config = RaymarchConfig::relaxed(&sphere);
        assert_eq!(config.lipschitz, 1.0);
        assert!((config.omega - 1.6).abs() < 1e-5);

        let twisted = SdfNode::Twist {
            child: std::sync::Arc::new(SdfNode::sphere(1.0)),
            strength: 0.5,
        };
        let config_twist = RaymarchConfig::relaxed(&twisted);
        assert!(
            config_twist.lipschitz > 1.0,
            "Twisted SDF should have L > 1, got {}",
            config_twist.lipschitz
        );
    }

    // ============ Compiled Backend Tests ============

    #[test]
    fn test_raymarch_compiled_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&sphere);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch_compiled(&compiled, origin, direction, 10.0);
        assert!(hit.is_some());

        let hit = hit.unwrap();
        assert!(
            (hit.distance - 4.0).abs() < 0.01,
            "Compiled sphere hit distance: {}",
            hit.distance
        );
    }

    #[test]
    fn test_raymarch_compiled_miss() {
        let sphere = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&sphere);
        let origin = Vec3::new(-5.0, 5.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch_compiled(&compiled, origin, direction, 20.0);
        assert!(hit.is_none());
    }

    #[test]
    fn test_raymarch_compiled_matches_interpreter() {
        let shape = SdfNode::sphere(1.0)
            .smooth_union(SdfNode::box3d(0.5, 0.5, 0.5).translate(2.0, 0.0, 0.0), 0.2);
        let compiled = CompiledSdf::compile(&shape);
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit_interp = raymarch(&shape, origin, direction, 20.0);
        let hit_compiled = raymarch_compiled(&compiled, origin, direction, 20.0);

        assert_eq!(hit_interp.is_some(), hit_compiled.is_some());
        if let (Some(h1), Some(h2)) = (hit_interp, hit_compiled) {
            assert!(
                (h1.distance - h2.distance).abs() < 0.01,
                "Interpreter={}, Compiled={}",
                h1.distance,
                h2.distance
            );
        }
    }

    #[test]
    fn test_render_depth_compiled() {
        let sphere = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&sphere);
        let camera_pos = Vec3::new(0.0, 0.0, -5.0);
        let camera_dir = Vec3::new(0.0, 0.0, 1.0);
        let camera_up = Vec3::new(0.0, 1.0, 0.0);

        let depth = render_depth_compiled(
            &compiled,
            camera_pos,
            camera_dir,
            camera_up,
            8,
            8,
            std::f32::consts::FRAC_PI_4,
            10.0,
        );

        assert_eq!(depth.len(), 64);
        let center_depth = depth[8 * 4 + 4];
        assert!(center_depth < f32::MAX);
        assert!(
            (center_depth - 4.0).abs() < 0.5,
            "Compiled depth center: {}",
            center_depth
        );
    }

    // ============ SIMD Packet Tests ============

    #[test]
    fn test_raymarch_simd_8_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&sphere);
        let config = RaymarchConfig::default();

        // 8 rays from different X positions, all heading +Z
        let origins = Vec3x8::new(
            [-2.0, -1.0, 0.0, 1.0, 2.0, 0.0, 0.0, 0.0],
            [0.0; 8],
            [-5.0; 8],
        );
        let directions = Vec3x8::new([0.0; 8], [0.0; 8], [1.0; 8]);

        let results = raymarch_simd_8(&compiled, origins, directions, 20.0, &config);

        // Center ray (index 2) should hit at ~4.0
        assert!(results[2].is_some(), "Center ray should hit sphere");
        let (dist, _, _) = results[2].unwrap();
        assert!((dist - 4.0).abs() < 0.1, "Center ray dist: {}", dist);

        // Edge rays (index 0, 4) should miss (|x| = 2 > radius 1)
        assert!(results[0].is_none(), "Ray at x=-2 should miss");
        assert!(results[4].is_none(), "Ray at x=2 should miss");
    }

    #[test]
    fn test_render_depth_compiled_simd() {
        let sphere = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&sphere);
        let camera_pos = Vec3::new(0.0, 0.0, -5.0);
        let camera_dir = Vec3::new(0.0, 0.0, 1.0);
        let camera_up = Vec3::new(0.0, 1.0, 0.0);

        let depth_scalar = render_depth_compiled(
            &compiled,
            camera_pos,
            camera_dir,
            camera_up,
            16,
            16,
            std::f32::consts::FRAC_PI_4,
            10.0,
        );
        let depth_simd = render_depth_compiled_simd(
            &compiled,
            camera_pos,
            camera_dir,
            camera_up,
            16,
            16,
            std::f32::consts::FRAC_PI_4,
            10.0,
        );

        assert_eq!(depth_scalar.len(), depth_simd.len());

        // Compare results: should be very close
        for i in 0..depth_scalar.len() {
            if depth_scalar[i] < f32::MAX && depth_simd[i] < f32::MAX {
                assert!(
                    (depth_scalar[i] - depth_simd[i]).abs() < 0.5,
                    "Pixel {} mismatch: scalar={}, simd={}",
                    i,
                    depth_scalar[i],
                    depth_simd[i]
                );
            }
        }
    }

    // ============ JIT Backend Tests ============

    #[cfg(feature = "jit")]
    #[test]
    fn test_raymarch_jit_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let jit = crate::compiled::jit::JitCompiledSdf::compile(&sphere).unwrap();
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit = raymarch_jit(&jit, origin, direction, 10.0);
        assert!(hit.is_some());

        let hit = hit.unwrap();
        assert!(
            (hit.distance - 4.0).abs() < 0.01,
            "JIT sphere hit distance: {}",
            hit.distance
        );
    }

    #[cfg(feature = "jit")]
    #[test]
    fn test_raymarch_jit_matches_compiled() {
        let shape =
            SdfNode::sphere(1.0).union(SdfNode::box3d(0.5, 0.5, 0.5).translate(2.0, 0.0, 0.0));
        let compiled = CompiledSdf::compile(&shape);
        let jit = crate::compiled::jit::JitCompiledSdf::compile(&shape).unwrap();
        let origin = Vec3::new(-5.0, 0.0, 0.0);
        let direction = Vec3::new(1.0, 0.0, 0.0);

        let hit_compiled = raymarch_compiled(&compiled, origin, direction, 20.0);
        let hit_jit = raymarch_jit(&jit, origin, direction, 20.0);

        assert_eq!(hit_compiled.is_some(), hit_jit.is_some());
        if let (Some(h1), Some(h2)) = (hit_compiled, hit_jit) {
            assert!(
                (h1.distance - h2.distance).abs() < 0.01,
                "Compiled={}, JIT={}",
                h1.distance,
                h2.distance
            );
        }
    }

    #[cfg(feature = "jit")]
    #[test]
    fn test_raymarch_jit_simd_8_sphere() {
        let sphere = SdfNode::sphere(1.0);
        let compiled = CompiledSdf::compile(&sphere);
        let jit_simd = crate::compiled::jit::JitSimdSdf::compile(&compiled).unwrap();
        let config = RaymarchConfig::default();

        let origins = Vec3x8::new([0.0; 8], [0.0; 8], [-5.0; 8]);
        let directions = Vec3x8::new([0.0; 8], [0.0; 8], [1.0; 8]);

        let results = raymarch_jit_simd_8(&jit_simd, &compiled, origins, directions, 20.0, &config);

        // All 8 rays from origin heading +Z should hit sphere at ~4.0
        for (i, result) in results.iter().enumerate() {
            assert!(result.is_some(), "Ray {} should hit", i);
            let (dist, _, _) = result.unwrap();
            assert!((dist - 4.0).abs() < 0.1, "Ray {} dist: {}", i, dist);
        }
    }
}