diff --git a/tools/neural_sim/README.md b/tools/neural_sim/README.md
new file mode 100644
index 0000000..c294e2d
--- /dev/null
+++ b/tools/neural_sim/README.md
@@ -0,0 +1,184 @@
+# neural_sim — FPGA-Neural V2 golden functional reference
+
+    Python simulator = golden functional reference
+    RTL / P&R        = hardware implementation
+
+This package is **not** a cycle-accurate FPGA simulator. It does not
+model clock cycles, SDRAM controller timing, or SPI transaction
+timing. Given weights, activations, and a network topology, it
+computes the mathematically correct result that the real V2 hardware
+(`hardware/v2/rtl/neural_processor.v`) must reproduce bit-for-bit. Its
+job is to be the reference everything else — RTL simulation, and
+eventually real hardware — is checked against.
+
+## 1. Numeric model
+
+INT8 in / INT8 weight / INT32 accumulate / INT8 out, matching
+`neural_processor.v` exactly (re-derived by reading that file, not
+assumed — see `numerics.py`'s own module docstring for the full,
+line-by-line derivation):
+
+| Stage | Width | Overflow behaviour |
+|---|---|---|
+| Multiply (INT8 × INT8) | 16-bit signed product | Cannot overflow for INT8 operands |
+| Per-tile adder tree (P_IN=8 products) | 32-bit signed | Cannot overflow for P_IN=8 |
+| Cross-tile accumulator | 32-bit signed | **Wraparound** (true two's-complement, matches a Verilog `reg signed [31:0]`'s silent overflow) |
+| Bias add | 32-bit signed | **Wraparound** |
+| Final activation/output | 8-bit signed | **Saturating** (the only saturating stage) |
+
+`numerics.wrap_acc()` implements genuine 32-bit wraparound (not
+Python's arbitrary-precision integers hiding the boundary) —
+`tests/test_numerics.py::test_wrap_acc_true_32bit_wraparound` proves
+`2**31` wraps to `-(2**31)`, exactly like the RTL.
+
+This package reuses `tools/validation/mac_oracle.py`'s own
+independently-derived two's-complement primitives (`to_signed`,
+`mac8_tree`) rather than duplicating them — that file already is this
+project's own hand-verified oracle for the identical wraparound-add
+semantics used by `rtl/mac_unit.v`/`rtl/mac8.v`, which
+`neural_processor.v`'s own header states is the SAME accumulation
+lineage.
+
+## 2. Neuron equation
+
+```
+y = activate(bias + sum(x[i] * w[i] for i in 0..N_INPUTS-1))
+```
+
+computed tile-by-tile in groups of `P_IN` (default 8, matching the
+frozen hardware reference), with the cross-tile accumulator carrying
+state (with wraparound) between tiles — exactly how
+`neural_processor.v`'s own pipeline works (one tile enters per cycle,
+`acc_reg` only clears at job load, activation is applied once at the
+end after `tile_last`).
+
+Two independent implementations are provided and cross-tested
+(`tests/test_neuron.py`):
+- `neuron.neuron_scalar` — plain Python, easiest to audit line-by-line
+  against the RTL.
+- `neuron.neuron_vectorized` — NumPy-based (int64 intra-tile dot
+  products, since that stage can never overflow; explicit Python
+  wraparound arithmetic across tiles, so NumPy's own dtype-wraparound
+  behaviour is never silently relied on).
+
+## 3. Layer equation
+
+```
+output[n] = activate(bias[n] + sum(input[i] * weight[n][i] for i in 0..N_INPUTS-1))
+```
+
+for `n` in `0..N_NEURONS-1` — `layer.FCLayer`, matching how a real job
+batch is submitted (one `WRITE_JOB` per neuron, same `x_base`
+activation tile, each with its own `w_base`/`result_addr`).
+
+## 4. Network model (deliberately small scope)
+
+`network.Network` is an ordered chain of `FCLayer`s. This is a real,
+deliberate scope limit, not a hardware limit: the real
+`dependency_manager.v` schedules an arbitrary DAG of neuron jobs via
+`producer_ids`/`required` fields, not just linear layer chains. A
+linear chain is what this first phase implements and verifies; see
+"Optional future extension" below for what's deferred.
+
+## 5. Memory model
+
+`memory.MemoryModel` reproduces the real, official V2 unified-SDRAM
+memory map's logical byte contents and addresses (weights@0x010000,
+activations@0x200000, results@0x300000, all within the single 8MB
+0x000000–0x7FFFFF space) — **not** SDRAM cycle timing (see
+`hardware/v2/nms/rtl/sdram_controller.v` / `sdram_model.v` for that).
+`write_weight`/`read_weight`/`write_activation`/`read_activation`/
+`write_result`/`read_result` are bounds-checked against each region's
+own real base address.
+
+## 6. Quantization / activation behaviour
+
+Exactly two activation encodings exist in `neural_processor.v`:
+
+- **`relu`** (the RTL's own `default` case — i.e. any activation code
+  other than exactly `ACT_NONE` also produces ReLU): `acc <= 0 -> 0`;
+  `0 < acc <= 127 -> acc` exactly; `acc > 127 -> 127`. **This is the
+  only activation the real, currently-exposed V2 SPI job protocol
+  applies** (`spi_host_bridge.v`'s `WRITE_JOB` opcode has no
+  activation-selection field) — matches the V2 LaTeX datasheet's own
+  "Activation: ReLU + INT8 saturate, fixed".
+- **`none`**: a two-sided saturating clamp to the full signed INT8
+  range `[-128, 127]` — implemented in the module and modeled here for
+  completeness, but not reachable via the real, currently-exposed
+  protocol.
+
+No other quantization/scaling/shift stage exists in the real RTL, and
+none is invented here.
+
+## 7. FPGA correspondence
+
+| This simulator | Real hardware |
+|---|---|
+| `numerics.py` | `hardware/v2/rtl/neural_processor.v` (bit-exact, line-referenced) |
+| `memory.py` | `hardware/v2/nms/rtl/sdram_controller.v`'s logical address space (not its timing) |
+| `network.py` (linear chains only) | `hardware/v2/rtl/dependency_manager.v` (arbitrary DAG — a superset, not yet modeled here) |
+| — (not modeled) | `spi_host_bridge.v` transaction timing, SDRAM refresh/burst timing, tile-scheduling latency |
+
+## 8. CLI usage
+
+```
+python -m tools.neural_sim random-network --n-inputs 8 --n-neurons 4 --seed 1
+python -m tools.neural_sim run --example 8to4 --seed 7
+python -m tools.neural_sim vectors --gen d_stress --out vec.json
+python -m tools.neural_sim vectors --all --out all_vectors.json
+python -m tools.neural_sim compare --expected vec.json --actual fpga_results.json
+```
+
+Built-in examples (`--example`): `8to1`, `8to4`, `8to16`, `8_8_1` (an
+8→8→1 two-layer network — the hidden width is 8, not some other
+number, specifically because every layer boundary must stay a
+multiple of `P_IN=8`, the same tiling constraint the real hardware
+has).
+
+## 9. Test-vector generation and export
+
+`vectors.py` provides six deterministic generators (`gen_simple_positive`,
+`gen_signed_mix`, `gen_extremes`, `gen_zero`, `gen_random`, `gen_dstress`
+— the last reproducing the existing RTL benchmark's own 256-neuron ×
+128-input dimensions). Every vector's `expected` field is computed by
+this package's own golden model, not hand-typed, and stores `inputs`,
+`weights`, `biases`, `activation`, `p_in`, `numeric_format`, and `seed`
+(when applicable) — everything a future RTL testbench needs to load a
+vector and check its own result, without re-deriving anything by hand:
+
+```
+Python (vectors.py) -> JSON golden vectors -> (future) Verilog testbench -> FPGA
+```
+
+## 10. FPGA-result comparison
+
+`compare.compare_results(expected, actual)` returns a `ComparisonReport`
+with `exact_match`, `num_mismatches`, `first_mismatch_index`,
+`first_mismatch_expected`/`_actual`, and `max_abs_diff`. **0 mismatches
+is the only passing criterion** — this package never hides a real
+numeric difference behind a tolerance. `compare.load_fpga_results(path)`
+accepts either a JSON list of ints or a plain whitespace/line-separated
+results file (a common shape for an RTL testbench's own dump).
+
+## Running the tests
+
+```
+python3 -m pytest tools/neural_sim/tests/ -q
+```
+
+96 tests, all passing as of this writing: signed-arithmetic edge
+cases (including a direct 32-bit wraparound check), neuron
+scalar-vs-vectorized cross-checks (parametrized across tile counts,
+activations, and seeds), layer tests, memory bounds/adjacency tests,
+golden-vector determinism and JSON round-trip tests, and comparison-
+utility tests.
+
+## Optional future extension (not implemented — by design)
+
+This first phase deliberately stops at "the mathematical golden model
+is unquestionably correct." A future phase could add, without
+changing anything above: a cycle-accurate scheduler model,
+`N_PROCESSORS` ∈ {1,2,4} scheduling, tile scheduling, SDRAM traffic
+estimation, SPI transaction modeling, and latency prediction. None of
+that exists yet, and this package does not claim to be
+"cycle-accurate" anywhere.
diff --git a/tools/neural_sim/__init__.py b/tools/neural_sim/__init__.py
new file mode 100644
index 0000000..f032bcf
--- /dev/null
+++ b/tools/neural_sim/__init__.py
@@ -0,0 +1,19 @@
+"""
+neural_sim -- the Python GOLDEN FUNCTIONAL REFERENCE for FPGA-Neural V2.
+
+Given weights, activations, and a network topology, this package
+computes the mathematically correct result that the real V2 hardware
+(hardware/v2/rtl/neural_processor.v, unmodified since before the
+single-SDRAM freeze) must reproduce bit-for-bit.
+
+This is NOT a cycle-accurate FPGA simulator: it models the numeric
+result only, not clock cycles, memory-controller timing, or SPI
+transaction timing. See README.md for the full scope statement and
+tools/neural_sim/network.py / tools/neural_sim/memory.py for what IS
+and is NOT modeled.
+
+    Python simulator = golden functional reference
+    RTL / P&R        = hardware implementation
+"""
+
+__version__ = "0.1.0"
diff --git a/tools/neural_sim/__main__.py b/tools/neural_sim/__main__.py
new file mode 100644
index 0000000..1b0b754
--- /dev/null
+++ b/tools/neural_sim/__main__.py
@@ -0,0 +1,122 @@
+"""
+CLI for the neural_sim golden simulator.
+
+    python -m tools.neural_sim random-network --n-inputs 8 --n-neurons 4 --seed 1
+    python -m tools.neural_sim run --example 8to4 --seed 7
+    python -m tools.neural_sim vectors --gen d_stress --out vec.json
+    python -m tools.neural_sim vectors --all --out all_vectors.json
+    python -m tools.neural_sim compare --expected golden.json --actual fpga_results.json
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import sys
+
+import numpy as np
+
+from .layer import FCLayer
+from .network import Network
+from . import examples as ex
+from . import vectors as vec
+from .compare import compare_results, load_fpga_results
+
+
+def cmd_random_network(args):
+    rng = np.random.default_rng(args.seed)
+    weights = rng.integers(-args.magnitude, args.magnitude + 1,
+                            size=(args.n_neurons, args.n_inputs), dtype=np.int64)
+    layer = FCLayer(weights, activation=args.activation)
+    inputs = rng.integers(-128, 128, size=args.n_inputs, dtype=np.int64)
+    outputs = layer.forward(inputs)
+    print(f"n_inputs={layer.n_inputs} n_neurons={layer.n_neurons} activation={layer.activation}")
+    print(f"inputs:  {inputs.tolist()}")
+    print(f"outputs: {outputs.tolist()}")
+    if args.save_weights:
+        with open(args.save_weights, "w") as f:
+            json.dump({"weights": weights.tolist(), "inputs": inputs.tolist(),
+                       "outputs": outputs.tolist()}, f, indent=2)
+        print(f"saved to {args.save_weights}")
+
+
+def cmd_run(args):
+    if args.example not in ex.ALL_EXAMPLES:
+        print(f"unknown example {args.example!r}; choices: {list(ex.ALL_EXAMPLES)}", file=sys.stderr)
+        sys.exit(1)
+    net = ex.ALL_EXAMPLES[args.example]()
+    rng = np.random.default_rng(args.seed)
+    inputs = rng.integers(-128, 128, size=net.n_inputs, dtype=np.int64)
+    outputs = net.forward(inputs)
+    print(f"example={args.example} n_inputs={net.n_inputs} n_outputs={net.n_outputs}")
+    print(f"inputs:  {inputs.tolist()}")
+    print(f"outputs: {outputs.tolist()}")
+
+
+def cmd_vectors(args):
+    if args.all:
+        vec.export_all_json(args.out)
+        print(f"exported all {len(vec.ALL_GENERATORS)} named vectors to {args.out}")
+        return
+    if args.gen not in vec.ALL_GENERATORS:
+        print(f"unknown generator {args.gen!r}; choices: {list(vec.ALL_GENERATORS)}", file=sys.stderr)
+        sys.exit(1)
+    v = vec.ALL_GENERATORS[args.gen]()
+    vec.export_json(v, args.out)
+    print(f"generated {v.name!r}: n_inputs={v.n_inputs} n_neurons={v.n_neurons} -> {args.out}")
+
+
+def cmd_compare(args):
+    expected = vec.load_json(args.expected).expected if args.expected.endswith(".json") and _is_vector_file(args.expected) \
+        else load_fpga_results(args.expected)
+    actual = load_fpga_results(args.actual)
+    report = compare_results(expected, actual)
+    print(report.summary())
+    sys.exit(0 if report.exact_match else 1)
+
+
+def _is_vector_file(path: str) -> bool:
+    try:
+        with open(path) as f:
+            d = json.load(f)
+        return isinstance(d, dict) and "expected" in d and "weights" in d
+    except Exception:
+        return False
+
+
+def main():
+    parser = argparse.ArgumentParser(prog="python -m tools.neural_sim",
+                                      description="FPGA-Neural V2 golden functional reference simulator")
+    sub = parser.add_subparsers(dest="command", required=True)
+
+    p_rand = sub.add_parser("random-network", help="generate and run a random network")
+    p_rand.add_argument("--n-inputs", type=int, default=8)
+    p_rand.add_argument("--n-neurons", type=int, default=4)
+    p_rand.add_argument("--seed", type=int, default=0)
+    p_rand.add_argument("--magnitude", type=int, default=20, help="max abs weight magnitude")
+    p_rand.add_argument("--activation", choices=["relu", "none"], default="relu")
+    p_rand.add_argument("--save-weights", metavar="PATH", default=None)
+    p_rand.set_defaults(func=cmd_random_network)
+
+    p_run = sub.add_parser("run", help="run one of the built-in example networks")
+    p_run.add_argument("--example", choices=list(ex.ALL_EXAMPLES), default="8to4")
+    p_run.add_argument("--seed", type=int, default=0, help="seed for the random input vector")
+    p_run.set_defaults(func=cmd_run)
+
+    p_vec = sub.add_parser("vectors", help="generate golden test vectors")
+    p_vec.add_argument("--gen", choices=list(vec.ALL_GENERATORS), default=None)
+    p_vec.add_argument("--all", action="store_true", help="export every named generator")
+    p_vec.add_argument("--out", required=True, metavar="PATH")
+    p_vec.set_defaults(func=cmd_vectors)
+
+    p_cmp = sub.add_parser("compare", help="compare FPGA results against Python golden results")
+    p_cmp.add_argument("--expected", required=True, metavar="PATH",
+                        help="a vectors-format JSON file (uses its 'expected' field) or a plain results file")
+    p_cmp.add_argument("--actual", required=True, metavar="PATH", help="FPGA-generated results file")
+    p_cmp.set_defaults(func=cmd_compare)
+
+    args = parser.parse_args()
+    args.func(args)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/neural_sim/compare.py b/tools/neural_sim/compare.py
new file mode 100644
index 0000000..cf60c9b
--- /dev/null
+++ b/tools/neural_sim/compare.py
@@ -0,0 +1,82 @@
+"""
+FPGA-vs-Python comparison utility (B10). Loads FPGA-generated results
+and compares them against this package's own Python golden results.
+The primary pass/fail criterion is 0 mismatches -- exact bit-exact
+match, never a tolerance-based "close enough" comparison (per this
+project's own explicit "for the true bit-exact model, do not hide
+numerical differences behind tolerances" instruction).
+"""
+from __future__ import annotations
+
+import json
+from dataclasses import dataclass
+from typing import List, Optional
+
+
+@dataclass
+class ComparisonReport:
+    exact_match: bool
+    total: int
+    num_mismatches: int
+    first_mismatch_index: Optional[int]
+    first_mismatch_expected: Optional[int]
+    first_mismatch_actual: Optional[int]
+    max_abs_diff: int
+
+    def summary(self) -> str:
+        if self.exact_match:
+            return f"EXACT MATCH: {self.total}/{self.total} outputs bit-exact, 0 mismatches"
+        return (
+            f"MISMATCH: {self.num_mismatches}/{self.total} outputs differ "
+            f"(first at index {self.first_mismatch_index}: "
+            f"expected={self.first_mismatch_expected} actual={self.first_mismatch_actual}, "
+            f"max_abs_diff={self.max_abs_diff})"
+        )
+
+
+def compare_results(expected: List[int], actual: List[int]) -> ComparisonReport:
+    if len(expected) != len(actual):
+        raise ValueError(
+            f"length mismatch: expected has {len(expected)} outputs, "
+            f"actual has {len(actual)} -- cannot compare index-by-index"
+        )
+
+    total = len(expected)
+    num_mismatches = 0
+    first_idx = None
+    first_exp = None
+    first_act = None
+    max_abs_diff = 0
+
+    for i, (e, a) in enumerate(zip(expected, actual)):
+        if e != a:
+            num_mismatches += 1
+            if first_idx is None:
+                first_idx, first_exp, first_act = i, e, a
+            max_abs_diff = max(max_abs_diff, abs(e - a))
+
+    return ComparisonReport(
+        exact_match=(num_mismatches == 0),
+        total=total,
+        num_mismatches=num_mismatches,
+        first_mismatch_index=first_idx,
+        first_mismatch_expected=first_exp,
+        first_mismatch_actual=first_act,
+        max_abs_diff=max_abs_diff,
+    )
+
+
+def load_fpga_results(path: str) -> List[int]:
+    """Loads FPGA-generated results from either a JSON list of ints, or
+    a plain text file with one (whitespace-separated) integer per
+    line -- a common shape for an RTL testbench's own $display/
+    $writememh-style dump."""
+    with open(path) as f:
+        text = f.read()
+    try:
+        data = json.loads(text)
+        if isinstance(data, dict) and "results" in data:
+            data = data["results"]
+        return [int(v) for v in data]
+    except json.JSONDecodeError:
+        return [int(tok) for tok in text.split()]
diff --git a/tools/neural_sim/examples.py b/tools/neural_sim/examples.py
new file mode 100644
index 0000000..b65ea87
--- /dev/null
+++ b/tools/neural_sim/examples.py
@@ -0,0 +1,60 @@
+"""
+Small example networks for experimentation (B12). Every example uses
+operations the real V2 accelerator actually implements (INT8 in/INT8
+weight/INT32-wraparound-accumulate/ReLU-saturate-out, P_IN=8 tiles) --
+none invent unsupported functionality. The 2-layer example is
+FPGA-compatible at the computational level via the real dependency-
+graph mechanism (each layer-2 neuron's producer_ids/required fields
+would gate it on all 4 layer-1 neurons completing first), but this
+simulator does not yet model the SPI job-construction/address-wiring
+needed to actually run it end-to-end on hardware (see network.py's
+own scope note).
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from .layer import FCLayer
+from .network import Network
+
+
+def example_8_to_1() -> Network:
+    """8 inputs -> 1 neuron, ReLU. Hand-picked, easy-to-verify weights."""
+    weights = [[1, 2, 3, 4, -1, -2, -3, -4]]
+    return Network([FCLayer(weights, activation="relu")])
+
+
+def example_8_to_4() -> Network:
+    """8 inputs -> 4 neurons, ReLU. Fixed-seed weights for reproducibility."""
+    rng = np.random.default_rng(1)
+    weights = rng.integers(-20, 21, size=(4, 8), dtype=np.int64)
+    return Network([FCLayer(weights, activation="relu")])
+
+
+def example_8_to_16() -> Network:
+    """8 inputs -> 16 neurons, ReLU. Fixed-seed weights."""
+    rng = np.random.default_rng(2)
+    weights = rng.integers(-20, 21, size=(16, 8), dtype=np.int64)
+    return Network([FCLayer(weights, activation="relu")])
+
+
+def example_8_8_1() -> Network:
+    """8 -> 8 -> 1, both layers ReLU. Fixed-seed weights. The hidden
+    layer is 8 wide (not, say, 4) because every layer boundary must
+    stay a multiple of P_IN=8 -- the real hardware always tiles in
+    groups of 8, so a hidden width that doesn't divide evenly would not
+    be a layer this simulator's own FCLayer (or the real accelerator)
+    can actually tile. See module docstring for the FPGA-compatibility
+    scope note re: multi-layer chaining."""
+    rng = np.random.default_rng(3)
+    w1 = rng.integers(-15, 16, size=(8, 8), dtype=np.int64)
+    w2 = rng.integers(-15, 16, size=(1, 8), dtype=np.int64)
+    return Network([FCLayer(w1, activation="relu"), FCLayer(w2, activation="relu")])
+
+
+ALL_EXAMPLES = {
+    "8to1": example_8_to_1,
+    "8to4": example_8_to_4,
+    "8to16": example_8_to_16,
+    "8_8_1": example_8_8_1,
+}
diff --git a/tools/neural_sim/layer.py b/tools/neural_sim/layer.py
new file mode 100644
index 0000000..a26bdca
--- /dev/null
+++ b/tools/neural_sim/layer.py
@@ -0,0 +1,73 @@
+"""
+Fully-connected layer: N_INPUTS -> N_NEURONS, each neuron with its own
+weight vector (and, optionally, its own bias), all sharing one
+activation setting -- matching how a real V2 job batch is submitted
+(one WRITE_JOB per neuron, all against the same x_base activation
+tile, each with its own w_base/result_addr).
+
+output[n] = activation(bias[n] + sum(input[i]*weight[n][i] for i in
+0..N_INPUTS-1)), computed with the EXACT same tile-wise, wraparound-
+accumulate-then-saturate semantics as numerics.neuron_reference.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from .numerics import ACT_RELU, DEFAULT_P_IN, check_int8
+from .neuron import neuron_vectorized
+
+
+class FCLayer:
+    """
+    weights: array-like, shape (n_neurons, n_inputs), signed INT8.
+    biases:  array-like, shape (n_neurons,), signed INT8 (default: all
+             zero, matching the real V2 protocol's own lack of a bias
+             field -- see numerics.py's own module docstring).
+    activation: 'relu' (default, matches the real, currently-exposed
+             V2 system) or 'none'.
+    p_in:    tile width (default 8, matches the frozen P_IN=8 reference).
+    """
+
+    def __init__(self, weights, biases=None, activation: str = ACT_RELU,
+                 p_in: int = DEFAULT_P_IN):
+        self.weights = np.asarray(weights, dtype=np.int64)
+        if self.weights.ndim != 2:
+            raise ValueError("weights must be 2D: (n_neurons, n_inputs)")
+        self.n_neurons, self.n_inputs = self.weights.shape
+        if self.n_inputs % p_in != 0:
+            raise ValueError(f"n_inputs={self.n_inputs} must be a multiple of p_in={p_in}")
+        for row in self.weights.tolist():
+            for w in row:
+                check_int8(w, "weight")
+
+        if biases is None:
+            self.biases = np.zeros(self.n_neurons, dtype=np.int64)
+        else:
+            self.biases = np.asarray(biases, dtype=np.int64)
+            if self.biases.shape != (self.n_neurons,):
+                raise ValueError("biases must have shape (n_neurons,)")
+            for b in self.biases.tolist():
+                check_int8(b, "bias")
+
+        self.activation = activation
+        self.p_in = p_in
+
+    def forward(self, inputs) -> np.ndarray:
+        """inputs: array-like, shape (n_inputs,), signed INT8. Returns
+        an int64 numpy array of shape (n_neurons,) -- each element is
+        an exact signed INT8 value (kept as int64 purely for easy
+        downstream composition; every value is guaranteed in
+        [-128, 127])."""
+        inputs = np.asarray(inputs, dtype=np.int64)
+        if inputs.shape != (self.n_inputs,):
+            raise ValueError(f"inputs must have shape ({self.n_inputs},), got {inputs.shape}")
+        for v in inputs.tolist():
+            check_int8(v, "input")
+
+        outputs = np.empty(self.n_neurons, dtype=np.int64)
+        for n in range(self.n_neurons):
+            outputs[n] = neuron_vectorized(
+                inputs, self.weights[n], bias=int(self.biases[n]),
+                activation=self.activation, p_in=self.p_in,
+            )
+        return outputs
diff --git a/tools/neural_sim/memory.py b/tools/neural_sim/memory.py
new file mode 100644
index 0000000..8db6e6a
--- /dev/null
+++ b/tools/neural_sim/memory.py
@@ -0,0 +1,110 @@
+"""
+Logical memory model of the real, official V2 unified-SDRAM memory map
+(hardware/v2/docs/DatasheetLatex/chapters/02-architecture.tex's own
+"Official V2 memory map" table, cross-checked against
+hardware/v2/nms/sim/tb_sdram_boundary.v's own use of the same
+addresses):
+
+    Region       Base address   Notes
+    Weights      0x010000       1MB-aligned
+    Activations  0x200000       1MB-aligned
+    Results      0x300000       1MB-aligned
+
+all three non-overlapping within the single 8MB (0x000000-0x7FFFFF)
+SDRAM address space. This model reproduces the LOGICAL byte contents
+and addresses only -- it does NOT reproduce SDRAM cycle timing,
+refresh, or burst behaviour (see hardware/v2/nms/sim/sdram_model.v /
+sdram_controller.v for that; this is a plain flat byte array).
+
+Region sizes are inferred from adjacency (each region's own end is the
+next region's own base) -- the real hardware does not enforce region
+size limits in the datapath itself (base addresses are host-
+programmable per job), so this is this model's own deliberately
+conservative bounds-checking convention, matching the same assumption
+tb_sdram_boundary.v's own "weights-last(pre-act)"/"activations-
+last(pre-res)" checks use.
+"""
+from __future__ import annotations
+
+SDRAM_SIZE = 8 * 1024 * 1024  # 8MB, 0x000000-0x7FFFFF
+
+WEIGHTS_BASE = 0x010000
+ACTIVATIONS_BASE = 0x200000
+RESULTS_BASE = 0x300000
+
+WEIGHTS_END = ACTIVATIONS_BASE       # exclusive
+ACTIVATIONS_END = RESULTS_BASE       # exclusive
+RESULTS_END = SDRAM_SIZE             # exclusive
+
+
+def _to_unsigned8(v: int) -> int:
+    return v & 0xFF
+
+
+def _to_signed8(v: int) -> int:
+    v &= 0xFF
+    return v - 256 if v >= 128 else v
+
+
+class MemoryModel:
+    """A flat 8MB byte array standing in for the real unified SDRAM,
+    with bounds-checked, region-aware, signed-INT8 read/write helpers."""
+
+    def __init__(self):
+        self._mem = bytearray(SDRAM_SIZE)
+
+    # ---- raw byte access (any address in the full 8MB space) ----
+    def read_byte(self, addr: int) -> int:
+        if not (0 <= addr < SDRAM_SIZE):
+            raise IndexError(f"address 0x{addr:06x} out of range [0, 0x{SDRAM_SIZE:06x})")
+        return _to_signed8(self._mem[addr])
+
+    def write_byte(self, addr: int, value: int) -> None:
+        if not (0 <= addr < SDRAM_SIZE):
+            raise IndexError(f"address 0x{addr:06x} out of range [0, 0x{SDRAM_SIZE:06x})")
+        if not (-128 <= value <= 127):
+            raise ValueError(f"value {value} out of signed INT8 range [-128, 127]")
+        self._mem[addr] = _to_unsigned8(value)
+
+    # ---- region-aware helpers (B7's own required helper names) ----
+    def _region_check(self, base: int, end: int, offset: int, label: str) -> int:
+        addr = base + offset
+        if not (base <= addr < end):
+            raise IndexError(
+                f"{label} offset {offset} (address 0x{addr:06x}) falls outside "
+                f"its own region [0x{base:06x}, 0x{end:06x})"
+            )
+        return addr
+
+    def write_weight(self, offset: int, value: int) -> None:
+        self.write_byte(self._region_check(WEIGHTS_BASE, WEIGHTS_END, offset, "weight"), value)
+
+    def read_weight(self, offset: int) -> int:
+        return self.read_byte(self._region_check(WEIGHTS_BASE, WEIGHTS_END, offset, "weight"))
+
+    def write_activation(self, offset: int, value: int) -> None:
+        self.write_byte(self._region_check(ACTIVATIONS_BASE, ACTIVATIONS_END, offset, "activation"), value)
+
+    def read_activation(self, offset: int) -> int:
+        return self.read_byte(self._region_check(ACTIVATIONS_BASE, ACTIVATIONS_END, offset, "activation"))
+
+    def write_result(self, offset: int, value: int) -> None:
+        self.write_byte(self._region_check(RESULTS_BASE, RESULTS_END, offset, "result"), value)
+
+    def read_result(self, offset: int) -> int:
+        return self.read_byte(self._region_check(RESULTS_BASE, RESULTS_END, offset, "result"))
+
+    # ---- bulk convenience (not a hardware concept, pure host-side sugar) ----
+    def write_weights(self, offset: int, values) -> None:
+        for i, v in enumerate(values):
+            self.write_weight(offset + i, int(v))
+
+    def read_weights(self, offset: int, count: int):
+        return [self.read_weight(offset + i) for i in range(count)]
+
+    def write_activations(self, offset: int, values) -> None:
+        for i, v in enumerate(values):
+            self.write_activation(offset + i, int(v))
+
+    def read_activations(self, offset: int, count: int):
+        return [self.read_activation(offset + i) for i in range(count)]
diff --git a/tools/neural_sim/network.py b/tools/neural_sim/network.py
new file mode 100644
index 0000000..10da76d
--- /dev/null
+++ b/tools/neural_sim/network.py
@@ -0,0 +1,61 @@
+"""
+Network model: an ordered chain of FCLayer instances, each layer's
+INT8 output feeding the next layer's INT8 input.
+
+Scope, deliberately kept small (per this project's own explicit "do
+not invent unsupported FPGA functionality" instruction): this models a
+linear chain of fully-connected+activation layers. The real V2
+hardware's dependency_manager.v is actually more general -- it
+schedules an arbitrary DAG of neuron "jobs" via producer_ids/required
+fields, so a layer boundary is not a hardware limitation, only a
+simulator scope limit for this first phase. A linear chain is exactly
+what a linear chain of dependency-graph layers computes, so this is
+faithful for the topologies it supports; it does not yet model
+arbitrary-DAG job graphs, SPI job submission, or scheduling -- see
+README.md "Optional future extension" and PRE_PCB_VERIFICATION.md's
+own dependency-graph description for what the real hardware supports
+beyond what this simulator currently models.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from .layer import FCLayer
+
+
+class Network:
+    def __init__(self, layers: list[FCLayer]):
+        if not layers:
+            raise ValueError("a Network needs at least one layer")
+        for i in range(1, len(layers)):
+            if layers[i].n_inputs != layers[i - 1].n_neurons:
+                raise ValueError(
+                    f"layer {i}'s n_inputs={layers[i].n_inputs} does not match "
+                    f"layer {i-1}'s n_neurons={layers[i-1].n_neurons}"
+                )
+        self.layers = layers
+
+    @property
+    def n_inputs(self) -> int:
+        return self.layers[0].n_inputs
+
+    @property
+    def n_outputs(self) -> int:
+        return self.layers[-1].n_neurons
+
+    def forward(self, inputs) -> np.ndarray:
+        """Runs `inputs` through every layer in order, returning the
+        final layer's INT8 output vector. Also available as
+        `forward_all` if every intermediate tensor is needed."""
+        return self.forward_all(inputs)[-1]
+
+    def forward_all(self, inputs) -> list[np.ndarray]:
+        """Returns [layer0_output, layer1_output, ..., layerN_output]
+        -- every intermediate tensor, not just the final one (useful
+        for debugging / per-layer golden-vector generation)."""
+        x = np.asarray(inputs, dtype=np.int64)
+        outputs = []
+        for layer in self.layers:
+            x = layer.forward(x)
+            outputs.append(x)
+        return outputs
diff --git a/tools/neural_sim/neuron.py b/tools/neural_sim/neuron.py
new file mode 100644
index 0000000..eda0425
--- /dev/null
+++ b/tools/neural_sim/neuron.py
@@ -0,0 +1,65 @@
+"""
+Neuron model: y = activation(bias + sum(x[i]*w[i] for i in 0..N_INPUTS-1)),
+matching neural_processor.v (P_IN=8 baseline, parameterizable). Two
+independent implementations are provided -- a plain-Python scalar
+reference (numerics.neuron_reference) and a NumPy-vectorized one below
+-- and tests/test_neuron.py proves they always agree, bit-exact,
+across every test-vector category.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from .numerics import (
+    DEFAULT_ACC_WIDTH, DEFAULT_DATA_WIDTH, DEFAULT_P_IN,
+    ACT_RELU, accumulate_tile, add_bias, activate_and_saturate,
+    check_int8, neuron_reference, wrap_acc,
+)
+
+
+def neuron_scalar(inputs, weights, bias: int = 0, activation: str = ACT_RELU,
+                   p_in: int = DEFAULT_P_IN) -> int:
+    """Plain-Python scalar reference -- a thin, explicit re-export of
+    numerics.neuron_reference (kept as its own name so neuron.py has a
+    clearly-named "scalar" counterpart to neuron_vectorized below)."""
+    return neuron_reference(inputs, weights, bias=bias, activation=activation, p_in=p_in)
+
+
+def neuron_vectorized(inputs, weights, bias: int = 0, activation: str = ACT_RELU,
+                       p_in: int = DEFAULT_P_IN,
+                       acc_width: int = DEFAULT_ACC_WIDTH,
+                       data_width: int = DEFAULT_DATA_WIDTH) -> int:
+    """
+    NumPy-vectorized equivalent of neuron_scalar. Uses int64 numpy
+    accumulation WITHIN one tile only (exact -- a single P_IN=8 tile's
+    product-sum can never exceed a few hundred thousand in magnitude,
+    nowhere near int64's range, so no numpy overflow risk there), then
+    applies the SAME explicit 32-bit-wraparound Python arithmetic as the
+    scalar path across tiles -- deliberately NOT letting NumPy's own
+    int32 dtype silently wrap with platform-dependent behaviour, per
+    this project's own "do not allow implicit integer promotion to hide
+    overflow" requirement.
+    """
+    x = np.asarray(inputs, dtype=np.int64)
+    w = np.asarray(weights, dtype=np.int64)
+    if x.shape != w.shape:
+        raise ValueError("inputs and weights must have the same shape")
+    if x.size == 0 or x.size % p_in != 0:
+        raise ValueError(f"len(inputs)={x.size} must be a nonzero multiple of p_in={p_in}")
+    for v in x.tolist():
+        check_int8(v, "input")
+    for v in w.tolist():
+        check_int8(v, "weight")
+
+    n_tiles = x.size // p_in
+    x_tiles = x.reshape(n_tiles, p_in)
+    w_tiles = w.reshape(n_tiles, p_in)
+    # exact per-tile dot products (int64, no overflow risk)
+    tile_sums = np.einsum("ij,ij->i", x_tiles, w_tiles)
+
+    acc = 0
+    for tile_sum in tile_sums.tolist():
+        acc = accumulate_tile(acc, int(tile_sum), acc_width=acc_width)
+
+    final_acc = add_bias(acc, bias, acc_width=acc_width, data_width=data_width)
+    return activate_and_saturate(final_acc, activation=activation, data_width=data_width)
diff --git a/tools/neural_sim/numerics.py b/tools/neural_sim/numerics.py
new file mode 100644
index 0000000..3c1a9f0
--- /dev/null
+++ b/tools/neural_sim/numerics.py
@@ -0,0 +1,193 @@
+"""
+Bit-exact numeric semantics of hardware/v2/rtl/neural_processor.v.
+
+The RTL is authoritative. Every function here was derived by reading
+neural_processor.v directly (not assumed), specifically:
+
+  Stage 1 (line ~139): product_comb[gm] = x0[gm] * w0[gm]
+      INT8 x INT8 signed multiply. PROD_WIDTH = 2*DATA_WIDTH = 16 bits
+      is always sufficient (min product -128*127=-16256, max
+      -128*-128=16384, both fit in signed 16 bits) -- the multiply
+      itself can NEVER overflow for DATA_WIDTH=8. Each product is then
+      sign-extended to ACC_WIDTH=32 bits (line ~151).
+
+  Stages 2..(1+TREE_LEVELS) (line ~166-206): a balanced binary adder
+      tree reduces the P_IN products to one `tile_sum`, entirely in
+      ACC_WIDTH=32-bit signed arithmetic. For P_IN=8 this can never
+      overflow either (max magnitude 8*16384=131072 << 2^31).
+
+  Stage (2+TREE_LEVELS) (line ~212-229): `acc_reg <= acc_reg +
+      tile_sum` -- a RUNNING accumulator across ALL TILES of one job
+      (cleared only at NP_LOAD_JOB), using plain Verilog `+` on a
+      32-bit signed reg. This is WRAPAROUND (modulo 2^32) arithmetic,
+      NOT saturating -- Verilog silently wraps a fixed-width `+`.
+      Practically never triggered for realistic tile counts (a job
+      would need on the order of 2^17 tiles for the running sum to
+      approach the 32-bit signed range), but modeled as true wraparound
+      here anyway, per this project's own explicit "do not let overflow
+      hide behind Python's arbitrary precision" requirement, and
+      because an adversarial/stress test vector may deliberately probe
+      this boundary.
+
+  Stage (3+TREE_LEVELS) (line ~235-260): `final_acc_reg <= acc_reg +
+      bias_ext` -- bias (sign-extended from an INT8 job_bias field) is
+      ALSO added with 32-bit wraparound semantics. NOTE: the real,
+      currently-exposed V2 SPI job protocol (spi_host_bridge.v's own
+      WRITE_JOB opcode) has no bias field at all -- bias is a
+      neural_processor.v MODULE-LEVEL capability, not something the
+      real V2 host can currently set. This model defaults bias=0 to
+      match the real, currently-exposed system behaviour, while still
+      implementing nonzero bias faithfully for anyone driving
+      neural_processor.v directly.
+
+  Stage (4+TREE_LEVELS) (line ~262-288): the ONLY saturating stage.
+      Two activation encodings exist in the RTL: ACT_NONE (a two-sided
+      saturating clamp to the full signed INT8 range [-128, 127]) and
+      ACT_RELU (the Verilog `case` statement's `default` branch, so
+      ANY activation code other than exactly ACT_NONE=0 also produces
+      ReLU behaviour). The real, currently-exposed V2 SPI job protocol
+      has no activation-selection field either -- ReLU is the only
+      activation the real system currently applies (matches the V2
+      LaTeX datasheet's own "Activation: ReLU + INT8 saturate, fixed").
+
+Overflow/signedness summary (as explicitly requested):
+  - multiplication: exact, cannot overflow for INT8 operands
+  - per-tile adder tree: exact for P_IN<=8, WRAPAROUND semantics modeled
+  - cross-tile accumulator: WRAPAROUND (32-bit, two's complement)
+  - bias add: WRAPAROUND (32-bit, two's complement)
+  - final activation/output: SATURATING (to INT8, either two-sided for
+    ACT_NONE or ReLU-then-saturate for ACT_RELU)
+
+Reuses tools/validation/mac_oracle.py's own independently-derived
+two's-complement primitives (`to_signed`, `mac8_tree`) rather than
+duplicating them -- that file is this project's own pre-existing,
+already-hand-verified oracle for the identical wraparound-add
+semantics (rtl/mac_unit.v / rtl/mac8.v), and neural_processor.v's own
+header states its accumulation is "same sign-extended INT32-style
+accumulation" as that same MAC lineage.
+"""
+from __future__ import annotations
+
+from tools.validation.mac_oracle import to_signed, to_unsigned, mac8_tree
+
+INT8_MIN = -128
+INT8_MAX = 127
+DEFAULT_DATA_WIDTH = 8
+DEFAULT_ACC_WIDTH = 32
+DEFAULT_P_IN = 8
+
+ACT_NONE = "none"
+ACT_RELU = "relu"
+
+
+def check_int8(val: int, name: str = "value") -> int:
+    """Assert `val` is a valid signed INT8 and return it unchanged."""
+    if not (INT8_MIN <= val <= INT8_MAX):
+        raise ValueError(f"{name}={val} out of signed INT8 range [{INT8_MIN}, {INT8_MAX}]")
+    return val
+
+
+def wrap_acc(val: int, acc_width: int = DEFAULT_ACC_WIDTH) -> int:
+    """Wrap a Python int to signed acc_width-bit two's complement -- the
+    exact semantics of a fixed-width Verilog `reg signed [acc_width-1:0]`
+    after a `+` that would otherwise overflow."""
+    return to_signed(val, acc_width)
+
+
+def tile_product_sum(x_tile, w_tile, acc_width: int = DEFAULT_ACC_WIDTH,
+                      data_width: int = DEFAULT_DATA_WIDTH) -> int:
+    """
+    One P_IN-wide tile: P_IN independent INT8xINT8 products, reduced by
+    the balanced adder tree (mac8_tree, acc_in=0). Matches
+    neural_processor.v's `tile_sum` (stages 1..1+TREE_LEVELS) exactly.
+    """
+    if len(x_tile) != len(w_tile):
+        raise ValueError("x_tile and w_tile must have the same length (P_IN)")
+    n = len(x_tile)
+    if n == 0 or (n & (n - 1)) != 0:
+        raise ValueError(f"tile length {n} must be a power of two (P_IN), matching the RTL's tree")
+    products = []
+    for x, w in zip(x_tile, w_tile):
+        check_int8(x, "input")
+        check_int8(w, "weight")
+        products.append(x * w)  # exact, INT8xINT8 never overflows PROD_WIDTH=16
+    return mac8_tree(products, acc_in=0, acc_width=acc_width)
+
+
+def accumulate_tile(acc_reg: int, tile_sum: int, acc_width: int = DEFAULT_ACC_WIDTH) -> int:
+    """`acc_reg <= acc_reg + tile_sum` -- one running-accumulator update
+    across tiles of the SAME job. WRAPAROUND, matching the RTL's plain
+    fixed-width `+` exactly (not saturating)."""
+    return wrap_acc(acc_reg + tile_sum, acc_width)
+
+
+def add_bias(acc_reg: int, bias: int, acc_width: int = DEFAULT_ACC_WIDTH,
+             data_width: int = DEFAULT_DATA_WIDTH) -> int:
+    """`final_acc_reg <= acc_reg + bias_ext` -- bias is sign-extended from
+    an INT8 value, then added with WRAPAROUND semantics (same as
+    accumulate_tile)."""
+    check_int8(bias, "bias")
+    return wrap_acc(acc_reg + bias, acc_width)
+
+
+def activate_and_saturate(final_acc: int, activation: str = ACT_RELU,
+                           data_width: int = DEFAULT_DATA_WIDTH) -> int:
+    """
+    The ONE saturating stage in the whole datapath -- neural_processor.v
+    lines ~256-288, reproduced exactly (not approximated):
+
+    ACT_NONE: two-sided saturating clamp to [INT8_MIN, INT8_MAX] -- if
+      final_acc fits in signed data_width bits, pass its exact truncated
+      value through; otherwise clamp to INT8_MIN (if negative) or
+      INT8_MAX (if positive).
+
+    ACT_RELU (default, and the RTL's own `case` default for ANY
+      activation code other than exactly ACT_NONE): final_acc<=0 -> 0;
+      0 < final_acc <= INT8_MAX -> final_acc exactly; final_acc >
+      INT8_MAX -> saturate to INT8_MAX. There is no negative saturation
+      branch for ReLU since negative values are already zeroed.
+    """
+    lo, hi = -(1 << (data_width - 1)), (1 << (data_width - 1)) - 1
+
+    if activation == ACT_NONE:
+        if lo <= final_acc <= hi:
+            return final_acc
+        return lo if final_acc < 0 else hi
+
+    # ACT_RELU (and, matching the RTL's `default:` case branch, any
+    # activation value that isn't exactly ACT_NONE)
+    if final_acc <= 0:
+        return 0
+    if final_acc > hi:
+        return hi
+    return final_acc
+
+
+def neuron_reference(inputs, weights, bias: int = 0, activation: str = ACT_RELU,
+                      p_in: int = DEFAULT_P_IN, acc_width: int = DEFAULT_ACC_WIDTH,
+                      data_width: int = DEFAULT_DATA_WIDTH) -> int:
+    """
+    Full, tile-by-tile, bit-exact reference for one neuron_processor.v
+    job: y = activation(bias + sum(x[i]*w[i] for i in 0..N_INPUTS-1)),
+    computed the SAME WAY the RTL computes it -- P_IN-wide tiles, each
+    reduced by the adder tree, accumulated across tiles with 32-bit
+    wraparound, THEN bias-added (also wraparound), THEN activated/
+    saturated exactly once at the end (matching tile_last/NP_FINISH).
+
+    len(inputs) must be a multiple of p_in (one real tile per group of
+    p_in inputs -- matches n_tiles*P_IN in the real WRITE_JOB protocol).
+    """
+    if len(inputs) != len(weights):
+        raise ValueError("inputs and weights must have the same length")
+    if len(inputs) == 0 or len(inputs) % p_in != 0:
+        raise ValueError(f"len(inputs)={len(inputs)} must be a nonzero multiple of p_in={p_in}")
+
+    acc = 0
+    for t in range(0, len(inputs), p_in):
+        x_tile = inputs[t:t + p_in]
+        w_tile = weights[t:t + p_in]
+        tile_sum = tile_product_sum(x_tile, w_tile, acc_width=acc_width, data_width=data_width)
+        acc = accumulate_tile(acc, tile_sum, acc_width=acc_width)
+
+    final_acc = add_bias(acc, bias, acc_width=acc_width, data_width=data_width)
+    return activate_and_saturate(final_acc, activation=activation, data_width=data_width)
diff --git a/tools/neural_sim/tests/test_compare.py b/tools/neural_sim/tests/test_compare.py
new file mode 100644
index 0000000..2a55063
--- /dev/null
+++ b/tools/neural_sim/tests/test_compare.py
@@ -0,0 +1,50 @@
+import json
+
+import pytest
+
+from tools.neural_sim.compare import compare_results, load_fpga_results
+
+
+def test_exact_match():
+    report = compare_results([1, 2, 3], [1, 2, 3])
+    assert report.exact_match
+    assert report.num_mismatches == 0
+
+
+def test_single_mismatch_reported_precisely():
+    report = compare_results([1, 2, 3], [1, 5, 3])
+    assert not report.exact_match
+    assert report.num_mismatches == 1
+    assert report.first_mismatch_index == 1
+    assert report.first_mismatch_expected == 2
+    assert report.first_mismatch_actual == 5
+    assert report.max_abs_diff == 3
+
+
+def test_multiple_mismatches_max_abs_diff():
+    report = compare_results([0, 0, 0], [10, -5, 0])
+    assert report.num_mismatches == 2
+    assert report.max_abs_diff == 10
+
+
+def test_length_mismatch_raises():
+    with pytest.raises(ValueError):
+        compare_results([1, 2], [1, 2, 3])
+
+
+def test_load_fpga_results_json_list(tmp_path):
+    path = tmp_path / "results.json"
+    path.write_text(json.dumps([1, -2, 3]))
+    assert load_fpga_results(str(path)) == [1, -2, 3]
+
+
+def test_load_fpga_results_json_dict_with_results_key(tmp_path):
+    path = tmp_path / "results.json"
+    path.write_text(json.dumps({"results": [4, 5, 6], "meta": "x"}))
+    assert load_fpga_results(str(path)) == [4, 5, 6]
+
+
+def test_load_fpga_results_plain_text(tmp_path):
+    path = tmp_path / "results.txt"
+    path.write_text("1 2 -3\n4 5")
+    assert load_fpga_results(str(path)) == [1, 2, -3, 4, 5]
diff --git a/tools/neural_sim/tests/test_layer.py b/tools/neural_sim/tests/test_layer.py
new file mode 100644
index 0000000..3204fb3
--- /dev/null
+++ b/tools/neural_sim/tests/test_layer.py
@@ -0,0 +1,44 @@
+import numpy as np
+import pytest
+
+from tools.neural_sim.layer import FCLayer
+from tools.neural_sim.neuron import neuron_scalar
+
+
+def test_layer_matches_per_neuron_scalar_reference():
+    rng = np.random.default_rng(7)
+    n_inputs, n_neurons = 32, 5
+    weights = rng.integers(-128, 128, size=(n_neurons, n_inputs), dtype=np.int64)
+    biases = rng.integers(-128, 128, size=n_neurons, dtype=np.int64)
+    inputs = rng.integers(-128, 128, size=n_inputs, dtype=np.int64)
+
+    layer = FCLayer(weights, biases=biases, activation="relu")
+    outputs = layer.forward(inputs)
+
+    for n in range(n_neurons):
+        expected = neuron_scalar(inputs.tolist(), weights[n].tolist(), bias=int(biases[n]), activation="relu")
+        assert outputs[n] == expected
+
+
+def test_layer_rejects_bad_input_shape():
+    layer = FCLayer([[1] * 8])
+    with pytest.raises(ValueError):
+        layer.forward([1] * 7)
+
+
+def test_layer_rejects_non_multiple_of_p_in():
+    with pytest.raises(ValueError):
+        FCLayer([[1] * 7])
+
+
+def test_layer_default_bias_is_zero():
+    layer = FCLayer([[1] * 8, [2] * 8])
+    assert layer.biases.tolist() == [0, 0]
+
+
+def test_layer_all_outputs_in_int8_range():
+    rng = np.random.default_rng(99)
+    weights = rng.integers(-128, 128, size=(20, 64), dtype=np.int64)
+    inputs = rng.integers(-128, 128, size=64, dtype=np.int64)
+    outputs = FCLayer(weights).forward(inputs)
+    assert all(-128 <= v <= 127 for v in outputs.tolist())
diff --git a/tools/neural_sim/tests/test_memory.py b/tools/neural_sim/tests/test_memory.py
new file mode 100644
index 0000000..e05cdcb
--- /dev/null
+++ b/tools/neural_sim/tests/test_memory.py
@@ -0,0 +1,63 @@
+import pytest
+
+from tools.neural_sim.memory import (
+    MemoryModel, WEIGHTS_BASE, ACTIVATIONS_BASE, RESULTS_BASE, SDRAM_SIZE,
+)
+
+
+def test_read_after_write_each_region():
+    mem = MemoryModel()
+    mem.write_weight(0, -5)
+    mem.write_activation(0, 42)
+    mem.write_result(0, -128)
+    assert mem.read_weight(0) == -5
+    assert mem.read_activation(0) == 42
+    assert mem.read_result(0) == -128
+
+
+def test_regions_are_at_the_real_v2_addresses():
+    mem = MemoryModel()
+    mem.write_weight(0, 1)
+    mem.write_activation(0, 2)
+    mem.write_result(0, 3)
+    assert mem.read_byte(WEIGHTS_BASE) == 1
+    assert mem.read_byte(ACTIVATIONS_BASE) == 2
+    assert mem.read_byte(RESULTS_BASE) == 3
+
+
+def test_adjacent_regions_do_not_corrupt_each_other():
+    mem = MemoryModel()
+    mem.write_weight(0x0FFFFF - WEIGHTS_BASE, 0x11)  # last word before activations
+    mem.write_activation(0, 0x22)                     # first word of activations
+    assert mem.read_weight(0x0FFFFF - WEIGHTS_BASE) == 0x11
+    assert mem.read_activation(0) == 0x22
+
+
+def test_out_of_range_byte_raises():
+    mem = MemoryModel()
+    with pytest.raises(IndexError):
+        mem.read_byte(-1)
+    with pytest.raises(IndexError):
+        mem.read_byte(SDRAM_SIZE)
+
+
+def test_value_out_of_int8_range_raises():
+    mem = MemoryModel()
+    with pytest.raises(ValueError):
+        mem.write_byte(0, 128)
+    with pytest.raises(ValueError):
+        mem.write_byte(0, -129)
+
+
+def test_region_bounds_checking_rejects_spillover():
+    mem = MemoryModel()
+    weights_size = ACTIVATIONS_BASE - WEIGHTS_BASE
+    with pytest.raises(IndexError):
+        mem.write_weight(weights_size, 0)  # one byte past the weights region
+
+
+def test_bulk_helpers_round_trip():
+    mem = MemoryModel()
+    values = list(range(-10, 10))
+    mem.write_weights(0, values)
+    assert mem.read_weights(0, len(values)) == values
diff --git a/tools/neural_sim/tests/test_network_examples.py b/tools/neural_sim/tests/test_network_examples.py
new file mode 100644
index 0000000..3150507
--- /dev/null
+++ b/tools/neural_sim/tests/test_network_examples.py
@@ -0,0 +1,36 @@
+import numpy as np
+import pytest
+
+from tools.neural_sim import examples as ex
+from tools.neural_sim.layer import FCLayer
+from tools.neural_sim.network import Network
+
+
+def test_network_rejects_empty():
+    with pytest.raises(ValueError):
+        Network([])
+
+
+def test_network_rejects_shape_mismatch():
+    layer1 = FCLayer([[1] * 8])           # 8 -> 1
+    layer2 = FCLayer([[1] * 8, [2] * 8])  # 8 -> 2, but layer1 outputs only 1
+    with pytest.raises(ValueError):
+        Network([layer1, layer2])
+
+
+@pytest.mark.parametrize("name", list(ex.ALL_EXAMPLES.keys()))
+def test_all_examples_run_and_stay_in_int8_range(name):
+    net = ex.ALL_EXAMPLES[name]()
+    rng = np.random.default_rng(123)
+    inputs = rng.integers(-128, 128, size=net.n_inputs, dtype=np.int64)
+    outputs = net.forward(inputs)
+    assert len(outputs) == net.n_outputs
+    assert all(-128 <= v <= 127 for v in outputs.tolist())
+
+
+def test_network_forward_all_matches_forward_final():
+    net = ex.example_8_8_1()
+    inputs = np.array([1, 2, 3, 4, -1, -2, -3, -4], dtype=np.int64)
+    all_outputs = net.forward_all(inputs)
+    assert all_outputs[-1].tolist() == net.forward(inputs).tolist()
+    assert len(all_outputs) == len(net.layers)
diff --git a/tools/neural_sim/tests/test_neuron.py b/tools/neural_sim/tests/test_neuron.py
new file mode 100644
index 0000000..c8b406a
--- /dev/null
+++ b/tools/neural_sim/tests/test_neuron.py
@@ -0,0 +1,30 @@
+import numpy as np
+import pytest
+
+from tools.neural_sim.neuron import neuron_scalar, neuron_vectorized
+from tools.neural_sim.numerics import ACT_NONE, ACT_RELU
+
+
+@pytest.mark.parametrize("n_tiles", [1, 2, 4, 16])
+@pytest.mark.parametrize("activation", [ACT_RELU, ACT_NONE])
+@pytest.mark.parametrize("seed", [0, 1, 2, 3, 4])
+def test_scalar_and_vectorized_agree(n_tiles, activation, seed):
+    rng = np.random.default_rng(seed)
+    n = n_tiles * 8
+    inputs = rng.integers(-128, 128, size=n, dtype=np.int64).tolist()
+    weights = rng.integers(-128, 128, size=n, dtype=np.int64).tolist()
+    bias = int(rng.integers(-128, 128))
+
+    scalar = neuron_scalar(inputs, weights, bias=bias, activation=activation)
+    vectorized = neuron_vectorized(inputs, weights, bias=bias, activation=activation)
+    assert scalar == vectorized
+
+
+def test_scalar_and_vectorized_agree_on_extremes():
+    inputs = [-128, 127] * 4
+    weights = [127, -128] * 4
+    assert neuron_scalar(inputs, weights) == neuron_vectorized(inputs, weights)
+
+
+def test_scalar_and_vectorized_agree_on_zero():
+    assert neuron_scalar([0] * 8, [0] * 8) == neuron_vectorized([0] * 8, [0] * 8) == 0
diff --git a/tools/neural_sim/tests/test_numerics.py b/tools/neural_sim/tests/test_numerics.py
new file mode 100644
index 0000000..e8c6e2d
--- /dev/null
+++ b/tools/neural_sim/tests/test_numerics.py
@@ -0,0 +1,116 @@
+import pytest
+
+from tools.neural_sim.numerics import (
+    check_int8, wrap_acc, tile_product_sum, accumulate_tile, add_bias,
+    activate_and_saturate, neuron_reference, ACT_NONE, ACT_RELU,
+)
+
+
+def test_check_int8_accepts_range():
+    assert check_int8(-128) == -128
+    assert check_int8(127) == 127
+    assert check_int8(0) == 0
+
+
+def test_check_int8_rejects_out_of_range():
+    with pytest.raises(ValueError):
+        check_int8(128)
+    with pytest.raises(ValueError):
+        check_int8(-129)
+
+
+def test_wrap_acc_no_overflow_is_identity():
+    assert wrap_acc(1000, acc_width=32) == 1000
+    assert wrap_acc(-1000, acc_width=32) == -1000
+
+
+def test_wrap_acc_true_32bit_wraparound():
+    # 2**31 is one past the max positive signed 32-bit value (2**31 - 1)
+    # -- must wrap to the most-negative value, exactly like a Verilog
+    # `reg signed [31:0]` silently overflowing.
+    assert wrap_acc(2**31, acc_width=32) == -(2**31)
+    assert wrap_acc(2**31 - 1, acc_width=32) == 2**31 - 1  # exact boundary, no wrap
+    assert wrap_acc(-(2**31) - 1, acc_width=32) == 2**31 - 1
+
+
+def test_tile_product_sum_exact_known_values():
+    # 1*1 + 2*1 + ... + 8*1 = 36
+    assert tile_product_sum(list(range(1, 9)), [1] * 8) == 36
+
+
+def test_tile_product_sum_extreme_product():
+    # -128 * -128 = 16384, the one INT8xINT8 case that does not fit
+    # symmetrically in magnitude terms
+    assert tile_product_sum([-128], [-128], acc_width=32) == 16384
+
+
+def test_tile_product_sum_rejects_non_power_of_two():
+    with pytest.raises(ValueError):
+        tile_product_sum([1, 2, 3], [1, 1, 1])
+
+
+def test_tile_product_sum_rejects_out_of_range_input():
+    with pytest.raises(ValueError):
+        tile_product_sum([200], [1])
+
+
+def test_accumulate_tile_matches_wrap_acc():
+    assert accumulate_tile(10, 20) == 30
+    assert accumulate_tile(2**31 - 1, 1) == -(2**31)
+
+
+def test_add_bias_wraparound():
+    assert add_bias(100, 27) == 127
+    assert add_bias(2**31 - 1, 127) == wrap_acc(2**31 - 1 + 127)
+
+
+def test_activate_relu_zeroes_non_positive():
+    assert activate_and_saturate(0, activation=ACT_RELU) == 0
+    assert activate_and_saturate(-1, activation=ACT_RELU) == 0
+    assert activate_and_saturate(-1000000, activation=ACT_RELU) == 0
+
+
+def test_activate_relu_passes_in_range():
+    assert activate_and_saturate(1, activation=ACT_RELU) == 1
+    assert activate_and_saturate(127, activation=ACT_RELU) == 127
+
+
+def test_activate_relu_saturates_positive():
+    assert activate_and_saturate(128, activation=ACT_RELU) == 127
+    assert activate_and_saturate(1000000, activation=ACT_RELU) == 127
+
+
+def test_activate_none_passes_full_signed_range():
+    assert activate_and_saturate(-128, activation=ACT_NONE) == -128
+    assert activate_and_saturate(127, activation=ACT_NONE) == 127
+    assert activate_and_saturate(0, activation=ACT_NONE) == 0
+
+
+def test_activate_none_saturates_both_sides():
+    assert activate_and_saturate(128, activation=ACT_NONE) == 127
+    assert activate_and_saturate(-129, activation=ACT_NONE) == -128
+    assert activate_and_saturate(1000000, activation=ACT_NONE) == 127
+    assert activate_and_saturate(-1000000, activation=ACT_NONE) == -128
+
+
+def test_neuron_reference_simple_positive():
+    y = neuron_reference(list(range(1, 9)), [1] * 8, activation=ACT_RELU)
+    assert y == 36
+
+
+def test_neuron_reference_multi_tile_accumulates_across_tiles():
+    # two tiles of 8, same weights -- accumulator must carry across tiles
+    inputs = [1] * 8 + [1] * 8
+    weights = [1] * 8 + [1] * 8
+    assert neuron_reference(inputs, weights, activation=ACT_RELU) == 16
+
+
+def test_neuron_reference_rejects_length_not_multiple_of_p_in():
+    with pytest.raises(ValueError):
+        neuron_reference([1] * 5, [1] * 5)
+
+
+def test_neuron_reference_bias_default_zero_matches_no_bias():
+    y_default = neuron_reference([1] * 8, [1] * 8)
+    y_explicit = neuron_reference([1] * 8, [1] * 8, bias=0)
+    assert y_default == y_explicit
diff --git a/tools/neural_sim/tests/test_vectors.py b/tools/neural_sim/tests/test_vectors.py
new file mode 100644
index 0000000..fc5db3a
--- /dev/null
+++ b/tools/neural_sim/tests/test_vectors.py
@@ -0,0 +1,66 @@
+import json
+import os
+
+from tools.neural_sim import vectors as vec
+from tools.neural_sim.layer import FCLayer
+
+
+def test_simple_positive_manually_predictable():
+    v = vec.gen_simple_positive()
+    # inputs=[1..8], weights=all 1 -> sum = 36, ReLU(36)=36
+    assert v.expected == [36]
+
+
+def test_zero_vector_is_all_zero_output():
+    v = vec.gen_zero()
+    assert v.expected == [0]
+
+
+def test_extremes_vector_is_self_consistent():
+    v = vec.gen_extremes()
+    layer = FCLayer(v.weights, biases=v.biases, activation=v.activation, p_in=v.p_in)
+    assert layer.forward(v.inputs).tolist() == v.expected
+
+
+def test_random_vector_is_deterministic_across_calls():
+    v1 = vec.gen_random(seed=555)
+    v2 = vec.gen_random(seed=555)
+    assert v1.inputs == v2.inputs
+    assert v1.weights == v2.weights
+    assert v1.expected == v2.expected
+
+
+def test_random_vector_different_seed_differs():
+    v1 = vec.gen_random(seed=1)
+    v2 = vec.gen_random(seed=2)
+    assert v1.inputs != v2.inputs
+
+
+def test_d_stress_dimensions_match_the_real_rtl_benchmark():
+    v = vec.gen_dstress()
+    assert v.n_neurons == 256
+    assert v.n_inputs == 128
+
+
+def test_d_stress_is_deterministic():
+    v1 = vec.gen_dstress(seed=42)
+    v2 = vec.gen_dstress(seed=42)
+    assert v1.expected == v2.expected
+
+
+def test_export_then_import_round_trip(tmp_path):
+    v = vec.gen_signed_mix()
+    path = str(tmp_path / "vec.json")
+    vec.export_json(v, path)
+    loaded = vec.load_json(path)
+    assert loaded == v
+
+
+def test_export_all_json_contains_every_generator(tmp_path):
+    path = str(tmp_path / "all.json")
+    vec.export_all_json(path)
+    with open(path) as f:
+        data = json.load(f)
+    assert set(data.keys()) == set(vec.ALL_GENERATORS.keys())
+    for name, d in data.items():
+        assert "expected" in d and "weights" in d and "inputs" in d
diff --git a/tools/neural_sim/vectors.py b/tools/neural_sim/vectors.py
new file mode 100644
index 0000000..1fbfd7f
--- /dev/null
+++ b/tools/neural_sim/vectors.py
@@ -0,0 +1,144 @@
+"""
+Deterministic FPGA-style test vectors + golden-vector JSON export
+(B8/B9). Every generator below returns a TestVector whose `expected`
+field is computed by this SAME package's own golden model
+(layer.FCLayer / numerics.neuron_reference) -- i.e. the vector is
+self-consistent and re-verifiable by construction (see
+tests/test_vectors.py's own round-trip test), not hand-typed.
+"""
+from __future__ import annotations
+
+import json
+from dataclasses import dataclass, field, asdict
+from typing import List, Optional
+
+import numpy as np
+
+from .numerics import ACT_RELU, DEFAULT_P_IN
+from .layer import FCLayer
+
+
+@dataclass
+class TestVector:
+    name: str
+    inputs: List[int]
+    weights: List[List[int]]     # shape (n_neurons, n_inputs)
+    biases: List[int]            # shape (n_neurons,)
+    activation: str
+    p_in: int
+    expected: List[int]          # shape (n_neurons,)
+    numeric_format: str = "int8 in / int8 weight / int32 accumulate (wraparound) / int8 out (saturate)"
+    seed: Optional[int] = None
+
+    @property
+    def n_neurons(self) -> int:
+        return len(self.expected)
+
+    @property
+    def n_inputs(self) -> int:
+        return len(self.inputs)
+
+
+def _make_vector(name: str, inputs, weights, biases=None, activation: str = ACT_RELU,
+                  p_in: int = DEFAULT_P_IN, seed: Optional[int] = None) -> TestVector:
+    layer = FCLayer(weights, biases=biases, activation=activation, p_in=p_in)
+    expected = layer.forward(inputs).tolist()
+    biases_list = layer.biases.tolist()
+    return TestVector(
+        name=name, inputs=list(int(v) for v in inputs),
+        weights=[[int(w) for w in row] for row in layer.weights.tolist()],
+        biases=biases_list, activation=activation, p_in=p_in,
+        expected=[int(e) for e in expected], seed=seed,
+    )
+
+
+# ---- Test 1: simple positive ----
+def gen_simple_positive(p_in: int = DEFAULT_P_IN) -> TestVector:
+    """Small positive inputs/weights, manually predictable:
+    inputs=[1..p_in], weights=all 1s -> sum = p_in*(p_in+1)/2."""
+    inputs = list(range(1, p_in + 1))
+    weights = [[1] * p_in]
+    return _make_vector("simple_positive", inputs, weights)
+
+
+# ---- Test 2: signed values ----
+def gen_signed_mix(p_in: int = DEFAULT_P_IN) -> TestVector:
+    """Alternating positive/negative INT8 inputs and weights."""
+    inputs = [((-1) ** i) * (10 + i) for i in range(p_in)]
+    weights = [[((-1) ** (i + 1)) * (5 + i) for i in range(p_in)]]
+    return _make_vector("signed_mix", inputs, weights)
+
+
+# ---- Test 3: extremes ----
+def gen_extremes(p_in: int = DEFAULT_P_IN) -> TestVector:
+    """-128/127 combinations designed to stress multiplication (the one
+    genuinely asymmetric INT8xINT8 case, -128*-128=16384, is
+    deliberately included) and accumulation across p_in terms."""
+    inputs = [(-128 if i % 2 == 0 else 127) for i in range(p_in)]
+    weights = [[(-128 if i % 2 == 1 else 127) for i in range(p_in)]]
+    return _make_vector("extremes", inputs, weights)
+
+
+# ---- Test 4: zero ----
+def gen_zero(p_in: int = DEFAULT_P_IN) -> TestVector:
+    inputs = [0] * p_in
+    weights = [[0] * p_in]
+    return _make_vector("zero", inputs, weights)
+
+
+# ---- Test 5: random (fixed seed) ----
+def gen_random(seed: int = 1234, n_inputs: int = 64, n_neurons: int = 4,
+               p_in: int = DEFAULT_P_IN) -> TestVector:
+    rng = np.random.default_rng(seed)
+    inputs = rng.integers(-128, 128, size=n_inputs, dtype=np.int64)
+    weights = rng.integers(-128, 128, size=(n_neurons, n_inputs), dtype=np.int64)
+    return _make_vector(f"random_seed{seed}", inputs, weights, p_in=p_in, seed=seed)
+
+
+# ---- Test 6: D-Stress (256 neurons x 128 inputs, reproducing the
+# existing RTL D-Stress benchmark's own dimensions) ----
+def gen_dstress(seed: int = 42, n_neurons: int = 256, n_inputs: int = 128,
+                 p_in: int = DEFAULT_P_IN) -> TestVector:
+    rng = np.random.default_rng(seed)
+    inputs = rng.integers(-128, 128, size=n_inputs, dtype=np.int64)
+    weights = rng.integers(-128, 128, size=(n_neurons, n_inputs), dtype=np.int64)
+    return _make_vector("d_stress", inputs, weights, p_in=p_in, seed=seed)
+
+
+ALL_GENERATORS = {
+    "simple_positive": gen_simple_positive,
+    "signed_mix": gen_signed_mix,
+    "extremes": gen_extremes,
+    "zero": gen_zero,
+    "random": gen_random,
+    "d_stress": gen_dstress,
+}
+
+
+def generate_all() -> dict:
+    return {name: fn() for name, fn in ALL_GENERATORS.items()}
+
+
+# ---- golden-vector export/import (machine-readable, JSON) ----
+def to_dict(vector: TestVector) -> dict:
+    return asdict(vector)
+
+
+def export_json(vector: TestVector, path: str) -> None:
+    with open(path, "w") as f:
+        json.dump(to_dict(vector), f, indent=2)
+
+
+def load_json(path: str) -> TestVector:
+    with open(path) as f:
+        d = json.load(f)
+    return TestVector(**d)
+
+
+def export_all_json(path: str) -> None:
+    """Export every named generator's vector into one JSON file, keyed
+    by name -- convenient for a future RTL testbench harness to load
+    once and iterate."""
+    vectors = generate_all()
+    with open(path, "w") as f:
+        json.dump({name: to_dict(v) for name, v in vectors.items()}, f, indent=2)