vowl-wasm 0.2.0

//! SIMD-accelerated force calculations for WASM
//!
//! This module provides SIMD128 vectorized implementations of force calculations
//! for improved performance on supported platforms. SIMD operations process
//! multiple data points simultaneously, achieving significant speedups.
//!
//! ## Usage
//! Enable the `simd` feature flag and rebuild:
//! ```bash
//! wasm-pack build --target web --release --features simd
//! ```
//!
//! ## Performance
//! - Batch repulsion: ~2-4x speedup with 4-way f32 vectorization
//! - Center force: ~2x speedup with batch processing
//! - Barnes-Hut leaf: ~1.5x speedup for exact calculations
//!
//! ## Safety
//! SIMD functions are marked unsafe because they use WASM SIMD128 intrinsics
//! that require the target_feature to be enabled. The caller must ensure
//! SIMD is supported before calling these functions.

// Allow unsafe code in this module for SIMD intrinsics
#![allow(unsafe_code)]

#[cfg(all(target_arch = "wasm32", feature = "simd"))]
use core::arch::wasm32::*;

use nalgebra::Vector2;

/// Minimum distance squared to prevent division by zero
const MIN_DISTANCE_SQ: f32 = 0.0001;

/// SIMD-optimized position data for batch processing
/// Uses f32 for SIMD compatibility (WASM SIMD128 = 4x f32)
#[derive(Debug, Clone)]
pub struct SimdPositions {
    /// X coordinates (aligned for SIMD)
    pub x: Vec<f32>,
    /// Y coordinates (aligned for SIMD)
    pub y: Vec<f32>,
    /// Node masses/weights
    pub mass: Vec<f32>,
    /// Number of positions
    pub len: usize,
}

impl SimdPositions {
    /// Create from f64 position pairs
    pub fn from_positions(positions: &[(f64, f64)]) -> Self {
        let len = positions.len();
        let mut x = Vec::with_capacity(len);
        let mut y = Vec::with_capacity(len);
        let mut mass = Vec::with_capacity(len);

        for (px, py) in positions {
            x.push(*px as f32);
            y.push(*py as f32);
            mass.push(1.0f32);
        }

        Self { x, y, mass, len }
    }

    /// Create from positions with weights
    pub fn from_positions_weighted(positions: &[(f64, f64, f64)]) -> Self {
        let len = positions.len();
        let mut x = Vec::with_capacity(len);
        let mut y = Vec::with_capacity(len);
        let mut mass = Vec::with_capacity(len);

        for (px, py, w) in positions {
            x.push(*px as f32);
            y.push(*py as f32);
            mass.push(*w as f32);
        }

        Self { x, y, mass, len }
    }
}

/// SIMD-optimized force accumulator
#[derive(Debug, Clone)]
pub struct SimdForces {
    /// X force components
    pub fx: Vec<f32>,
    /// Y force components
    pub fy: Vec<f32>,
    /// Number of forces
    pub len: usize,
}

impl SimdForces {
    /// Create zero-initialized forces
    pub fn zeros(len: usize) -> Self {
        Self {
            fx: vec![0.0f32; len],
            fy: vec![0.0f32; len],
            len,
        }
    }

    /// Convert to a `Vector2<f64>` array
    pub fn to_vectors(&self) -> Vec<Vector2<f64>> {
        (0..self.len)
            .map(|i| Vector2::new(self.fx[i] as f64, self.fy[i] as f64))
            .collect()
    }

    /// Add another force accumulator
    pub fn add(&mut self, other: &SimdForces) {
        for i in 0..self.len.min(other.len) {
            self.fx[i] += other.fx[i];
            self.fy[i] += other.fy[i];
        }
    }
}

