// Package flock is the boid-flock workload's pure Odin subsystem: a
// deterministic, stateful 3D boid simulation with a uniform spatial hash,
// plus the staged MultiMesh transform writer the authored class submits.
// It imports neither Godot, Gododin, nor Raylib, and after `flock_create`
// returns it allocates nothing: every array lives in one arena-backed block
// carved at creation and reused by every frame forever.
//
// Determinism is a property of the whole pipeline: per-agent state is seeded
// from (seed, index) by integer bit mixing exactly in the spirit of the
// multimesh fixture's `stream_parameter` — multiply-xor rounds with rational
// constants so every seed and index bit reaches every state byte — and
// `flock_step` performs only f32 arithmetic in a fixed order over fixed
// arrays with a caller-provided dt and no clock reads. The same seed,
// parameters, and step count therefore replay to the identical state on a
// host, which `state_digest` pins. There is no from-first-principles
// instance oracle here the way the stateless multimesh producer has one:
// every boid's next state depends on its neighbors, so replay requires the
// full configuration and history. Small all-pairs reference tests check physics. The verification surface is
// honest about that:
//
//   - `state_digest` pins the simulation across runs and transports;
//   - `clone` lets a test replay a second flock over the same history;
//   - `verify_instance_against` recomputes one boid's twelve written floats
//     from the flock's own current state, which checks the *layout writer*
//     against the state it claims to encode, not the physics;
//   - `verify_readback` permutes the engine's Transform3D getter shape back
//     to the written buffer shape and compares bitwise, which checks the
//     engine read path against exactly the floats this package wrote.
//
// The buffer layout is the locked engine contract, not a local assumption,
// carried over unchanged from the multimesh fixture: RenderingServer
// multimesh data is one f32 per component with `stride_cache = 12` floats
// per instance for MULTIMESH_TRANSFORM_3D with no colors and no custom data
// (`servers/rendering/renderer_rd/storage_rd/mesh_storage.cpp`,
// `_multimesh_allocate_data` at line 1577: color offset 12, custom-data
// offset 12, stride 12, re-verified in the pinned 4.8-dev5 source). The
// twelve components are interleaved, not a plain Transform3D memory
// copy (`_multimesh_instance_get_transform` at line 2036): basis row 0,
// then origin.x, basis row 1, then origin.y, basis row 2, then origin.z.
// `_multimesh_set_buffer` validates the submitted size against
// `instances * stride` and hands the caller's floats straight to
// `RD::buffer_update`, so the engine consumes the submission synchronously
// and retains nothing.
package flock

import "core:math"
import "core:mem"

// FLOATS_PER_INSTANCE is the locked TRANSFORM_3D stride without colors or
// custom data: three basis rows plus one origin row, interleaved.
FLOATS_PER_INSTANCE :: 12

// MIN_PERCEPTION is the smallest perception radius any step accepts. The
// spatial hash's backing storage is sized once at creation for this floor
// (see `min_cell_size`), so a live parameter change can never outgrow the
// arena.
MIN_PERCEPTION :: 0.5

// MAX_AGENTS is the largest flock `flock_create` accepts, mirroring the
// multimesh fixture's capacity ceiling.
MAX_AGENTS :: 1_000_000

// Supported numerical domain. These generous bounds exceed the demo controls
// while keeping squared distances, weighted steering and worst-case million-
// agent sums finite in f32. This kernel deliberately refuses unbounded physics.
MAX_HALF_EXTENT :: 1_000_000.0
MAX_PERCEPTION :: 1_000_000.0
MAX_SPEED :: 10_000.0
MAX_WEIGHT :: 10_000.0
MAX_PREDATOR_COORDINATE :: 10_000_000.0
MAX_DT :: 1.0

// GRID_MAX_DIM bounds the spatial hash's cells per axis so a small
// perception radius over a large volume cannot explode the cell arrays; the
// effective cell size grows instead. Cells only ever get *larger* than the
// perception radius, never smaller, so the 27-cell neighborhood always
// covers every boid within perception.
GRID_MAX_DIM :: 64

// Internal steering constants. They are fixed, rational, and f32-exact so
// they take no part in any claim: the reviewed knobs are the Params fields.
// The rules are Reynolds steering: each rule builds a desired velocity of
// magnitude max_speed, forms the steering (desired - velocity), clamps it to
// a per-rule force bound, and scales it by its weight. Proportional
// attraction (a raw (center - position) * weight force) is deliberately
// avoided: it has no speed bound at range and condenses the whole flock into
// a point, which the workload's measurements exposed and this form removes.
MAX_FORCE_FACTOR       :: 0.5 // per-rule force bound = max_speed * this
STEERING_TOTAL_FACTOR  :: 4.0 // total acceleration cap = max_speed * this
BOUNDS_GAIN            :: 24.0 // wall pushback stiffness (1/s^2)
BOUNDS_MARGIN_FACTOR   :: 0.5 // wall margin = perception * this
PREDATOR_RADIUS_FACTOR :: 2.0 // flee radius = perception * this
PREDATOR_STRENGTH      :: 2.0 // flee gain = separation weight * this
FACING_EPSILON2        :: 1e-8 // min speed^2 that refreshes the facing vector

