#!/usr/bin/env python3
"""Deterministic, standard-library reference oracle for the SmoothLife contract.

The implementation deliberately favors explicit scalar loops and fixed operation order over
speed. Model-semantic arrays use Python binary64. The retained packed-unitary FFT path uses
``f32`` after every arithmetic operation so its stages can be compared bit-for-bit.
"""

from __future__ import annotations

import math
import struct
from dataclasses import dataclass
from typing import Any, Callable, Iterable, Iterator, Sequence, TypedDict

TAU = 2.0 * math.pi
U32_MASK = (1 << 32) - 1
U64_MODULUS = 1 << 64


def clamp(value: float, low: float = 0.0, high: float = 1.0) -> float:
    return min(high, max(low, value))


def product(values: Iterable[int]) -> int:
    result = 1
    for value in values:
        result *= value
    return result


def flat_index(coords: Sequence[int], shape: Sequence[int]) -> int:
    """Canonical x-fast row-major index."""
    if len(coords) != len(shape) or not shape:
        raise ValueError("coordinates and non-empty shape must have equal rank")
    stride = 1
    index = 0
    for coord, extent in zip(coords, shape, strict=True):
        if extent <= 0 or not 0 <= coord < extent:
            raise IndexError((tuple(coords), tuple(shape)))
        index += coord * stride
        stride *= extent
    return index


def unflatten_index(index: int, shape: Sequence[int]) -> tuple[int, ...]:
    if not 0 <= index < product(shape):
        raise IndexError((index, tuple(shape)))
    coords: list[int] = []
    for extent in shape:
        coords.append(index % extent)
        index //= extent
    return tuple(coords)


def iter_coords(shape: Sequence[int]) -> Iterator[tuple[int, ...]]:
    for index in range(product(shape)):
        yield unflatten_index(index, shape)


def signed_offset(index: int, extent: int) -> int:
    """Canonical periodic offset. The contract requires this for even extents."""
    if not 0 <= index < extent:
        raise IndexError((index, extent))
    return index if index < extent / 2 else index - extent


def wrap_coords(coords: Sequence[int], shape: Sequence[int]) -> tuple[int, ...]:
    return tuple(coord % extent for coord, extent in zip(coords, shape, strict=True))


# Shared scalar rules -------------------------------------------------------


def transition_l(r: float, center: float, width: float) -> float:
    if width <= 0.0:
        raise ValueError("transition width must be positive")
    if r < center - width / 2.0:
        return 0.0
    if r > center + width / 2.0:
        return 1.0
    return (r - center) / width + 0.5


def logistic(x: float, center: float, width: float) -> float:
    if width <= 0.0:
        raise ValueError("curve width must be positive")
    z = 4.0 * (x - center) / width
    if z >= 0.0:
        return 1.0 / (1.0 + math.exp(-z))
    ez = math.exp(z)
    return ez / (1.0 + ez)


def rising_curve(curve_type: int, x: float, center: float, width: float) -> float:
    if curve_type == 0:
        return 1.0 if x >= center else 0.0
    if width <= 0.0:
        raise ValueError("curve width must be positive")
    if curve_type in (1, 2, 3):
        if x < center - width / 2.0:
            return 0.0
        if x > center + width / 2.0:
            return 1.0
        u = (x - center + width / 2.0) / width
        if curve_type == 1:
            return u
        if curve_type == 2:
            return u * u * (3.0 - 2.0 * u)
        return 0.5 * math.sin((TAU / 2.0) * (x - center) / width) + 0.5
    if curve_type == 4:
        return logistic(x, center, width)
    if curve_type == 5:
        return math.atan((x - center) * (TAU / 2.0) / width) / (TAU / 2.0) + 0.5
    if curve_type == 6:
        return (
            1.1
            * math.atan((x - center) / width)
            / (TAU / 4.0)
            * math.cos(1.4 * (x - center))
            + 1.0
        ) / 2.0
    if curve_type == 7:
        return (logistic(x, center, width) - 0.5) * (
            1.0 + math.exp(-((x - center) ** 2) / (width * width))
        ) + 0.5
    raise ValueError(f"invalid rising curve type {curve_type}")


def window_curve(window_type: int, n: float, a: float, b: float, width: float) -> float:
    if not 0 <= window_type <= 9:
        raise ValueError(f"invalid window type {window_type}")
    if window_type <= 7:
        return rising_curve(window_type, n, a, width) * (
            1.0 - rising_curve(window_type, n, b, width)
        )
    base = logistic(n, a, width) * (1.0 - logistic(n, b, width))
    mid = (a + b) / 2.0
    notch = 0.2 * math.exp(-((20.0 * (n - mid)) ** 2))
    return base * (1.0 - notch if window_type == 8 else 1.0 + notch)


def mix_curve(
    mix_type: int, left: float, right: float, m: float, width: float
) -> float:
    q = rising_curve(mix_type, m, 0.5, width)
    return left * (1.0 - q) + right * q


@dataclass(frozen=True)
class Rule:
    b1: float = 0.278
    b2: float = 0.365
    d1: float = 0.267
    d2: float = 0.445
    sn: float = 0.065
    sm: float = 0.110
    sigmode: int = 2
    sigtype: int = 4
    mixtype: int = 4


def rule_target(n: float, m: float, rule: Rule) -> float:
    if not 1 <= rule.sigmode <= 4:
        raise ValueError(f"invalid sigmode {rule.sigmode}")
    if not 0 <= rule.sigtype <= 9 or not 0 <= rule.mixtype <= 7:
        raise ValueError("invalid rule curve enum")
    birth = window_curve(rule.sigtype, n, rule.b1, rule.b2, rule.sn)
    death = window_curve(rule.sigtype, n, rule.d1, rule.d2, rule.sn)
    q = rising_curve(rule.mixtype, m, 0.5, rule.sm)
    if rule.sigmode == 1:
        return birth * (1.0 - m) + death * m
    if rule.sigmode == 2:
        return birth * (1.0 - q) + death * q
    if rule.sigmode == 3:
        low = rule.b1 * (1.0 - m) + rule.d1 * m
        high = rule.b2 * (1.0 - m) + rule.d2 * m
        return window_curve(rule.sigtype, n, low, high, rule.sn)
    low = rule.b1 * (1.0 - q) + rule.d1 * q
    high = rule.b2 * (1.0 - q) + rule.d2 * q
    return window_curve(rule.sigtype, n, low, high, rule.sn)


# Sampled periodic kernels -------------------------------------------------