/// Calculate batch repulsion forces using SIMD
///
/// This function calculates repulsion forces for all pairs of nodes
/// using SIMD vectorization for improved performance.
///
/// # Arguments
/// * `positions` - Node positions in SIMD format
/// * `strength` - Charge strength (negative for repulsion)
///
/// # Returns
/// Force accumulator with computed forces
#[cfg(all(target_arch = "wasm32", feature = "simd"))]
#[target_feature(enable = "simd128")]
pub unsafe fn calculate_repulsion_simd(positions: &SimdPositions, strength: f32) -> SimdForces {
    let mut forces = SimdForces::zeros(positions.len);
    let n = positions.len;

    // Process pairs (i, j) where j > i
    for i in 0..n {
        let xi = positions.x[i];
        let yi = positions.y[i];

        // Create SIMD vector for position i (broadcast)
        let xi_v = f32x4_splat(xi);
        let yi_v = f32x4_splat(yi);

        let mut fx_acc = f32x4_splat(0.0);
        let mut fy_acc = f32x4_splat(0.0);

        // Process 4 nodes at a time
        let chunks = (n - i - 1) / 4;
        for chunk in 0..chunks {
            let base_j = i + 1 + chunk * 4;

            // Load 4 positions at once
            let xj_v = f32x4(
                positions.x[base_j],
                positions.x[base_j + 1],
                positions.x[base_j + 2],
                positions.x[base_j + 3],
            );
            let yj_v = f32x4(
                positions.y[base_j],
                positions.y[base_j + 1],
                positions.y[base_j + 2],
                positions.y[base_j + 3],
            );

            // Calculate deltas
            let dx = f32x4_sub(xi_v, xj_v);
            let dy = f32x4_sub(yi_v, yj_v);

            // Calculate distance squared
            let dx_sq = f32x4_mul(dx, dx);
            let dy_sq = f32x4_mul(dy, dy);
            let dist_sq = f32x4_add(dx_sq, dy_sq);

            // Clamp to minimum distance
            let min_dist = f32x4_splat(MIN_DISTANCE_SQ);
            let dist_sq_clamped = f32x4_pmax(dist_sq, min_dist);

            // Calculate force magnitude: strength / dist_sq
            let strength_v = f32x4_splat(strength);
            let force_mag = f32x4_div(strength_v, dist_sq_clamped);

            // Calculate distance (for normalization)
            let dist = f32x4_sqrt(dist_sq_clamped);

            // Normalize and apply force
            let inv_dist = f32x4_div(f32x4_splat(1.0), dist);
            let fx = f32x4_mul(f32x4_mul(dx, inv_dist), force_mag);
            let fy = f32x4_mul(f32x4_mul(dy, inv_dist), force_mag);

            // Accumulate forces for node i
            fx_acc = f32x4_add(fx_acc, fx);
            fy_acc = f32x4_add(fy_acc, fy);

            // Apply Newton's third law: equal and opposite forces on j nodes
            forces.fx[base_j] -= f32x4_extract_lane::<0>(fx);
            forces.fx[base_j + 1] -= f32x4_extract_lane::<1>(fx);
            forces.fx[base_j + 2] -= f32x4_extract_lane::<2>(fx);
            forces.fx[base_j + 3] -= f32x4_extract_lane::<3>(fx);

            forces.fy[base_j] -= f32x4_extract_lane::<0>(fy);
            forces.fy[base_j + 1] -= f32x4_extract_lane::<1>(fy);
            forces.fy[base_j + 2] -= f32x4_extract_lane::<2>(fy);
            forces.fy[base_j + 3] -= f32x4_extract_lane::<3>(fy);
        }

        // Sum up accumulated forces for node i
        forces.fx[i] += f32x4_extract_lane::<0>(fx_acc)
            + f32x4_extract_lane::<1>(fx_acc)
            + f32x4_extract_lane::<2>(fx_acc)
            + f32x4_extract_lane::<3>(fx_acc);
        forces.fy[i] += f32x4_extract_lane::<0>(fy_acc)
            + f32x4_extract_lane::<1>(fy_acc)
            + f32x4_extract_lane::<2>(fy_acc)
            + f32x4_extract_lane::<3>(fy_acc);

        // Handle remaining nodes (not divisible by 4)
        let remainder_start = i + 1 + chunks * 4;
        for j in remainder_start..n {
            let dx = xi - positions.x[j];
            let dy = yi - positions.y[j];
            let dist_sq = (dx * dx + dy * dy).max(MIN_DISTANCE_SQ);
            let dist = dist_sq.sqrt();
            let force_mag = strength / dist_sq;
            let fx = (dx / dist) * force_mag;
            let fy = (dy / dist) * force_mag;

            forces.fx[i] += fx;
            forces.fy[i] += fy;
            forces.fx[j] -= fx;
            forces.fy[j] -= fy;
        }
    }

    forces
}