// Params holds one flock's tunable behavior. Weights, perception, and
// max_speed are live-updatable between steps; the world half-extents are a
// creation fact fixed into the arena sizing. Both `half_extents` and
// `predator` are plain [3]f32 records in x, y, z order.
Params :: struct {
	separation:      f32,
	alignment:       f32,
	cohesion:        f32,
	perception:      f32,
	max_speed:       f32,
	half_extents:    [3]f32,
	predator:        [3]f32,
	predator_active: bool,
}

// DEFAULT_PARAMS is the example's starting tuning: a 60x30x60 world whose
// density at 10,000 agents (~0.0926 agents/unit^3) matches the 3.5-unit
// perception radius for a healthy neighbor count, and balanced weights.
DEFAULT_PARAMS :: Params{
	separation   = 1.5,
	alignment    = 1.0,
	cohesion     = 0.9,
	perception   = 3.5,
	max_speed    = 6.0,
	half_extents = {30.0, 15.0, 30.0},
}

// Flock is the sole allocation owner for the simulation. Its struct header
// and every SoA and spatial-hash array live in one arena-backed allocation
// created by `flock_create` and released by `flock_destroy`; no step, write,
// or query allocates.
Flock :: struct {
	live:          bool,
	allocator:     mem.Allocator,
	count:         i32,
	frame:         i64,
	max_neighbors: i32,
	params:        Params, // extents and arena sizing fixed at creation

	// Structure-of-arrays agent state, views into the arena block.
	px, py, pz: []f32, // positions
	vx, vy, vz: []f32, // velocities
	fx, fy, fz: []f32, // last-known unit facing, refreshed only when speed > epsilon
	ax, ay, az: []f32, // per-step acceleration scratch (two-pass integration)

	// Uniform spatial hash. cell_size >= perception always holds, so the 27
	// neighboring cells cover every boid within perception. Rebuilt every
	// step by a two-pass counting sort into these preallocated arrays.
	// Test-and-guard hook: when validate_hash is set, rebuild_hash checks
	// the freshly built hash (structure plus per-item placement against the
	// positions it was built from) and records the verdict. Reading it after
	// flock_step returns is sound; reading the raw arrays' placement after
	// integration is not, because the agents have moved since the build.
	validate_hash:     bool,
	hash_valid_observed: bool,

	min_cell_size: f32,
	cell_size:     f32,
	grid_x:        i32,
	grid_y:        i32,
	grid_z:        i32,
	cells:         i32,
	cell_counts:   []i32, // per-cell counting cursor; equals per-cell totals after the build
	cell_starts:   []i32, // exclusive prefix sums, length cells+1
	cell_items:    []i32, // agent indices sorted by cell, ascending within a cell
}

// params_valid reports whether `params` is admissible for creation: finite
// positive extents, a perception at or above the arena floor, a positive
// max_speed, and non-negative finite weights, all inside the supported domain.
params_valid :: proc(params: Params) -> bool {
	if !(params.half_extents.x > 0) || !(params.half_extents.y > 0) || !(params.half_extents.z > 0) {
		return false
	}
	if !finite(params.half_extents.x) || !finite(params.half_extents.y) ||
	   !finite(params.half_extents.z) || params.half_extents.x > MAX_HALF_EXTENT ||
	   params.half_extents.y > MAX_HALF_EXTENT || params.half_extents.z > MAX_HALF_EXTENT {
		return false
	}
	if !finite(params.predator[0]) || !finite(params.predator[1]) ||
	   !finite(params.predator[2]) || math.abs(params.predator[0]) > MAX_PREDATOR_COORDINATE ||
	   math.abs(params.predator[1]) > MAX_PREDATOR_COORDINATE || math.abs(params.predator[2]) > MAX_PREDATOR_COORDINATE {
		return false
	}
	if !finite(params.perception) || params.perception < MIN_PERCEPTION || params.perception > MAX_PERCEPTION {
		return false
	}
	if !finite(params.max_speed) || params.max_speed <= 0 || params.max_speed > MAX_SPEED {
		return false
	}
	if !finite(params.separation) || params.separation < 0 || params.separation > MAX_WEIGHT {
		return false
	}
	if !finite(params.alignment) || params.alignment < 0 || params.alignment > MAX_WEIGHT {
		return false
	}
	if !finite(params.cohesion) || params.cohesion < 0 || params.cohesion > MAX_WEIGHT {
		return false
	}
	return true
}

// steering_valid checks every live field before a step changes state. Extents
// remain creation facts in f.params; all finite weights and predator coordinates
// still matter when an inactive feature is enabled on a later frame.
steering_valid :: proc(params: Params) -> bool {
	return params_valid(params)
}