def sampled_kernels(
    shape: Sequence[int], ra: float, rr: float, rb: float
) -> dict[str, Any]:
    if ra <= 0.0 or rr <= 0.0 or rb <= 0.0:
        raise ValueError("kernel geometry must be positive")
    ri = ra / rr
    width = ra / rb
    disk: list[float] = []
    ring: list[float] = []
    support: list[dict[str, Any]] = []
    for coords in iter_coords(shape):
        offsets = tuple(
            signed_offset(coord, extent)
            for coord, extent in zip(coords, shape, strict=True)
        )
        radius = math.sqrt(sum(offset * offset for offset in offsets))
        inner_l = transition_l(radius, ri, width)
        kd = 1.0 - inner_l
        kr = inner_l * (1.0 - transition_l(radius, ra, width))
        disk.append(kd)
        ring.append(kr)
        if kd != 0.0 or kr != 0.0:
            support.append(
                {
                    "coords": list(coords),
                    "offset": list(offsets),
                    "radius": radius,
                    "disk": kd,
                    "ring": kr,
                }
            )
    disk_sum = sum(disk)
    ring_sum = sum(ring)
    if disk_sum == 0.0 or ring_sum == 0.0:
        raise ValueError("sampled kernel has zero normalization")
    warnings = []
    for axis, extent in enumerate(shape):
        if extent % 2 == 0 and any(
            abs(entry["offset"][axis]) == extent // 2
            and (entry["disk"] > 0.0 or entry["ring"] > 0.0)
            for entry in support
        ):
            warnings.append(
                f"nonzero support touches periodic Nyquist offset on axis {axis}"
            )
    if len(support) == product(shape):
        warnings.append("nonzero disk/ring support covers every periodic offset")
    return {
        "shape": list(shape),
        "ra": ra,
        "rr": rr,
        "rb": rb,
        "ri": ri,
        "width": width,
        "raw_disk": disk,
        "raw_ring": ring,
        "disk_sum": disk_sum,
        "ring_sum": ring_sum,
        "normalized_disk": [value / disk_sum for value in disk],
        "normalized_ring": [value / ring_sum for value in ring],
        "support": support,
        "warnings": warnings,
    }


def circular_convolution(
    field: Sequence[float], kernel: Sequence[float], shape: Sequence[int]
) -> list[float]:
    if len(field) != product(shape) or len(kernel) != product(shape):
        raise ValueError("field and kernel must match shape")
    result: list[float] = []
    for output in iter_coords(shape):
        total = 0.0
        for kernel_coords in iter_coords(shape):
            source = tuple(
                (out_coord - kernel_coord) % extent
                for out_coord, kernel_coord, extent in zip(
                    output, kernel_coords, shape, strict=True
                )
            )
            total += (
                field[flat_index(source, shape)]
                * kernel[flat_index(kernel_coords, shape)]
            )
        result.append(total)
    return result


def neighborhoods(
    field: Sequence[float], shape: Sequence[int], ra: float, rr: float, rb: float
) -> tuple[list[float], list[float], dict[str, Any]]:
    kernels = sampled_kernels(shape, ra, rr, rb)
    m = [
        value / kernels["disk_sum"]
        for value in circular_convolution(field, kernels["raw_disk"], shape)
    ]
    n = [
        value / kernels["ring_sum"]
        for value in circular_convolution(field, kernels["raw_ring"], shape)
    ]
    return m, n, kernels


# Exact ChaCha12 stream ----------------------------------------------------


def _rotl32(value: int, count: int) -> int:
    return ((value << count) & U32_MASK) | (value >> (32 - count))


def _quarter_round(state: list[int], a: int, b: int, c: int, d: int) -> None:
    state[a] = (state[a] + state[b]) & U32_MASK
    state[d] = _rotl32(state[d] ^ state[a], 16)
    state[c] = (state[c] + state[d]) & U32_MASK
    state[b] = _rotl32(state[b] ^ state[c], 12)
    state[a] = (state[a] + state[b]) & U32_MASK
    state[d] = _rotl32(state[d] ^ state[a], 8)
    state[c] = (state[c] + state[d]) & U32_MASK
    state[b] = _rotl32(state[b] ^ state[c], 7)


class ChaCha12:
    """ChaCha with 256-bit key, original 64-bit counter, and 64-bit stream.

    The key is little-endian ``seed: u64`` followed by 24 zero bytes. Counter and stream
    start at zero. ``u64`` consumes low then high sequential ``u32`` words.
    """

    def __init__(self, seed: int):
        if not 0 <= seed < U64_MODULUS:
            raise ValueError("seed must fit u64")
        self._key = [seed & U32_MASK, (seed >> 32) & U32_MASK] + [0] * 6
        self._counter = 0
        self._stream = 0
        self._buffer: list[int] = []
        self._position = 0

    def _block(self) -> list[int]:
        initial = [
            0x61707865,
            0x3320646E,
            0x79622D32,
            0x6B206574,
            *self._key,
            self._counter & U32_MASK,
            (self._counter >> 32) & U32_MASK,
            self._stream & U32_MASK,
            (self._stream >> 32) & U32_MASK,
        ]
        state = initial.copy()
        for _ in range(6):
            _quarter_round(state, 0, 4, 8, 12)
            _quarter_round(state, 1, 5, 9, 13)
            _quarter_round(state, 2, 6, 10, 14)
            _quarter_round(state, 3, 7, 11, 15)
            _quarter_round(state, 0, 5, 10, 15)
            _quarter_round(state, 1, 6, 11, 12)
            _quarter_round(state, 2, 7, 8, 13)
            _quarter_round(state, 3, 4, 9, 14)
        output = [
            (value + origin) & U32_MASK
            for value, origin in zip(state, initial, strict=True)
        ]
        self._counter = (self._counter + 1) & ((1 << 64) - 1)
        return output

    def u32(self) -> int:
        if self._position == len(self._buffer):
            self._buffer = self._block()
            self._position = 0
        result = self._buffer[self._position]
        self._position += 1
        return result

    def u64(self) -> int:
        low = self.u32()
        high = self.u32()
        return low | (high << 32)

    def float53(self) -> float:
        return (self.u64() >> 11) * (1.0 / (1 << 53))

    def integer(self, low: int, high: int) -> int:
        """Uniform integer in [low, high), rejecting the incomplete high residue."""
        if high <= low:
            raise ValueError("empty integer range")
        span = high - low
        if span > U64_MODULUS:
            raise ValueError("range is wider than u64")
        limit = U64_MODULUS - (U64_MODULUS % span)
        while True:
            value = self.u64()
            if value < limit:
                return low + value % span


# Deterministic initializers -----------------------------------------------


def _periodic_continuous_distance(sample: int, center: float, extent: int) -> float:
    distance = abs(sample - center)
    return min(distance, extent - distance)