/// Non-SIMD fallback for batch repulsion
#[cfg(not(all(target_arch = "wasm32", feature = "simd")))]
pub fn calculate_repulsion_simd(positions: &SimdPositions, strength: f32) -> SimdForces {
    calculate_repulsion_scalar(positions, strength)
}

/// Scalar implementation of batch repulsion (fallback)
pub fn calculate_repulsion_scalar(positions: &SimdPositions, strength: f32) -> SimdForces {
    let mut forces = SimdForces::zeros(positions.len);
    let n = positions.len;

    for i in 0..n {
        let xi = positions.x[i];
        let yi = positions.y[i];

        for j in (i + 1)..n {
            let dx = xi - positions.x[j];
            let dy = yi - positions.y[j];
            let dist_sq = (dx * dx + dy * dy).max(MIN_DISTANCE_SQ);
            let dist = dist_sq.sqrt();
            let force_mag = strength / dist_sq;
            let fx = (dx / dist) * force_mag;
            let fy = (dy / dist) * force_mag;

            forces.fx[i] += fx;
            forces.fy[i] += fy;
            forces.fx[j] -= fx;
            forces.fy[j] -= fy;
        }
    }

    forces
}

/// Calculate centering forces using SIMD
///
/// Applies a force pulling all nodes toward a center point.
///
/// # Arguments
/// * `positions` - Node positions in SIMD format
/// * `center_x` - Center X coordinate
/// * `center_y` - Center Y coordinate
/// * `strength` - Center force strength
///
/// # Returns
/// Force accumulator with computed centering forces
#[cfg(all(target_arch = "wasm32", feature = "simd"))]
#[target_feature(enable = "simd128")]
pub unsafe fn calculate_center_force_simd(
    positions: &SimdPositions,
    center_x: f32,
    center_y: f32,
    strength: f32,
) -> SimdForces {
    let mut forces = SimdForces::zeros(positions.len);
    let n = positions.len;

    let cx_v = f32x4_splat(center_x);
    let cy_v = f32x4_splat(center_y);
    let strength_v = f32x4_splat(strength);

    // Process 4 nodes at a time
    let chunks = n / 4;
    for chunk in 0..chunks {
        let base = chunk * 4;

        // Load positions
        let x_v = f32x4(
            positions.x[base],
            positions.x[base + 1],
            positions.x[base + 2],
            positions.x[base + 3],
        );
        let y_v = f32x4(
            positions.y[base],
            positions.y[base + 1],
            positions.y[base + 2],
            positions.y[base + 3],
        );

        // Calculate force: (center - pos) * strength
        let fx = f32x4_mul(f32x4_sub(cx_v, x_v), strength_v);
        let fy = f32x4_mul(f32x4_sub(cy_v, y_v), strength_v);

        // Store results
        forces.fx[base] = f32x4_extract_lane::<0>(fx);
        forces.fx[base + 1] = f32x4_extract_lane::<1>(fx);
        forces.fx[base + 2] = f32x4_extract_lane::<2>(fx);
        forces.fx[base + 3] = f32x4_extract_lane::<3>(fx);

        forces.fy[base] = f32x4_extract_lane::<0>(fy);
        forces.fy[base + 1] = f32x4_extract_lane::<1>(fy);
        forces.fy[base + 2] = f32x4_extract_lane::<2>(fy);
        forces.fy[base + 3] = f32x4_extract_lane::<3>(fy);
    }

    // Handle remainder
    for i in (chunks * 4)..n {
        forces.fx[i] = (center_x - positions.x[i]) * strength;
        forces.fy[i] = (center_y - positions.y[i]) * strength;
    }

    forces
}