// min_cell_size is the smallest cell size the arena was sized for: the
// perception floor, unless the volume's longest axis would need more than
// GRID_MAX_DIM cells, in which case the floor rises so the grid fits.
min_cell_size :: proc(params: Params) -> f32 {
	longest := max(params.half_extents.x, max(params.half_extents.y, params.half_extents.z))
	return max(MIN_PERCEPTION, 2.0 * longest / f32(GRID_MAX_DIM))
}

// grid_dims computes the cell count per axis for a cell size, clamped into
// [1, GRID_MAX_DIM].
grid_dims :: proc(params: Params, cell: f32) -> (i32, i32, i32) {
	dx := clamp_i32(i32(math.ceil(2.0 * params.half_extents.x / cell)), 1, GRID_MAX_DIM)
	dy := clamp_i32(i32(math.ceil(2.0 * params.half_extents.y / cell)), 1, GRID_MAX_DIM)
	dz := clamp_i32(i32(math.ceil(2.0 * params.half_extents.z / cell)), 1, GRID_MAX_DIM)
	return dx, dy, dz
}

// finite reports whether the f32 value is neither NaN nor infinite (this
// Odin core names is_nan/is_inf but no is_finite).
finite :: proc "contextless" (value: f32) -> bool {
	return !math.is_nan(value) && !math.is_inf(value)
}

clamp_i32 :: proc(value, low, high: i32) -> i32 {
	return min(high, max(low, value))
}

clamp_f32 :: proc(value, low, high: f32) -> f32 {
	return min(high, max(low, value))
}

mix_round :: proc(value: u32) -> u32 {
	return value * 2654435761 ~ (value >> 16)
}

// agent_seed derives one agent's bit stream from the population seed, the
// agent index, and a per-stream odd multiplier: a multiply-xor cascade so
// every seed and index bit reaches every parameter byte, and a one-bit seed
// change moves every agent (the stream_parameter discipline).
agent_seed :: proc(index: i32, seed: i64, stream: u32) -> u32 {
	mixed := u32(seed) ~ (u32(seed >> 32) * 0x9E3779B9)
	mixed = mixed ~ (u32(index) * stream)
	mixed = mix_round(mixed)
	mixed = mix_round(mixed)
	return mixed
}

// unit01 maps 16 bits onto [0, 1) exactly: the value is k/65536 for an
// integer k, so every arithmetic that follows starts from an exact f32.
unit01 :: proc(bits: u32) -> f32 {
	return f32((bits >> 8) & 0xFFFF) * (1.0 / 65536.0)
}

// Create_Status names the last flock_create failure reason for diagnostics:
// 0 means the last call succeeded, 1..4 are the refused branches. Main-thread
// diagnostics only; no step reads it.
Create_Status :: enum u8 {
	Ok,
	Bad_Count,
	Bad_Params,
	Bad_Grid,
	Allocation,
}
create_status: Create_Status = .Ok