def planar_splats(shape: Sequence[int], ra: float, seed: int) -> dict[str, Any]:
    if not 1 <= len(shape) <= 3 or ra <= 0.0:
        raise ValueError("planar splats require 1D, 2D, or 3D positive geometry")
    denominator = math.prod(min(2.0 * ra, extent) for extent in shape)
    count = math.floor(product(shape) / denominator) + 1
    rng = ChaCha12(seed)
    field = [0.0] * product(shape)
    splats: list[dict[str, Any]] = []
    for _ in range(count):
        center = [rng.float53() * extent for extent in shape]
        radius = (0.5 + 0.5 * rng.float53()) * ra
        splats.append({"center": center, "radius": radius})
        for coords in iter_coords(shape):
            distance = math.sqrt(
                sum(
                    _periodic_continuous_distance(coord, axis_center, extent) ** 2
                    for coord, axis_center, extent in zip(
                        coords, center, shape, strict=True
                    )
                )
            )
            if distance < radius:
                field[flat_index(coords, shape)] = 1.0
    return {
        "shape": list(shape),
        "ra": ra,
        "seed": seed,
        "count": count,
        "draw_order": [
            *[f"center_{axis}_float53" for axis in "xyz"[: len(shape)]],
            "radius_float53",
        ],
        "radius_interval": "[0.5*ra, ra)",
        "paint": "integer lattice samples with periodic Euclidean distance strictly less than radius",
        "splats": splats,
        "field": field,
    }


# Packed-unitary IEEE-f32 FFT ---------------------------------------------


def f32(value: float) -> float:
    try:
        return struct.unpack("<f", struct.pack("<f", value))[0]
    except (OverflowError, struct.error) as error:
        raise ValueError(
            f"value cannot be represented as IEEE binary32: {value!r}"
        ) from error


def f32_bits(value: float) -> int:
    try:
        return struct.unpack("<I", struct.pack("<f", f32(value)))[0]
    except struct.error as error:
        raise ValueError(
            f"value cannot be encoded as IEEE binary32: {value!r}"
        ) from error


def fadd(left: float, right: float) -> float:
    return f32(f32(left) + f32(right))


def fsub(left: float, right: float) -> float:
    return f32(f32(left) - f32(right))


def fmul(left: float, right: float) -> float:
    return f32(f32(left) * f32(right))


def _cmul(left: tuple[float, float], right: tuple[float, float]) -> tuple[float, float]:
    real = fsub(fmul(left[0], right[0]), fmul(left[1], right[1]))
    imag = fadd(fmul(left[0], right[1]), fmul(left[1], right[0]))
    return real, imag


def bit_reverse(value: int, bits: int) -> int:
    result = 0
    for bit in range(bits):
        result = (result << 1) | ((value >> bit) & 1)
    return result


def _power_bits(extent: int) -> int:
    if extent <= 0 or extent & (extent - 1):
        raise ValueError(f"extent {extent} is not a positive power of two")
    return extent.bit_length() - 1


class LegacyPlanEntry(TypedDict):
    output: int
    source: list[int]
    twiddle: list[float]


def legacy_butterfly_plan(length: int, stage: int, sign: int) -> list[LegacyPlanEntry]:
    bits = _power_bits(length)
    if not 1 <= stage <= bits or sign not in (-1, 1):
        raise ValueError("invalid butterfly plan")
    span = 1 << stage
    entries: list[LegacyPlanEntry] = []
    for output in range(length):
        j = output % span
        if j < span // 2:
            source_a, source_b = output, output + span // 2
        else:
            source_a, source_b = output - span // 2, output
        if stage == 1:
            source_a = bit_reverse(source_a, bits)
            source_b = bit_reverse(source_b, bits)
        angle = sign * TAU * j / span
        entries.append(
            {
                "output": output,
                "source": [source_a, source_b],
                "twiddle": [f32(math.cos(angle)), f32(math.sin(angle))],
            }
        )
    return entries


def legacy_x_conversion_plan(nx: int, sign: int) -> list[LegacyPlanEntry]:
    _power_bits(nx)
    half = nx // 2
    entries: list[LegacyPlanEntry] = []
    for output in range(half + 1):
        if sign == -1 and output in (0, half):
            source_a = source_b = 0
        else:
            source_a, source_b = output, half - output
        angle = sign * TAU * (output / nx + 0.25)
        entries.append(
            {
                "output": output,
                "source": [source_a, source_b],
                "twiddle": [f32(math.cos(angle)), f32(math.sin(angle))],
            }
        )
    return entries