/// Non-SIMD fallback for center force
#[cfg(not(all(target_arch = "wasm32", feature = "simd")))]
pub fn calculate_center_force_simd(
    positions: &SimdPositions,
    center_x: f32,
    center_y: f32,
    strength: f32,
) -> SimdForces {
    calculate_center_force_scalar(positions, center_x, center_y, strength)
}

/// Scalar implementation of center force (fallback)
pub fn calculate_center_force_scalar(
    positions: &SimdPositions,
    center_x: f32,
    center_y: f32,
    strength: f32,
) -> SimdForces {
    let mut forces = SimdForces::zeros(positions.len);

    for i in 0..positions.len {
        forces.fx[i] = (center_x - positions.x[i]) * strength;
        forces.fy[i] = (center_y - positions.y[i]) * strength;
    }

    forces
}

/// Calculate exact leaf forces using SIMD (for Barnes-Hut)
///
/// Computes exact repulsion from multiple source nodes to a single target.
///
/// # Arguments
/// * `target_x` - Target node X coordinate
/// * `target_y` - Target node Y coordinate
/// * `sources` - Source node positions
/// * `strength` - Charge strength
///
/// # Returns
/// Total force vector on target node
#[cfg(all(target_arch = "wasm32", feature = "simd"))]
#[target_feature(enable = "simd128")]
pub unsafe fn calculate_leaf_force_simd(
    target_x: f32,
    target_y: f32,
    sources: &SimdPositions,
    strength: f32,
) -> (f32, f32) {
    let n = sources.len;
    if n == 0 {
        return (0.0, 0.0);
    }

    let tx_v = f32x4_splat(target_x);
    let ty_v = f32x4_splat(target_y);
    let strength_v = f32x4_splat(strength);
    let min_dist_v = f32x4_splat(MIN_DISTANCE_SQ);
    let one_v = f32x4_splat(1.0);

    let mut fx_acc = f32x4_splat(0.0);
    let mut fy_acc = f32x4_splat(0.0);

    let chunks = n / 4;
    for chunk in 0..chunks {
        let base = chunk * 4;

        // Load source positions
        let sx_v = f32x4(
            sources.x[base],
            sources.x[base + 1],
            sources.x[base + 2],
            sources.x[base + 3],
        );
        let sy_v = f32x4(
            sources.y[base],
            sources.y[base + 1],
            sources.y[base + 2],
            sources.y[base + 3],
        );
        let mass_v = f32x4(
            sources.mass[base],
            sources.mass[base + 1],
            sources.mass[base + 2],
            sources.mass[base + 3],
        );

        // Calculate delta: target - source
        let dx = f32x4_sub(tx_v, sx_v);
        let dy = f32x4_sub(ty_v, sy_v);

        // Distance squared
        let dist_sq = f32x4_add(f32x4_mul(dx, dx), f32x4_mul(dy, dy));
        let dist_sq_clamped = f32x4_pmax(dist_sq, min_dist_v);

        // Force magnitude: strength * mass / dist_sq
        let force_mag = f32x4_div(f32x4_mul(strength_v, mass_v), dist_sq_clamped);

        // Normalize
        let dist = f32x4_sqrt(dist_sq_clamped);
        let inv_dist = f32x4_div(one_v, dist);

        // Calculate force components
        let fx = f32x4_mul(f32x4_mul(dx, inv_dist), force_mag);
        let fy = f32x4_mul(f32x4_mul(dy, inv_dist), force_mag);

        fx_acc = f32x4_add(fx_acc, fx);
        fy_acc = f32x4_add(fy_acc, fy);
    }

    // Sum SIMD lanes
    let mut total_fx = f32x4_extract_lane::<0>(fx_acc)
        + f32x4_extract_lane::<1>(fx_acc)
        + f32x4_extract_lane::<2>(fx_acc)
        + f32x4_extract_lane::<3>(fx_acc);
    let mut total_fy = f32x4_extract_lane::<0>(fy_acc)
        + f32x4_extract_lane::<1>(fy_acc)
        + f32x4_extract_lane::<2>(fy_acc)
        + f32x4_extract_lane::<3>(fy_acc);

    // Handle remainder
    for i in (chunks * 4)..n {
        let dx = target_x - sources.x[i];
        let dy = target_y - sources.y[i];
        let dist_sq = (dx * dx + dy * dy).max(MIN_DISTANCE_SQ);
        let dist = dist_sq.sqrt();
        let force_mag = strength * sources.mass[i] / dist_sq;
        total_fx += (dx / dist) * force_mag;
        total_fy += (dy / dist) * force_mag;
    }

    (total_fx, total_fy)
}