// flock_create allocates the whole simulation in one arena-backed
// allocation: the struct header plus every SoA and spatial-hash array, the
// hash sized for the worst grid the live perception floor can produce. The
// per-agent seeds mix the population seed and the agent index so a one-bit
// seed change moves every agent, and initial positions land inside 80% of
// the world bounds moving at half to full max speed.
flock_create :: proc(count: i32, seed: i64, params: Params, allocator: mem.Allocator = context.allocator) -> (^Flock, bool) {
	create_status = .Ok
	if count <= 0 || count > MAX_AGENTS {
		create_status = .Bad_Count
		return nil, false
	}
	if !params_valid(params) {
		create_status = .Bad_Params
		return nil, false
	}

	cell := min_cell_size(params)
	gx, gy, gz := grid_dims(params, cell)
	cells := gx * gy * gz
	if cells <= 0 {
		create_status = .Bad_Grid
		return nil, false
	}

	n := int(count)
	float_values := 12 * n                            // px..vz (6), facing (3), acceleration (3)
	i32_values := 2*int(cells) + 1 + n                // counts + starts + items
	header_bytes := (size_of(Flock) + 7) / 8 * 8      // header, 8-aligned
	float_bytes := (4*float_values + 7) / 8 * 8       // f32 arrays, 8-aligned
	storage, allocation_error := mem.alloc(header_bytes + float_bytes + 4*i32_values, 8, allocator)
	if allocation_error != nil || storage == nil {
		create_status = .Allocation
		return nil, false
	}

	// Carve the typed views out of the single block: f32 arrays first, then
	// the i32 hash arrays. mem.alloc zeroed the storage; only the non-zero
	// owner facts are stated below.
	block := cast([^]byte)storage
	fbase := cast([^]f32)rawptr(&block[header_bytes])
	ibase := cast([^]i32)rawptr(&block[header_bytes + float_bytes])

	f := cast(^Flock)storage
	f.px = fbase[0*n:n]
	f.py = fbase[1*n:2*n]
	f.pz = fbase[2*n:3*n]
	f.vx = fbase[3*n:4*n]
	f.vy = fbase[4*n:5*n]
	f.vz = fbase[5*n:6*n]
	f.fx = fbase[6*n:7*n]
	f.fy = fbase[7*n:8*n]
	f.fz = fbase[8*n:9*n]
	f.ax = fbase[9*n:10*n]
	f.ay = fbase[10*n:11*n]
	f.az = fbase[11*n:12*n]
	c := int(cells)
	f.cell_counts = ibase[0:c]
	f.cell_starts = ibase[c : 2*c+1]
	f.cell_items = ibase[2*c+1 : 2*c+1+n]

	f.live = true
	f.allocator = allocator
	f.count = count
	f.frame = 0
	f.max_neighbors = 0
	f.params = params
	f.min_cell_size = cell
	f.cell_size = cell
	f.grid_x = gx
	f.grid_y = gy
	f.grid_z = gz
	f.cells = cells

	for index in 0..<count {
		seed_x := agent_seed(index, seed, 2654435761)
		seed_v := agent_seed(index, seed, 2246822519)
		seed_s := agent_seed(index, seed, 3266489917)

		// Position: uniform inside 80% of the world box from exact f32
		// inputs.
		f.px[index] = (2.0 * unit01(seed_x) - 1.0) * params.half_extents.x * 0.8
		f.py[index] = (2.0 * unit01(seed_x * 65539) - 1.0) * params.half_extents.y * 0.8
		f.pz[index] = (2.0 * unit01(seed_x * 40503) - 1.0) * params.half_extents.z * 0.8

		// Velocity: a uniform direction on the sphere (the vertical from one
		// uniform, the ring angle from another) at half to full max speed.
		z := 2.0 * unit01(seed_v) - 1.0
		ring := math.sqrt(max(0.0, 1.0 - z * z))
		theta := 6.2831853 * unit01(seed_v * 65539)
		speed := params.max_speed * (0.5 + 0.5 * unit01(seed_s))
		dir_x := ring * math.cos(theta)
		dir_y := z
		dir_z := ring * math.sin(theta)
		f.vx[index] = dir_x * speed
		f.vy[index] = dir_y * speed
		f.vz[index] = dir_z * speed

		// Facing starts at the initial heading; it refreshes only when the
		// agent actually moves, so a stopped agent keeps its orientation.
		f.fx[index] = dir_x
		f.fy[index] = dir_y
		f.fz[index] = dir_z
	}
	return f, true
}

// flock_destroy releases the arena and clears the caller's pointer. It is
// safe against double destruction and tolerates any partial state.
flock_destroy :: proc(f: ^^Flock) -> bool {
	if f == nil || f^ == nil {
		return false
	}
	owner := f^
	if !owner.live {
		return false
	}
	allocator := owner.allocator
	owner.live = false
	if mem.free(rawptr(owner), allocator) != nil {
		owner.live = true
		return false
	}
	f^ = nil
	return true
}

// cell_of maps a world position onto its grid cell, clamped into the grid.
// Positions stay inside the bounds by construction, so the clamp is a
// robustness net, not a behavior.
cell_of :: proc(x, y, z: f32, params: Params, cell: f32, gx, gy, gz: i32) -> i32 {
	cx := clamp_i32(i32(math.floor((x + params.half_extents.x) / cell)), 0, gx - 1)
	cy := clamp_i32(i32(math.floor((y + params.half_extents.y) / cell)), 0, gy - 1)
	cz := clamp_i32(i32(math.floor((z + params.half_extents.z) / cell)), 0, gz - 1)
	return (cz * gy + cy) * gx + cx
}

// steer_rule accumulates one Reynolds steering: `direction` (any nonzero
// vector) becomes a desired velocity of magnitude `speed`, the steering
// (desired - current velocity) is clamped to `force`, scaled by `weight`,
// and added to the caller's accumulator. Returns false and adds nothing when
// the direction degenerates. Deterministic f32 throughout.
steer_rule :: proc(acc_x, acc_y, acc_z: ^f32, dx, dy, dz, vx, vy, vz, speed, force, weight: f32) -> bool {
	length2 := dx*dx + dy*dy + dz*dz
	if !(length2 > 1e-12) {
		return false
	}
	scale := speed / math.sqrt(length2)
	sx := dx * scale - vx
	sy := dy * scale - vy
	sz := dz * scale - vz
	s2 := sx*sx + sy*sy + sz*sz
	if s2 > force*force {
		shrink := force / math.sqrt(s2)
		sx *= shrink
		sy *= shrink
		sz *= shrink
	}
	acc_x^ += sx * weight
	acc_y^ += sy * weight
	acc_z^ += sz * weight
	return true
}