def _packed_shape(shape: Sequence[int]) -> tuple[int, ...]:
    return (shape[0] // 2 + 1, *shape[1:])


def legacy_pack_real(
    field: Sequence[float], shape: Sequence[int]
) -> list[tuple[float, float]]:
    nx = shape[0]
    _power_bits(nx)
    if len(field) != product(shape):
        raise ValueError("field must match shape")
    packed_shape = _packed_shape(shape)
    result = [(f32(0.0), f32(0.0)) for _ in range(product(packed_shape))]
    for coords in iter_coords(shape[1:]):
        for packed_x in range(nx // 2):
            even = (2 * packed_x, *coords)
            odd = (2 * packed_x + 1, *coords)
            destination = (packed_x, *coords)
            result[flat_index(destination, packed_shape)] = (
                f32(field[flat_index(even, shape)]),
                f32(field[flat_index(odd, shape)]),
            )
    return result


def legacy_unpack_real(
    packed: Sequence[tuple[float, float]], shape: Sequence[int]
) -> list[float]:
    nx = shape[0]
    packed_shape = _packed_shape(shape)
    result = [f32(0.0)] * product(shape)
    for coords in iter_coords(shape[1:]):
        for packed_x in range(nx // 2):
            value = packed[flat_index((packed_x, *coords), packed_shape)]
            result[flat_index((2 * packed_x, *coords), shape)] = f32(value[0])
            result[flat_index((2 * packed_x + 1, *coords), shape)] = f32(value[1])
    return result


def _legacy_axis_stage(
    values: Sequence[tuple[float, float]],
    packed_shape: Sequence[int],
    axis: int,
    plan: Sequence[LegacyPlanEntry],
    active_length: int | None = None,
) -> list[tuple[float, float]]:
    active_length = active_length or packed_shape[axis]
    scale = f32(1.0 / math.sqrt(2.0))
    output = [(f32(0.0), f32(0.0)) for _ in values]
    for coords in iter_coords(packed_shape):
        axis_coord = coords[axis]
        if axis_coord >= active_length:
            continue
        entry = plan[axis_coord]
        source_pair = entry["source"]
        twiddle_values = entry["twiddle"]
        a_coords = list(coords)
        b_coords = list(coords)
        a_coords[axis] = source_pair[0]
        b_coords[axis] = source_pair[1]
        a = values[flat_index(a_coords, packed_shape)]
        b = values[flat_index(b_coords, packed_shape)]
        twiddle = (twiddle_values[0], twiddle_values[1])
        # Matches: (a.r + cos*b.r - sin*b.i) / sqrt(2), then the imag expression.
        real = fsub(fadd(a[0], fmul(twiddle[0], b[0])), fmul(twiddle[1], b[1]))
        imag = fadd(fadd(a[1], fmul(twiddle[0], b[1])), fmul(twiddle[1], b[0]))
        output[flat_index(coords, packed_shape)] = (
            fmul(real, scale),
            fmul(imag, scale),
        )
    return output


def _legacy_x_conversion(
    values: Sequence[tuple[float, float]], shape: Sequence[int], sign: int
) -> list[tuple[float, float]]:
    packed_shape = _packed_shape(shape)
    nx = shape[0]
    plan = legacy_x_conversion_plan(nx, sign)
    active_length = nx // 2 + 1 if sign == -1 else nx // 2
    scale = f32(0.5 / math.sqrt(2.0) if sign == -1 else 0.5 * math.sqrt(2.0))
    output = [(f32(0.0), f32(0.0)) for _ in values]
    for coords in iter_coords(packed_shape):
        x = coords[0]
        if x >= active_length:
            continue
        entry = plan[x]
        source_pair = entry["source"]
        twiddle_values = entry["twiddle"]
        a_coords = (source_pair[0], *coords[1:])
        b_coords = (source_pair[1], *coords[1:])
        a = values[flat_index(a_coords, packed_shape)]
        raw_b = values[flat_index(b_coords, packed_shape)]
        b = (raw_b[0], f32(-raw_b[1]))
        sum_value = (fadd(a[0], b[0]), fadd(a[1], b[1]))
        difference = (fsub(a[0], b[0]), fsub(a[1], b[1]))
        twiddle = (twiddle_values[0], twiddle_values[1])
        rotated = _cmul(difference, twiddle)
        value = (
            fmul(fadd(sum_value[0], rotated[0]), scale),
            fmul(fadd(sum_value[1], rotated[1]), scale),
        )
        output[flat_index(coords, packed_shape)] = value
    return output


def legacy_forward(
    field: Sequence[float], shape: Sequence[int], capture_stages: bool = False
) -> tuple[list[tuple[float, float]], list[dict[str, Any]]]:
    if not 1 <= len(shape) <= 3:
        raise ValueError("legacy FFT rank must be 1, 2, or 3")
    for extent in shape:
        _power_bits(extent)
    packed_shape = _packed_shape(shape)
    values = legacy_pack_real(field, shape)
    stages: list[dict[str, Any]] = [
        {"name": "adjacent_real_pack", "values": values.copy()}
    ]
    x_length = shape[0] // 2
    for stage in range(1, _power_bits(shape[0])):
        values = _legacy_axis_stage(
            values,
            packed_shape,
            0,
            legacy_butterfly_plan(x_length, stage, -1),
            x_length,
        )
        stages.append({"name": f"x_butterfly_{stage}", "values": values.copy()})
    values = _legacy_x_conversion(values, shape, -1)
    stages.append({"name": "x_real_to_half_complex", "values": values.copy()})
    for axis in range(1, len(shape)):
        for stage in range(1, _power_bits(shape[axis]) + 1):
            values = _legacy_axis_stage(
                values,
                packed_shape,
                axis,
                legacy_butterfly_plan(shape[axis], stage, -1),
            )
            stages.append(
                {"name": f"{'xyz'[axis]}_butterfly_{stage}", "values": values.copy()}
            )
    return values, stages if capture_stages else []


def legacy_inverse(
    spectrum: Sequence[tuple[float, float]],
    shape: Sequence[int],
    capture_stages: bool = False,
) -> tuple[list[float], list[dict[str, Any]]]:
    packed_shape = _packed_shape(shape)
    if len(spectrum) != product(packed_shape):
        raise ValueError("spectrum must match packed shape")
    values = [(f32(real), f32(imag)) for real, imag in spectrum]
    stages: list[dict[str, Any]] = []
    for axis in range(len(shape) - 1, 0, -1):
        for stage in range(1, _power_bits(shape[axis]) + 1):
            values = _legacy_axis_stage(
                values, packed_shape, axis, legacy_butterfly_plan(shape[axis], stage, 1)
            )
            stages.append(
                {
                    "name": f"{'xyz'[axis]}_inverse_butterfly_{stage}",
                    "values": values.copy(),
                }
            )
    values = _legacy_x_conversion(values, shape, 1)
    stages.append({"name": "x_half_complex_to_packed_real", "values": values.copy()})
    x_length = shape[0] // 2
    for stage in range(1, _power_bits(shape[0])):
        values = _legacy_axis_stage(
            values, packed_shape, 0, legacy_butterfly_plan(x_length, stage, 1), x_length
        )
        stages.append({"name": f"x_inverse_butterfly_{stage}", "values": values.copy()})
    real = legacy_unpack_real(values, shape)
    stages.append({"name": "adjacent_real_unpack", "values": real.copy()})
    return real, stages if capture_stages else []


def legacy_spectral_product(
    field_spectrum: Sequence[tuple[float, float]],
    kernel_spectrum: Sequence[tuple[float, float]],
    sample_count: int,
    kernel_sum: float,
) -> list[tuple[float, float]]:
    if len(field_spectrum) != len(kernel_spectrum) or kernel_sum == 0.0:
        raise ValueError("incompatible spectra or zero kernel sum")
    correction = f32(math.sqrt(sample_count) / kernel_sum)
    result: list[tuple[float, float]] = []
    for field_value, kernel_value in zip(field_spectrum, kernel_spectrum, strict=True):
        scaled_kernel = (
            fmul(kernel_value[0], correction),
            fmul(kernel_value[1], correction),
        )
        result.append(_cmul(field_value, scaled_kernel))
    return result


def legacy_convolution(
    field: Sequence[float],
    kernel: Sequence[float],
    shape: Sequence[int],
    capture_stages: bool = False,
) -> dict[str, Any]:
    field_spectrum, field_stages = legacy_forward(field, shape, capture_stages)
    kernel_spectrum, kernel_stages = legacy_forward(kernel, shape, capture_stages)
    spectral_product = legacy_spectral_product(
        field_spectrum, kernel_spectrum, product(shape), sum(kernel)
    )
    output, inverse_stages = legacy_inverse(spectral_product, shape, capture_stages)
    direct = [
        value / sum(kernel) for value in circular_convolution(field, kernel, shape)
    ]
    return {
        "field_spectrum": field_spectrum,
        "kernel_spectrum": kernel_spectrum,
        "spectral_product": spectral_product,
        "output": output,
        "direct_normalized": direct,
        "field_stages": field_stages,
        "kernel_stages": kernel_stages,
        "inverse_stages": inverse_stages,
        "correction": f32(math.sqrt(product(shape)) / sum(kernel)),
    }


# Integrators ---------------------------------------------------------------

VectorDerivative = Callable[[Sequence[float], float], Sequence[float]]


def _clamped_axpy(
    origin: Sequence[float], scale: float, delta: Sequence[float]
) -> list[float]:
    return [
        clamp(value + scale * change)
        for value, change in zip(origin, delta, strict=True)
    ]


def euler_step(
    state: Sequence[float], time: float, dt: float, derivative: VectorDerivative
) -> dict[str, Any]:
    k1 = list(derivative(state, time))
    return {"k1": k1, "next": _clamped_axpy(state, dt, k1)}


def ab3_startup(
    initial: Sequence[float], time: float, dt: float, derivative: VectorDerivative
) -> dict[str, Any]:
    states = [list(initial)]
    derivatives: list[list[float]] = [list(derivative(states[0], time))]
    methods = ["Euler"]
    states.append(_clamped_axpy(states[-1], dt, derivatives[-1]))
    derivatives.append(list(derivative(states[-1], time + dt)))
    methods.append("AB2")
    ab2 = [
        (3.0 * current - previous) / 2.0
        for current, previous in zip(derivatives[-1], derivatives[-2], strict=True)
    ]
    states.append(_clamped_axpy(states[-1], dt, ab2))
    derivatives.append(list(derivative(states[-1], time + 2.0 * dt)))
    methods.append("AB3")
    ab3 = [
        (23.0 * current - 16.0 * previous + 5.0 * oldest) / 12.0
        for current, previous, oldest in zip(
            derivatives[-1], derivatives[-2], derivatives[-3], strict=True
        )
    ]
    states.append(_clamped_axpy(states[-1], dt, ab3))
    return {
        "methods": methods,
        "states": states,
        "derivatives": derivatives,
        "ab2_combination": ab2,
        "ab3_combination": ab3,
    }


def rk4_step(
    state: Sequence[float], time: float, dt: float, derivative: VectorDerivative
) -> dict[str, Any]:
    origin = list(state)
    k1 = list(derivative(origin, time))
    stage2 = _clamped_axpy(origin, dt / 2.0, k1)
    k2 = list(derivative(stage2, time + dt / 2.0))
    stage3 = _clamped_axpy(origin, dt / 2.0, k2)
    k3 = list(derivative(stage3, time + dt / 2.0))
    stage4 = _clamped_axpy(origin, dt, k3)
    k4 = list(derivative(stage4, time + dt))
    combined = [
        (first + 2.0 * second + 2.0 * third + fourth) / 6.0
        for first, second, third, fourth in zip(k1, k2, k3, k4, strict=True)
    ]
    return {
        "stage_states": [origin, stage2, stage3, stage4],
        "derivatives": [k1, k2, k3, k4],
        "combined": combined,
        "next": _clamped_axpy(origin, dt, combined),
    }


TargetFunction = Callable[[Sequence[float]], Sequence[float]]


def rk4_relaxation(
    state: Sequence[float], dt: float, target: TargetFunction, reference: str
) -> dict[str, Any]:
    if reference not in ("stage_state", "step_origin"):
        raise ValueError("unknown RK4 relaxation reference")
    origin = list(state)
    target1 = list(target(origin))
    k1 = [desired - value for desired, value in zip(target1, origin, strict=True)]
    stage2 = _clamped_axpy(origin, dt / 2.0, k1)
    target2 = list(target(stage2))
    base2 = stage2 if reference == "stage_state" else origin
    k2 = [desired - value for desired, value in zip(target2, base2, strict=True)]
    stage3 = _clamped_axpy(origin, dt / 2.0, k2)
    target3 = list(target(stage3))
    base3 = stage3 if reference == "stage_state" else origin
    k3 = [desired - value for desired, value in zip(target3, base3, strict=True)]
    stage4 = _clamped_axpy(origin, dt, k3)
    target4 = list(target(stage4))
    base4 = stage4 if reference == "stage_state" else origin
    k4 = [desired - value for desired, value in zip(target4, base4, strict=True)]
    combined = [
        (first + 2.0 * second + 2.0 * third + fourth) / 6.0
        for first, second, third, fourth in zip(k1, k2, k3, k4, strict=True)
    ]
    return {
        "reference": reference,
        "stage_states": [origin, stage2, stage3, stage4],
        "targets": [target1, target2, target3, target4],
        "derivatives": [k1, k2, k3, k4],
        "combined": combined,
        "next": _clamped_axpy(origin, dt, combined),
    }


# Corrected multiscale ------------------------------------------------------


@dataclass(frozen=True)
class Scale:
    ra: float
    rr: float
    rb: float
    dt: float
    rule: Rule


def _scale_input(
    state: Sequence[float],
    shape: Sequence[int],
    scales: Sequence[Scale],
    index: int,
    interpretation: str,
) -> tuple[list[float], list[float]]:
    _, ring_i, _ = neighborhoods(
        state, shape, scales[index].ra, scales[index].rr, scales[index].rb
    )
    if interpretation == "independent":
        disk_i, _, _ = neighborhoods(
            state, shape, scales[index].ra, scales[index].rr, scales[index].rb
        )
        return ring_i, disk_i
    if interpretation != "chained":
        raise ValueError("unknown kernel interpretation")
    if index < len(scales) - 1:
        _, next_ring, _ = neighborhoods(
            state,
            shape,
            scales[index + 1].ra,
            scales[index + 1].rr,
            scales[index + 1].rb,
        )
        return ring_i, next_ring
    disk_i, _, _ = neighborhoods(
        state, shape, scales[index].ra, scales[index].rr, scales[index].rb
    )
    return ring_i, disk_i


def multiscale_step(
    initial: Sequence[float],
    shape: Sequence[int],
    scales: Sequence[Scale],
    interpretation: str,
    composition: str,
    dynamics: str,
) -> dict[str, Any]:
    if len(scales) != 3 or dynamics not in ("growth", "relaxation"):
        raise ValueError("multiscale requires three growth/relaxation scales")

    def response(
        state: Sequence[float], reference: Sequence[float], index: int
    ) -> tuple[list[float], list[float]]:
        n, m = _scale_input(state, shape, scales, index, interpretation)
        target = [
            rule_target(nv, mv, scales[index].rule) for nv, mv in zip(n, m, strict=True)
        ]
        if dynamics == "growth":
            increment = [scales[index].dt * (2.0 * value - 1.0) for value in target]
        else:
            increment = [
                scales[index].dt * (value - base)
                for value, base in zip(target, reference, strict=True)
            ]
        return target, increment

    targets: list[list[float]] = []
    increments: list[list[float]] = []
    clamp_stages: list[list[float]] = []
    if composition == "sequential":
        state = list(initial)
        for index in range(3):
            target, increment = response(state, state, index)
            targets.append(target)
            increments.append(increment)
            state = [
                clamp(value + change)
                for value, change in zip(state, increment, strict=True)
            ]
            clamp_stages.append(state.copy())
        result = state
    else:
        for index in range(3):
            target, increment = response(initial, initial, index)
            targets.append(target)
            increments.append(increment)
        if composition == "ordered_clamped_sum":
            result = list(initial)
            for increment in increments:
                result = [
                    clamp(value + change)
                    for value, change in zip(result, increment, strict=True)
                ]
                clamp_stages.append(result.copy())
        elif composition == "mean_increment":
            combined = [sum(changes) / 3.0 for changes in zip(*increments, strict=True)]
            result = [
                clamp(value + change)
                for value, change in zip(initial, combined, strict=True)
            ]
            clamp_stages.append(result.copy())
        else:
            raise ValueError("unknown multiscale composition")
    return {
        "interpretation": interpretation,
        "composition": composition,
        "dynamics": dynamics,
        "targets": targets,
        "increments": increments,
        "clamp_stages": clamp_stages,
        "next": result,
    }


# Cube sphere ---------------------------------------------------------------

Vec3 = tuple[float, float, float]
SPHERE_FRAMES: tuple[tuple[Vec3, Vec3, Vec3], ...] = (
    ((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)),
    ((0.0, 1.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 0.0, 1.0)),
    ((-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, 1.0)),
    ((0.0, -1.0, 0.0), (1.0, 0.0, 0.0), (0.0, 0.0, 1.0)),
    ((0.0, 0.0, 1.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0)),
    ((0.0, 0.0, -1.0), (1.0, 0.0, 0.0), (0.0, -1.0, 0.0)),
)


def _dot(left: Vec3, right: Vec3) -> float:
    return left[0] * right[0] + left[1] * right[1] + left[2] * right[2]


def _cross(left: Vec3, right: Vec3) -> Vec3:
    return (
        left[1] * right[2] - left[2] * right[1],
        left[2] * right[0] - left[0] * right[2],
        left[0] * right[1] - left[1] * right[0],
    )


def _normalize(value: Vec3) -> Vec3:
    length = math.sqrt(_dot(value, value))
    return value[0] / length, value[1] / length, value[2] / length


def _face_direction(face: int, u: float, v: float) -> Vec3:
    normal, axis_u, axis_v = SPHERE_FRAMES[face]
    tu = math.tan(u * math.pi / 4.0)
    tv = math.tan(v * math.pi / 4.0)
    return _normalize(
        (
            normal[0] + tu * axis_u[0] + tv * axis_v[0],
            normal[1] + tu * axis_u[1] + tv * axis_v[1],
            normal[2] + tu * axis_u[2] + tv * axis_v[2],
        )
    )


def sphere_direction(face: int, x: int, y: int, k: int) -> Vec3:
    u = 2.0 * (x + 0.5) / k - 1.0
    v = 2.0 * (y + 0.5) / k - 1.0
    return _face_direction(face, u, v)


def _spherical_triangle_area(a: Vec3, b: Vec3, c: Vec3) -> float:
    numerator = abs(_dot(a, _cross(b, c)))
    denominator = 1.0 + _dot(a, b) + _dot(b, c) + _dot(c, a)
    return 2.0 * math.atan2(numerator, denominator)


def sphere_cell_area(face: int, x: int, y: int, k: int, radius: float) -> float:
    u0, u1 = 2.0 * x / k - 1.0, 2.0 * (x + 1) / k - 1.0
    v0, v1 = 2.0 * y / k - 1.0, 2.0 * (y + 1) / k - 1.0
    a = _face_direction(face, u0, v0)
    b = _face_direction(face, u1, v0)
    c = _face_direction(face, u1, v1)
    d = _face_direction(face, u0, v1)
    return (
        radius
        * radius
        * (_spherical_triangle_area(a, b, c) + _spherical_triangle_area(a, c, d))
    )


def sphere_geometry(k: int) -> dict[str, Any]:
    if k <= 0 or k % 2:
        raise ValueError("sphere fixture requires positive even K")
    radius = k / 2.0
    directions: list[list[float]] = []
    areas: list[float] = []
    for face in range(6):
        for y in range(k):
            for x in range(k):
                directions.append(list(sphere_direction(face, x, y, k)))
                areas.append(sphere_cell_area(face, x, y, k, radius))
    return {
        "k": k,
        "radius": radius,
        "directions": directions,
        "areas": areas,
        "total_area": sum(areas),
        "analytic_area": 4.0 * math.pi * radius * radius,
    }


# Each row is (source face, destination face, destination quad corners). The source
# corners are always active-face [(1,1),(2,1),(2,2),(1,2)]. This is the original
# 24-side atlas table; diagonal gutters are intentionally absent/masked.
_LEGACY_SIDE_ROWS: tuple[tuple[int, int, tuple[tuple[int, int], ...]], ...] = (
    (0, 1, ((0, 1), (1, 1), (1, 2), (0, 2))),
    (0, 3, ((2, 1), (3, 1), (3, 2), (2, 2))),
    (0, 4, ((3, 1), (3, 2), (2, 2), (2, 1))),
    (0, 5, ((2, 2), (2, 1), (3, 1), (3, 2))),
    (1, 0, ((2, 1), (3, 1), (3, 2), (2, 2))),
    (1, 2, ((0, 1), (1, 1), (1, 2), (0, 2))),
    (1, 4, ((2, 3), (1, 3), (1, 2), (2, 2))),
    (1, 5, ((2, 1), (1, 1), (1, 0), (2, 0))),
    (2, 1, ((2, 1), (3, 1), (3, 2), (2, 2))),
    (2, 3, ((0, 1), (1, 1), (1, 2), (0, 2))),
    (2, 4, ((0, 2), (0, 1), (1, 1), (1, 2))),
    (2, 5, ((1, 1), (1, 2), (0, 2), (0, 1))),
    (3, 0, ((0, 1), (1, 1), (1, 2), (0, 2))),
    (3, 2, ((2, 1), (3, 1), (3, 2), (2, 2))),
    (3, 4, ((1, 0), (2, 0), (2, 1), (1, 1))),
    (3, 5, ((1, 2), (2, 2), (2, 3), (1, 3))),
    (4, 0, ((1, 3), (1, 2), (2, 2), (2, 3))),
    (4, 1, ((2, 3), (1, 3), (1, 2), (2, 2))),
    (4, 2, ((2, 2), (2, 3), (1, 3), (1, 2))),
    (4, 3, ((1, 2), (2, 2), (2, 3), (1, 3))),
    (5, 0, ((2, 0), (2, 1), (1, 1), (1, 0))),
    (5, 1, ((2, 1), (1, 1), (1, 0), (2, 0))),
    (5, 2, ((1, 1), (1, 0), (2, 0), (2, 1))),
    (5, 3, ((1, 0), (2, 0), (2, 1), (1, 1))),
)


def legacy_sphere_map(face: int, x: int, y: int, k: int) -> tuple[int, int, int] | None:
    if 0 <= x < k and 0 <= y < k:
        return face, x, y
    if (x < 0 or x >= k) and (y < 0 or y >= k):
        return None
    point = (1.0 + (x + 0.5) / k, 1.0 + (y + 0.5) / k)
    if x < 0:
        bounds = (0, 1, 1, 2)
    elif x >= k:
        bounds = (2, 3, 1, 2)
    elif y < 0:
        bounds = (1, 2, 0, 1)
    else:
        bounds = (1, 2, 2, 3)
    for source_face, destination_face, corners in _LEGACY_SIDE_ROWS:
        xs = [corner[0] for corner in corners]
        ys = [corner[1] for corner in corners]
        if destination_face != face or (min(xs), max(xs), min(ys), max(ys)) != bounds:
            continue
        d0, d1, _, d3 = corners
        e1 = (d1[0] - d0[0], d1[1] - d0[1])
        e2 = (d3[0] - d0[0], d3[1] - d0[1])
        relative = (point[0] - d0[0], point[1] - d0[1])
        determinant = e1[0] * e2[1] - e1[1] * e2[0]
        alpha = (relative[0] * e2[1] - relative[1] * e2[0]) / determinant
        beta = (e1[0] * relative[1] - e1[1] * relative[0]) / determinant
        source_x = min(k - 1, max(0, math.floor(alpha * k)))
        source_y = min(k - 1, max(0, math.floor(beta * k)))
        return source_face, source_x, source_y
    raise AssertionError((face, x, y, bounds))


def planar_to_geodesic(planar_radius: float, sphere_radius: float) -> float:
    argument = 1.0 - planar_radius * planar_radius / (
        2.0 * sphere_radius * sphere_radius
    )
    if not -1.0 <= argument <= 1.0:
        raise ValueError("planar chord radius is invalid for this sphere")
    return sphere_radius * math.acos(clamp(argument, -1.0, 1.0))


def spherical_cap_area(radius: float, sphere_radius: float) -> float:
    return (
        2.0
        * math.pi
        * sphere_radius
        * sphere_radius
        * (1.0 - math.cos(radius / sphere_radius))
    )


def sphere_neighborhoods(
    field: Sequence[float], k: int, ra_planar: float, model: str
) -> dict[str, Any]:
    if len(field) != 6 * k * k or model not in ("corrected", "legacy"):
        raise ValueError("invalid sphere field or model")
    geometry = sphere_geometry(k)
    radius = k / 2.0
    ri_geo = planar_to_geodesic(ra_planar / 3.0, radius)
    ra_geo = planar_to_geodesic(ra_planar, radius)
    directions = [tuple(value) for value in geometry["directions"]]
    areas = list(geometry["areas"])
    m_values: list[float] = []
    n_values: list[float] = []
    disk_denominators: list[float] = []
    ring_denominators: list[float] = []
    visited_disk_sums: list[float] = []
    visited_ring_sums: list[float] = []
    search = math.ceil(2.0 * ra_planar)
    legacy_disk = spherical_cap_area(ri_geo, radius)
    legacy_ring = spherical_cap_area(ra_geo, radius) - legacy_disk
    for face in range(6):
        for y in range(k):
            for x in range(k):
                center_index = face * k * k + y * k + x
                center = directions[center_index]
                disk_numerator = ring_numerator = 0.0
                disk_denominator = ring_denominator = 0.0
                if model == "corrected":
                    candidates = (
                        (candidate_face, candidate_x, candidate_y)
                        for candidate_face in range(6)
                        for candidate_y in range(k)
                        for candidate_x in range(k)
                    )
                else:
                    mapped: list[tuple[int, int, int]] = []
                    for dy in range(-search, search + 1):
                        for dx in range(-search, search + 1):
                            candidate = legacy_sphere_map(face, x + dx, y + dy, k)
                            if candidate is not None:
                                mapped.append(candidate)
                    candidates = iter(mapped)
                for candidate_face, candidate_x, candidate_y in candidates:
                    candidate_index = (
                        candidate_face * k * k + candidate_y * k + candidate_x
                    )
                    candidate = directions[candidate_index]
                    distance = radius * math.acos(
                        clamp(_dot(center, candidate), -1.0, 1.0)
                    )
                    disk_weight = 1.0 - transition_l(distance, ri_geo, 1.0)
                    ring_weight = transition_l(distance, ri_geo, 1.0) * (
                        1.0 - transition_l(distance, ra_geo, 1.0)
                    )
                    area = areas[candidate_index]
                    weighted_value = field[candidate_index] * area
                    disk_numerator += weighted_value * disk_weight
                    ring_numerator += weighted_value * ring_weight
                    disk_denominator += area * disk_weight
                    ring_denominator += area * ring_weight
                visited_disk, visited_ring = disk_denominator, ring_denominator
                if model == "legacy":
                    disk_denominator, ring_denominator = legacy_disk, legacy_ring
                m_values.append(disk_numerator / disk_denominator)
                n_values.append(ring_numerator / ring_denominator)
                disk_denominators.append(disk_denominator)
                ring_denominators.append(ring_denominator)
                visited_disk_sums.append(visited_disk)
                visited_ring_sums.append(visited_ring)
    return {
        "model": model,
        "ri_geodesic": ri_geo,
        "ra_geodesic": ra_geo,
        "search_bound": None if model == "corrected" else search,
        "disk_denominators": disk_denominators,
        "ring_denominators": ring_denominators,
        "visited_disk_sums": visited_disk_sums,
        "visited_ring_sums": visited_ring_sums,
        "m": m_values,
        "n": n_values,
    }


def sphere_step(
    field: Sequence[float], k: int, ra_planar: float, rule: Rule, model: str
) -> dict[str, Any]:
    result = sphere_neighborhoods(field, k, ra_planar, model)
    target = [
        rule_target(n, m, rule) for n, m in zip(result["n"], result["m"], strict=True)
    ]
    result.update(
        {
            "s": target,
            "next_discrete": [clamp(value) for value in target],
            "next_smooth": [
                clamp(value + 0.1 * (2.0 * desired - 1.0))
                for value, desired in zip(field, target, strict=True)
            ],
        }
    )
    return result


def sphere_overlays(k: int, seed: int, draws: int = 1000) -> dict[str, Any]:
    if k <= 0 or k % 2:
        raise ValueError("sphere initializer requires positive even K")
    rng = ChaCha12(seed)
    radius = k / 2.0
    directions = [
        sphere_direction(face, x, y, k)
        for face in range(6)
        for y in range(k)
        for x in range(k)
    ]
    field = [0.0] * (6 * k * k)
    draw_records: list[dict[str, Any]] = []
    for _ in range(draws):
        face = rng.integer(0, 6)
        x = rng.integer(0, k)
        y = rng.integer(0, k)
        paint_radius = rng.integer(2, 8)
        value = 1.0 if rng.integer(0, 2) else 0.0
        if len(draw_records) < 8:
            draw_records.append(
                {"face": face, "x": x, "y": y, "radius": paint_radius, "value": value}
            )
        center = directions[face * k * k + y * k + x]
        for index, candidate in enumerate(directions):
            distance = radius * math.acos(clamp(_dot(center, candidate), -1.0, 1.0))
            if distance < paint_radius:
                field[index] = value
    return {
        "k": k,
        "seed": seed,
        "draw_count": draws,
        "draw_order": ["face", "x", "y", "radius_integer_2_through_7", "binary_value"],
        "first_draws": draw_records,
        "field": field,
    }


# Delayed-time --------------------------------------------------------------


class DelayedStencilEntry(TypedDict):
    dx: int
    dy: int
    distance: float
    delay: int
    disk: float
    ring: float


class DelayedStencil(TypedDict):
    ra: float
    ri: float
    depth: int
    search: int
    disk_sum: float
    ring_sum: float
    entries: list[DelayedStencilEntry]


def delayed_stencil(ra: float, depth: int = 16) -> DelayedStencil:
    if ra <= 0.0 or depth <= 0:
        raise ValueError("invalid delayed-time geometry")
    ri = ra / 3.0
    search = math.ceil(ra + 0.5)
    entries: list[DelayedStencilEntry] = []
    disk_sum = ring_sum = 0.0
    for dy in range(-search, search + 1):
        for dx in range(-search, search + 1):
            distance = math.sqrt(dx * dx + dy * dy)
            disk = 1.0 - transition_l(distance, ri, 1.0)
            ring = transition_l(distance, ri, 1.0) * (
                1.0 - transition_l(distance, ra, 1.0)
            )
            delay = math.floor(distance + 0.5)
            entries.append(
                {
                    "dx": dx,
                    "dy": dy,
                    "distance": distance,
                    "delay": delay,
                    "disk": disk,
                    "ring": ring,
                }
            )
            disk_sum += disk
            ring_sum += ring
    if disk_sum == 0.0 or ring_sum == 0.0:
        raise ValueError("delayed stencil has zero normalization")
    return {
        "ra": ra,
        "ri": ri,
        "depth": depth,
        "search": search,
        "disk_sum": disk_sum,
        "ring_sum": ring_sum,
        "entries": entries,
    }


def delayed_step(
    history: Sequence[Sequence[float]],
    shape: Sequence[int],
    head: int,
    ra: float,
    rule: Rule,
) -> dict[str, Any]:
    if (
        len(shape) != 2
        or not history
        or any(len(layer) != product(shape) for layer in history)
    ):
        raise ValueError("delayed time requires a rectangular 2D history")
    depth = len(history)
    if not 0 <= head < depth:
        raise ValueError("head outside history")
    latest = (head - 1) % depth
    stencil = delayed_stencil(ra, depth)
    m_values: list[float] = []
    n_values: list[float] = []
    for x, y in iter_coords(shape):
        disk_numerator = ring_numerator = 0.0
        for entry in stencil["entries"]:
            layer = (latest - entry["delay"]) % depth
            source = ((x - entry["dx"]) % shape[0], (y - entry["dy"]) % shape[1])
            value = history[layer][flat_index(source, shape)]
            disk_numerator += value * entry["disk"]
            ring_numerator += value * entry["ring"]
        m_values.append(disk_numerator / stencil["disk_sum"])
        n_values.append(ring_numerator / stencil["ring_sum"])
    target = [rule_target(n, m, rule) for n, m in zip(n_values, m_values, strict=True)]
    discrete = [clamp(value) for value in target]
    smooth = [
        clamp(base + 0.1 * (2.0 * desired - 1.0))
        for base, desired in zip(history[latest], target, strict=True)
    ]
    return {
        "head": head,
        "latest": latest,
        "next_head": (head + 1) % depth,
        "stencil": stencil,
        "m": m_values,
        "n": n_values,
        "s": target,
        "next_discrete": discrete,
        "next_smooth": smooth,
    }


def delayed_boxes(
    shape: Sequence[int], seed: int, depth: int = 16, draws: int = 1000
) -> dict[str, Any]:
    if len(shape) != 2 or any(extent <= 0 for extent in shape) or depth <= 0:
        raise ValueError("invalid delayed box shape/depth")
    rng = ChaCha12(seed)
    field = [0.0] * product(shape)
    first_boxes: list[dict[str, Any]] = []
    for _ in range(draws):
        x = rng.integer(0, shape[0])
        y = rng.integer(0, shape[1])
        width = rng.integer(10, 20)
        height = rng.integer(10, 20)
        value = 1.0 if rng.integer(0, 2) else 0.0
        if len(first_boxes) < 8:
            first_boxes.append(
                {"x": x, "y": y, "width": width, "height": height, "value": value}
            )
        for dy in range(height):
            for dx in range(width):
                destination = ((x + dx) % shape[0], (y + dy) % shape[1])
                field[flat_index(destination, shape)] = value
    return {
        "shape": list(shape),
        "seed": seed,
        "depth": depth,
        "box_count": draws,
        "boundary": "half-open [x,x+width) x [y,y+height), periodic per sample",
        "draw_order": [
            "x",
            "y",
            "width_integer_10_through_19",
            "height_integer_10_through_19",
            "binary_value",
        ],
        "first_boxes": first_boxes,
        "frame": field,
        "history": [field.copy() for _ in range(depth)],
    }