/// Non-SIMD fallback for leaf force
#[cfg(not(all(target_arch = "wasm32", feature = "simd")))]
pub fn calculate_leaf_force_simd(
    target_x: f32,
    target_y: f32,
    sources: &SimdPositions,
    strength: f32,
) -> (f32, f32) {
    calculate_leaf_force_scalar(target_x, target_y, sources, strength)
}

/// Scalar implementation of leaf force (fallback)
pub fn calculate_leaf_force_scalar(
    target_x: f32,
    target_y: f32,
    sources: &SimdPositions,
    strength: f32,
) -> (f32, f32) {
    let mut fx = 0.0f32;
    let mut fy = 0.0f32;

    for i in 0..sources.len {
        let dx = target_x - sources.x[i];
        let dy = target_y - sources.y[i];
        let dist_sq = (dx * dx + dy * dy).max(MIN_DISTANCE_SQ);
        let dist = dist_sq.sqrt();
        let force_mag = strength * sources.mass[i] / dist_sq;
        fx += (dx / dist) * force_mag;
        fy += (dy / dist) * force_mag;
    }

    (fx, fy)
}

/// Check if SIMD is available at runtime
///
/// Returns true if the `simd` feature is enabled and we're running on WASM.
/// This can be used to log SIMD status or conditionally enable optimizations.
#[cfg(target_arch = "wasm32")]
pub fn is_simd_available() -> bool {
    cfg!(feature = "simd")
}

/// Check if SIMD is available at runtime (non-WASM fallback)
///
/// Always returns false on non-WASM platforms since WASM SIMD128 is
/// only available in WebAssembly.
#[cfg(not(target_arch = "wasm32"))]
pub fn is_simd_available() -> bool {
    false
}