// rebuild_hash refills the spatial hash for the current positions with a
// two-pass counting sort: count agents per cell, exclusive-prefix the
// starts, then scatter ascending agent indices through a cursor reused from
// cell_counts. Within every cell the items therefore ascend by agent index,
// which makes every neighbor query a deterministic sequence no matter where
// the agents wandered.
rebuild_hash :: proc(f: ^Flock, perception: f32) {
	cell := max(perception, f.min_cell_size)
	if cell != f.cell_size {
		f.cell_size = cell
		f.grid_x, f.grid_y, f.grid_z = grid_dims(f.params, cell)
		f.cells = f.grid_x * f.grid_y * f.grid_z
	}
	counts := f.cell_counts
	starts := f.cell_starts
	items := f.cell_items
	for c in 0..<f.cells {
		counts[c] = 0
	}
	starts[0] = 0
	for index in 0..<f.count {
		c := cell_of(f.px[index], f.py[index], f.pz[index], f.params, cell, f.grid_x, f.grid_y, f.grid_z)
		counts[c] += 1
	}
	for c in 0..<f.cells {
		starts[c+1] = starts[c] + counts[c]
	}
	for c in 0..<f.cells {
		counts[c] = 0 // reuse as the per-cell scatter cursor
	}
	for index in 0..<f.count {
		c := cell_of(f.px[index], f.py[index], f.pz[index], f.params, cell, f.grid_x, f.grid_y, f.grid_z)
		items[starts[c] + counts[c]] = index
		counts[c] += 1
	}
	if f.validate_hash {
		f.hash_valid_observed = flock_hash_valid(f)
	}
}

// flock_hash_structure_valid checks the hash invariants that hold at any
// time, integration or not: slice lengths, an exclusive prefix sum of starts
// ending exactly at count, and per-cell counts consistent with their
// windows. It allocates nothing.
flock_hash_structure_valid :: proc(f: ^Flock) -> bool {
	if f == nil || !f.live {
		return false
	}
	counts := f.cell_counts
	starts := f.cell_starts
	if len(counts) < int(f.cells) || len(starts) < int(f.cells)+1 || len(f.cell_items) < int(f.count) {
		return false
	}
	if starts[0] != 0 {
		return false
	}
	for c in 0..<f.cells {
		if starts[c+1] < starts[c] {
			return false
		}
		if counts[c] != starts[c+1] - starts[c] {
			return false
		}
	}
	return starts[f.cells] == f.count
}

// flock_hash_valid adds the placement invariant to the structural one:
// every item sits inside its cell's window, in range, and actually belongs
// to the cell it appears under. The placement half is only meaningful
// against the positions the hash was built from, so tests observe it through
// `hash_valid_observed` (armed by `validate_hash` inside rebuild_hash)
// rather than calling this after a step has already integrated the flock.
flock_hash_valid :: proc(f: ^Flock) -> bool {
	if !flock_hash_structure_valid(f) {
		return false
	}
	starts := f.cell_starts
	items := f.cell_items
	for c in 0..<f.cells {
		for k in starts[c]..<starts[c+1] {
			j := items[k]
			if j < 0 || j >= f.count {
				return false
			}
			actual := cell_of(f.px[j], f.py[j], f.pz[j], f.params, f.cell_size, f.grid_x, f.grid_y, f.grid_z)
			if actual != c {
				return false
			}
		}
	}
	return true
}