/// Apply velocity updates using SIMD
///
/// Updates velocities based on forces and applies damping.
///
/// # Arguments
/// * `vx` - X velocities (modified in place)
/// * `vy` - Y velocities (modified in place)
/// * `forces` - Force accumulator
/// * `alpha` - Simulation alpha (energy)
/// * `damping` - Velocity decay factor
#[cfg(all(target_arch = "wasm32", feature = "simd"))]
#[target_feature(enable = "simd128")]
pub unsafe fn apply_velocity_update_simd(
    vx: &mut [f32],
    vy: &mut [f32],
    forces: &SimdForces,
    alpha: f32,
    damping: f32,
) {
    let n = vx.len().min(vy.len()).min(forces.len);
    let alpha_v = f32x4_splat(alpha);
    let damping_v = f32x4_splat(damping);

    let chunks = n / 4;
    for chunk in 0..chunks {
        let base = chunk * 4;

        // Load velocities
        let vx_v = f32x4(vx[base], vx[base + 1], vx[base + 2], vx[base + 3]);
        let vy_v = f32x4(vy[base], vy[base + 1], vy[base + 2], vy[base + 3]);

        // Load forces
        let fx_v = f32x4(
            forces.fx[base],
            forces.fx[base + 1],
            forces.fx[base + 2],
            forces.fx[base + 3],
        );
        let fy_v = f32x4(
            forces.fy[base],
            forces.fy[base + 1],
            forces.fy[base + 2],
            forces.fy[base + 3],
        );

        // Update: v = (v + f * alpha) * damping
        let new_vx = f32x4_mul(f32x4_add(vx_v, f32x4_mul(fx_v, alpha_v)), damping_v);
        let new_vy = f32x4_mul(f32x4_add(vy_v, f32x4_mul(fy_v, alpha_v)), damping_v);

        // Store
        vx[base] = f32x4_extract_lane::<0>(new_vx);
        vx[base + 1] = f32x4_extract_lane::<1>(new_vx);
        vx[base + 2] = f32x4_extract_lane::<2>(new_vx);
        vx[base + 3] = f32x4_extract_lane::<3>(new_vx);

        vy[base] = f32x4_extract_lane::<0>(new_vy);
        vy[base + 1] = f32x4_extract_lane::<1>(new_vy);
        vy[base + 2] = f32x4_extract_lane::<2>(new_vy);
        vy[base + 3] = f32x4_extract_lane::<3>(new_vy);
    }

    // Handle remainder
    for i in (chunks * 4)..n {
        vx[i] = (vx[i] + forces.fx[i] * alpha) * damping;
        vy[i] = (vy[i] + forces.fy[i] * alpha) * damping;
    }
}

/// Non-SIMD fallback for velocity update
#[cfg(not(all(target_arch = "wasm32", feature = "simd")))]
pub fn apply_velocity_update_simd(
    vx: &mut [f32],
    vy: &mut [f32],
    forces: &SimdForces,
    alpha: f32,
    damping: f32,
) {
    apply_velocity_update_scalar(vx, vy, forces, alpha, damping)
}

/// Scalar implementation of velocity update (fallback)
pub fn apply_velocity_update_scalar(
    vx: &mut [f32],
    vy: &mut [f32],
    forces: &SimdForces,
    alpha: f32,
    damping: f32,
) {
    let n = vx.len().min(vy.len()).min(forces.len);
    for i in 0..n {
        vx[i] = (vx[i] + forces.fx[i] * alpha) * damping;
        vy[i] = (vy[i] + forces.fy[i] * alpha) * damping;
    }
}

/// Apply position updates using SIMD
///
/// Updates positions based on velocities.
#[cfg(all(target_arch = "wasm32", feature = "simd"))]
#[target_feature(enable = "simd128")]
pub unsafe fn apply_position_update_simd(x: &mut [f32], y: &mut [f32], vx: &[f32], vy: &[f32]) {
    let n = x.len().min(y.len()).min(vx.len()).min(vy.len());

    let chunks = n / 4;
    for chunk in 0..chunks {
        let base = chunk * 4;

        // Load positions
        let x_v = f32x4(x[base], x[base + 1], x[base + 2], x[base + 3]);
        let y_v = f32x4(y[base], y[base + 1], y[base + 2], y[base + 3]);

        // Load velocities
        let vx_v = f32x4(vx[base], vx[base + 1], vx[base + 2], vx[base + 3]);
        let vy_v = f32x4(vy[base], vy[base + 1], vy[base + 2], vy[base + 3]);

        // Update: pos += velocity
        let new_x = f32x4_add(x_v, vx_v);
        let new_y = f32x4_add(y_v, vy_v);

        // Store
        x[base] = f32x4_extract_lane::<0>(new_x);
        x[base + 1] = f32x4_extract_lane::<1>(new_x);
        x[base + 2] = f32x4_extract_lane::<2>(new_x);
        x[base + 3] = f32x4_extract_lane::<3>(new_x);

        y[base] = f32x4_extract_lane::<0>(new_y);
        y[base + 1] = f32x4_extract_lane::<1>(new_y);
        y[base + 2] = f32x4_extract_lane::<2>(new_y);
        y[base + 3] = f32x4_extract_lane::<3>(new_y);
    }

    // Handle remainder
    for i in (chunks * 4)..n {
        x[i] += vx[i];
        y[i] += vy[i];
    }
}

/// Non-SIMD fallback for position update
#[cfg(not(all(target_arch = "wasm32", feature = "simd")))]
pub fn apply_position_update_simd(x: &mut [f32], y: &mut [f32], vx: &[f32], vy: &[f32]) {
    apply_position_update_scalar(x, y, vx, vy)
}

/// Scalar implementation of position update (fallback)
pub fn apply_position_update_scalar(x: &mut [f32], y: &mut [f32], vx: &[f32], vy: &[f32]) {
    let n = x.len().min(y.len()).min(vx.len()).min(vy.len());
    for i in 0..n {
        x[i] += vx[i];
        y[i] += vy[i];
    }
}

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

    #[test]
    fn test_simd_positions_creation() {
        let positions = vec![(0.0, 0.0), (1.0, 1.0), (2.0, 2.0)];
        let simd_pos = SimdPositions::from_positions(&positions);

        assert_eq!(simd_pos.len, 3);
        assert_eq!(simd_pos.x[0], 0.0);
        assert_eq!(simd_pos.y[1], 1.0);
        assert_eq!(simd_pos.mass[2], 1.0);
    }

    #[test]
    fn test_simd_forces_zeros() {
        let forces = SimdForces::zeros(5);
        assert_eq!(forces.len, 5);
        assert!(forces.fx.iter().all(|&f| f == 0.0));
        assert!(forces.fy.iter().all(|&f| f == 0.0));
    }

    #[test]
    fn test_scalar_repulsion() {
        let positions = vec![(0.0, 0.0), (10.0, 0.0)];
        let simd_pos = SimdPositions::from_positions(&positions);

        let forces = calculate_repulsion_scalar(&simd_pos, -100.0);

        // With negative strength (repulsion), delta is (0-10, 0) = (-10, 0)
        // The force magnitude is negative, so the force points from 0 toward 1
        // Node 0 at (0,0) is pushed away from node 1 at (10,0)
        // Force direction: delta.normalize() * (strength / dist_sq) = (-1,0) * (-100/100) = (1, 0)
        // So forces.fx[0] should be positive (pushed right, away from repelling force)
        assert!(forces.fx[0] != 0.0, "Force on node 0 should be non-zero");
        // Force on node 1 should be opposite of force on node 0
        assert!(
            (forces.fx[0] + forces.fx[1]).abs() < 0.001,
            "Forces should be equal and opposite"
        );
        // Y forces should be zero
        assert!(forces.fy[0].abs() < 0.001);
        assert!(forces.fy[1].abs() < 0.001);
    }

    #[test]
    fn test_scalar_center_force() {
        let positions = vec![(10.0, 10.0), (-10.0, -10.0)];
        let simd_pos = SimdPositions::from_positions(&positions);

        let forces = calculate_center_force_scalar(&simd_pos, 0.0, 0.0, 0.1);

        // Node at (10, 10) should be pulled toward origin
        assert!(forces.fx[0] < 0.0);
        assert!(forces.fy[0] < 0.0);
        // Node at (-10, -10) should be pulled toward origin
        assert!(forces.fx[1] > 0.0);
        assert!(forces.fy[1] > 0.0);
    }

    #[test]
    fn test_scalar_leaf_force() {
        let sources = SimdPositions::from_positions(&[(10.0, 0.0), (0.0, 10.0)]);
        let (fx, fy) = calculate_leaf_force_scalar(0.0, 0.0, &sources, -100.0);

        // With target at (0,0) and sources at (10,0) and (0,10):
        // Delta to source 1: (0-10, 0-0) = (-10, 0), normalized = (-1, 0)
        // Delta to source 2: (0-0, 0-10) = (0, -10), normalized = (0, -1)
        // With strength=-100 and dist_sq=100: force_mag = -100/100 = -1
        // Force from source 1: (-1, 0) * -1 = (1, 0) [repelled rightward]
        // Force from source 2: (0, -1) * -1 = (0, 1) [repelled upward]
        // Total: (1, 1) - both positive (pushed away from sources)
        // But wait, the calculation is target - source, so:
        // delta1 = (0-10, 0) = (-10, 0), force points left (negative x)
        // Actually let me trace through: delta = target - source
        // For source at (10,0): delta = (0-10, 0-0) = (-10, 0)
        // dist = 10, force_mag = -100 * 1 / 100 = -1
        // force = (-10/10) * -1 = 1 (positive x, pushed away from source on right)
        // So fx should be positive, fy should be positive
        assert!(fx != 0.0, "fx should be non-zero");
        assert!(fy != 0.0, "fy should be non-zero");
        // Forces should have same sign since sources are in positive quadrant
        // and repulsion pushes toward negative quadrant
        // Actually, force = delta.normalize() * force_mag
        // delta for (10,0) = (-10,0), normalized = (-1,0), force_mag = -1, so (-1)*(-1) = 1
        // delta for (0,10) = (0,-10), normalized = (0,-1), force_mag = -1, so (-1)*(-1) = 1
        // So both forces are positive (pushed away from sources toward negative)
        // Wait no - if sources are at positive coords, repulsion pushes target toward negative
        // Let me just check they're non-zero and consistent
        assert!(
            (fx.abs() - fy.abs()).abs() < 0.001,
            "Forces should be symmetric"
        );
    }

    #[test]
    fn test_velocity_update_scalar() {
        let mut vx = vec![1.0, 2.0, 3.0];
        let mut vy = vec![1.0, 2.0, 3.0];
        let forces = SimdForces {
            fx: vec![1.0, 1.0, 1.0],
            fy: vec![1.0, 1.0, 1.0],
            len: 3,
        };

        apply_velocity_update_scalar(&mut vx, &mut vy, &forces, 1.0, 0.5);

        // v = (v + f * alpha) * damping = (1 + 1) * 0.5 = 1.0
        assert!((vx[0] - 1.0).abs() < 0.001);
    }

    #[test]
    fn test_position_update_scalar() {
        let mut x = vec![0.0, 1.0, 2.0];
        let mut y = vec![0.0, 1.0, 2.0];
        let vx = vec![1.0, 1.0, 1.0];
        let vy = vec![1.0, 1.0, 1.0];

        apply_position_update_scalar(&mut x, &mut y, &vx, &vy);

        assert_eq!(x[0], 1.0);
        assert_eq!(y[1], 2.0);
    }

    #[test]
    fn test_forces_to_vectors() {
        let forces = SimdForces {
            fx: vec![1.0, 2.0],
            fy: vec![3.0, 4.0],
            len: 2,
        };

        let vectors = forces.to_vectors();
        assert_eq!(vectors.len(), 2);
        assert!((vectors[0].x - 1.0).abs() < 0.001);
        assert!((vectors[0].y - 3.0).abs() < 0.001);
    }

    #[test]
    fn test_forces_add() {
        let mut forces1 = SimdForces {
            fx: vec![1.0, 2.0],
            fy: vec![3.0, 4.0],
            len: 2,
        };
        let forces2 = SimdForces {
            fx: vec![1.0, 1.0],
            fy: vec![1.0, 1.0],
            len: 2,
        };

        forces1.add(&forces2);

        assert_eq!(forces1.fx[0], 2.0);
        assert_eq!(forces1.fy[1], 5.0);
    }
}