// flock_step advances the simulation one caller-provided dt. The pass
// structure is fixed: rebuild the spatial hash, gather every agent's
// neighbors from the 27 surrounding cells in deterministic order into
// acceleration scratch, then integrate the whole flock with semi-implicit
// Euler and soft-clamp positions into the bounds. Nothing allocates and no
// clock is read; dt arrives from the caller, and the same seed, parameters,
// and step count replay to the identical state.
flock_step :: proc(f: ^Flock, dt: f32, params: Params) -> bool {
	if f == nil || !f.live {
		return false
	}
	if !finite(dt) || dt <= 0 || dt > MAX_DT {
		return false
	}
	if !steering_valid(params) {
		return false
	}

	rebuild_hash(f, params.perception)

	perception2 := params.perception * params.perception
	margin := params.perception * BOUNDS_MARGIN_FACTOR
	force := params.max_speed * MAX_FORCE_FACTOR
	total_limit2 := params.max_speed * STEERING_TOTAL_FACTOR
	total_limit2 = total_limit2 * total_limit2
	flee_radius := params.perception * PREDATOR_RADIUS_FACTOR
	flee_radius2 := flee_radius * flee_radius
	flee_gain := params.separation * PREDATOR_STRENGTH

	max_neighbors := i32(0)

	// Pass one: steering. Neighbors come from the 27 cells around the
	// agent's own cell, iterated in fixed z/y/x order (x varies fastest), each cell's
	// items in ascending agent index; every candidate is filtered by the
	// actual squared distance against the perception radius.
	for i in 0..<f.count {
		px := f.px[i]
		py := f.py[i]
		pz := f.pz[i]
		sepx, sepy, sepz := f32(0), f32(0), f32(0)
		alix, aliy, aliz := f32(0), f32(0), f32(0)
		cohx, cohy, cohz := f32(0), f32(0), f32(0)
		neighbors := i32(0)

		c := cell_of(px, py, pz, f.params, f.cell_size, f.grid_x, f.grid_y, f.grid_z)
		gx, gy, gz := int(f.grid_x), int(f.grid_y), int(f.grid_z)
		cz := int(c) / (gx * gy)
		rem := int(c) - cz * gx * gy
		cy := rem / gx
		cx := rem - cy * gx

		first_x := max(0, cx - 1)
		after_x := min(gx, cx + 2)

		for dz in -1..=1 {
			nz := cz + dz
			if nz < 0 || nz >= gz {
				continue
			}
			for dy in -1..=1 {
				ny := cy + dy
				if ny < 0 || ny >= gy {
					continue
				}
				// Adjacent x cells occupy one contiguous counting-sort span.
				// Joining their ranges preserves the old candidate and f32 sum
				// order, including empty cells and clamped edge neighborhoods.
				row := (nz * gy + ny) * gx
				for k in f.cell_starts[row+first_x]..<f.cell_starts[row+after_x] {
					j := f.cell_items[k]
					if j == i {
						continue
					}
					dxp := px - f.px[j]
					dyp := py - f.py[j]
					dzp := pz - f.pz[j]
					d2 := dxp*dxp + dyp*dyp + dzp*dzp
					if d2 >= perception2 {
						continue
					}
					neighbors += 1
					d := math.sqrt(d2)
					if d > 1e-6 {
						weight := 1.0 / d // separation falls off as 1/distance
						sepx += dxp * weight
						sepy += dyp * weight
						sepz += dzp * weight
					}
					alix += f.vx[j]
					aliy += f.vy[j]
					aliz += f.vz[j]
					cohx += f.px[j]
					cohy += f.py[j]
					cohz += f.pz[j]
				}
			}
		}
		if neighbors > max_neighbors {
			max_neighbors = neighbors
		}

		// Reynolds steering per rule: a unit desired direction scaled to
		// max_speed, steered toward by a force clamped to `force` and scaled
		// by the rule's weight. Each rule therefore pulls with at most a
		// bounded, speed-shaped force no matter how far the target is.
		acx, acy, acz := f32(0), f32(0), f32(0)
		if neighbors > 0 {
			steer_rule(&acx, &acy, &acz, sepx, sepy, sepz, f.vx[i], f.vy[i], f.vz[i], params.max_speed, force, params.separation)
			n := f32(neighbors)
			steer_rule(&acx, &acy, &acz, alix / n, aliy / n, aliz / n, f.vx[i], f.vy[i], f.vz[i], params.max_speed, force, params.alignment)
			steer_rule(&acx, &acy, &acz, cohx / n - px, cohy / n - py, cohz / n - pz, f.vx[i], f.vy[i], f.vz[i], params.max_speed, force, params.cohesion)
		}

		// Bounds steering: a smooth pushback that ramps up over the margin
		// band inside each wall, never a hard clamp at this stage.
		ex := f.params.half_extents.x - margin
		ey := f.params.half_extents.y - margin
		ez := f.params.half_extents.z - margin
		if px > ex {
			acx -= BOUNDS_GAIN * (px - ex)
		}
		if px < -ex {
			acx += BOUNDS_GAIN * (-ex - px)
		}
		if py > ey {
			acy -= BOUNDS_GAIN * (py - ey)
		}
		if py < -ey {
			acy += BOUNDS_GAIN * (-ey - py)
		}
		if pz > ez {
			acz -= BOUNDS_GAIN * (pz - ez)
		}
		if pz < -ez {
			acz += BOUNDS_GAIN * (-ez - pz)
		}

		// Predator flee: the same steering shape, weighted by the
		// separation weight and fading over the flee radius.
		if params.predator_active {
			fdx := px - params.predator[0]
			fdy := py - params.predator[1]
			fdz := pz - params.predator[2]
			fd2 := fdx*fdx + fdy*fdy + fdz*fdz
			if fd2 < flee_radius2 {
				gain := flee_gain * (1.0 - math.sqrt(fd2) / flee_radius)
				steer_rule(&acx, &acy, &acz, fdx, fdy, fdz, f.vx[i], f.vy[i], f.vz[i], params.max_speed, force, gain)
			}
		}

		// Clamp the total steering so no frame can inject an unbounded
		// impulse into an otherwise stable one.
		a2 := acx*acx + acy*acy + acz*acz
		if a2 > total_limit2 {
			scale := params.max_speed * STEERING_TOTAL_FACTOR / math.sqrt(a2)
			acx *= scale
			acy *= scale
			acz *= scale
		}
		f.ax[i] = acx
		f.ay[i] = acy
		f.az[i] = acz
	}

	// Pass two: integration. Semi-implicit Euler (velocity first, then
	// position from the new velocity), a max-speed clamp, and a final soft
	// position clamp into the bounds as the safety net behind the smooth
	// wall steering. Facing refreshes only above the epsilon, so an agent
	// that stops keeps pointing where it last moved.
	for i in 0..<f.count {
		vx := f.vx[i] + f.ax[i] * dt
		vy := f.vy[i] + f.ay[i] * dt
		vz := f.vz[i] + f.az[i] * dt
		speed2 := vx*vx + vy*vy + vz*vz
		max2 := params.max_speed * params.max_speed
		if speed2 > max2 {
			scale := params.max_speed / math.sqrt(speed2)
			vx *= scale
			vy *= scale
			vz *= scale
			speed2 = max2
		}
		f.px[i] = clamp_f32(f.px[i] + vx * dt, -f.params.half_extents.x, f.params.half_extents.x)
		f.py[i] = clamp_f32(f.py[i] + vy * dt, -f.params.half_extents.y, f.params.half_extents.y)
		f.pz[i] = clamp_f32(f.pz[i] + vz * dt, -f.params.half_extents.z, f.params.half_extents.z)
		f.vx[i] = vx
		f.vy[i] = vy
		f.vz[i] = vz
		if speed2 > FACING_EPSILON2 {
			inverse_speed := 1.0 / math.sqrt(speed2)
			f.fx[i] = vx * inverse_speed
			f.fy[i] = vy * inverse_speed
			f.fz[i] = vz * inverse_speed
		}
	}

	f.max_neighbors = max_neighbors
	f.frame += 1
	return true
}

// write_basis fills one instance's twelve floats in the engine's interleaved
// multimesh buffer layout (`mesh_storage.cpp`,
// `_multimesh_instance_get_transform`): basis row 0, then origin.x, basis
// row 1, then origin.y, basis row 2, then origin.z. The basis is built
// right-handed and orthonormal from the agent's unit facing: back = -facing
// (engine meshes face local -Z), side = world_up x back normalized with a
// (0, 0, 1) seed-up fallback when the agent flies nearly vertically, and up
// = back x side. The rotation stays orthonormal, so a downstream engine that
// renormalizes still reads back the same numbers.
write_basis :: proc(sink: []f32, px, py, pz, fx, fy, fz: f32) {
	bx, by, bz := -fx, -fy, -fz
	// side = (0, 1, 0) x back = (bz, 0, -bx)
	sx, sy, sz := bz, f32(0), -bx
	sl := math.sqrt(sx*sx + sy*sy + sz*sz)
	if sl < 1e-4 {
		// side = (0, 0, 1) x back = (-by, bx, 0)
		sx, sy, sz = -by, bx, f32(0)
		sl = math.sqrt(sx*sx + sy*sy + sz*sz)
	}
	sx /= sl
	sy /= sl
	sz /= sl
	// up = back x side completes the right-handed orthonormal triad
	// (side x up = back).
	ux := by * sz - bz * sy
	uy := bz * sx - bx * sz
	uz := bx * sy - by * sx
	sink[0] = sx
	sink[1] = ux
	sink[2] = bx
	sink[3] = px
	sink[4] = sy
	sink[5] = uy
	sink[6] = by
	sink[7] = py
	sink[8] = sz
	sink[9] = uz
	sink[10] = bz
	sink[11] = pz
}

// flock_write_transforms writes `count * FLOATS_PER_INSTANCE` floats into
// `sink` in the locked interleaved layout, one instance per agent: the basis
// from the agent's facing and the agent's position as the origin. The
// procedure touches exactly count*12 floats and allocates nothing.
flock_write_transforms :: proc(f: ^Flock, sink: []f32) -> bool {
	if f == nil || !f.live {
		return false
	}
	if i64(len(sink)) < i64(f.count)*FLOATS_PER_INSTANCE {
		return false
	}
	for index in 0..<f.count {
		write_basis(sink[index*FLOATS_PER_INSTANCE:index*FLOATS_PER_INSTANCE+FLOATS_PER_INSTANCE],
			f.px[index], f.py[index], f.pz[index],
			f.fx[index], f.fy[index], f.fz[index])
	}
	return true
}

// verify_instance_against recomputes instance `index`'s twelve interleaved
// floats from the flock's own current state and compares them bitwise to the
// caller's sink. This checks the layout writer against the state it claims
// to encode — it is deliberately NOT a physics oracle: the simulation is
// neighbor-coupled and stateful, so an independent oracle must replay the whole configuration/history.
// Small all-pairs tests check neighbor physics. Determinism is pinned by `state_digest` equality
// across runs; layout correctness is pinned by this check plus the engine
// readback path through `verify_readback`.
verify_instance_against :: proc(f: ^Flock, index: i32, sink: []f32) -> bool {
	if f == nil || !f.live {
		return false
	}
	if index < 0 || index >= f.count {
		return false
	}
	if len(sink) < FLOATS_PER_INSTANCE {
		return false
	}
	expected: [FLOATS_PER_INSTANCE]f32
	write_basis(expected[:], f.px[index], f.py[index], f.pz[index], f.fx[index], f.fy[index], f.fz[index])
	for component in 0..<FLOATS_PER_INSTANCE {
		if transmute(u32)sink[component] != transmute(u32)expected[component] {
			return false
		}
	}
	return true
}

// verify_readback permutes the engine getter's Transform3D memory shape —
// basis rows 0..2 then origin, exactly what
// `rendering_server_multimesh_instance_get_transform` returns after the
// server decodes its interleaved storage — back into the written buffer
// order, and compares it bitwise against the caller's written window. It
// validates the engine read path against exactly the floats this package
// wrote, with no simulation state involved, so it stays valid even after the
// flock has stepped on.
verify_readback :: proc(written: []f32, stored: []f32) -> bool {
	if len(written) < FLOATS_PER_INSTANCE || len(stored) < FLOATS_PER_INSTANCE {
		return false
	}
	permuted: [FLOATS_PER_INSTANCE]f32
	permuted[0] = stored[0]
	permuted[1] = stored[1]
	permuted[2] = stored[2]
	permuted[3] = stored[9] // origin.x
	permuted[4] = stored[3]
	permuted[5] = stored[4]
	permuted[6] = stored[5]
	permuted[7] = stored[10] // origin.y
	permuted[8] = stored[6]
	permuted[9] = stored[7]
	permuted[10] = stored[8]
	permuted[11] = stored[11] // origin.z
	for component in 0..<FLOATS_PER_INSTANCE {
		if transmute(u32)written[component] != transmute(u32)permuted[component] {
			return false
		}
	}
	return true
}

// clone duplicates the whole simulation into a fresh arena: same count,
// same parameters, same every float. Test-only; used to replay an identical
// history beside the original and compare digests.
clone :: proc(f: ^Flock) -> (^Flock, bool) {
	if f == nil || !f.live {
		return nil, false
	}
	copied, ok := flock_create(f.count, 0, f.params)
	if !ok {
		return nil, false
	}
	copied.frame = f.frame
	copied.max_neighbors = f.max_neighbors
	copied.cell_size = f.cell_size
	copied.grid_x = f.grid_x
	copied.grid_y = f.grid_y
	copied.grid_z = f.grid_z
	copied.cells = f.cells
	copy(copied.px, f.px)
	copy(copied.py, f.py)
	copy(copied.pz, f.pz)
	copy(copied.vx, f.vx)
	copy(copied.vy, f.vy)
	copy(copied.vz, f.vz)
	copy(copied.fx, f.fx)
	copy(copied.fy, f.fy)
	copy(copied.fz, f.fz)
	return copied, true
}

// state_digest folds every position, velocity, and facing into one stable
// 64-bit value (FNV-1a over the f32 bit patterns) so a whole simulation
// state can be compared across runs, transports, and flocks without holding
// the arrays. Same seed, parameters, and step count fold to the same
// digest; flipping one seed bit does not.
state_digest :: proc(f: ^Flock) -> u64 {
	if f == nil || !f.live {
		return 0
	}
	hash := u64(0xcbf29ce484222325)
	fold_u64(&hash, u64(f.count))
	for array in ([9][]f32{f.px, f.py, f.pz, f.vx, f.vy, f.vz, f.fx, f.fy, f.fz}) {
		for value in array {
			fold_u32(&hash, transmute(u32) value)
		}
	}
	return hash
}

fold_u32 :: proc(hash: ^u64, value: u32) {
	for shift in 0..<4 {
		hash^ = (hash^ ~ u64((value >> (u32(shift) * 8)) & 0xFF)) * 0x100000001b3
	}
}

fold_u64 :: proc(hash: ^u64, value: u64) {
	for shift in 0..<8 {
		hash^ = (hash^ ~ u64((value >> (u32(shift) * 8)) & 0xFF)) * 0x100000001b3
	}
}

// max_neighbors_observed reports the largest per-agent neighbor count seen
// in the most recent step — the workload's measured spatial-hash bound. It
// is zero before the first step.
max_neighbors_observed :: proc(f: ^Flock) -> i32 {
	if f == nil || !f.live {
		return 0
	}
	return f.max_neighbors
}

// flock_count reports the agent count.
flock_count :: proc(f: ^Flock) -> i32 {
	if f == nil {
		return 0
	}
	return f.count
}

// flock_frame reports the number of completed steps.
flock_frame :: proc(f: ^Flock) -> i64 {
	if f == nil {
		return 0
	}
	return f.frame
}