Perf: vectorised element-wise buffer loops keeping per-element atomicity (#226)

benchmarks/typedpython/results/2026-10-25-buffer-vec-scale.txt

@@ -0,0 +1,8 @@
+# scale_arr (a[i] = a[i] * k), n=2,000,000 doubles, 3.15.0rc1t, aarch64, 11 interleaved rounds, median; ONE quick run
+# load before 13.15, after 12.82 (NOT quiet; the 2026-10-17 quiet-window run measured atomic/plain 2.54)
+plain 0.990 ms, atomic scalar (-DTP_NO_BUF_VEC) 1.162 ms, vectorised 0.900 ms
+atomic/plain 1.173, vectorised/plain 0.909
+# re-run with the inline-asm helper (ldr/str q), load 8.9: plain 0.782, atomic scalar 0.981, vectorised 0.962 ms
+atomic/plain 1.255, vectorised/plain 1.230  (the 0.91 above was the NEON-intrinsic helper, withdrawn)
+# third run: precise "Q" memory operands, no clobber, 4 elements per step; load 8.73 -> 8.73, one run
+plain 0.684 ms, atomic scalar 0.853 ms, vectorised 0.646 ms; atomic/plain 1.246, vectorised/plain 0.945

docs/design/typedpython.md

@@ -402,6 +402,32 @@ i64/f64 지역)는 **컴파일된 구간 안의 임시 값**만 줄인다. 컨
 **선행 조건:** 스텁 생성기가 픽스처 대역이 아닌 실제 `pythonx-map.toml` 을 읽고, 생성된 스텁이
 IDE·타입 체커가 찾는 자리(pythonx-compose 배포본)에 설치돼야 한다. 2026-10 현재 미완료.
 
+### 4.x 벡터화된 버퍼 루프 (#226)
+
+`array.array[float]` 버퍼 인자(#149)의 원소 접근은 원소 하나당 relaxed atomic 8바이트 접근이다. 컴파일러는 atomic 접근을
+벡터화하지 않으므로 `a[i] = a[i] * k` 형태의 저장 루프가 느렸다(측정: 평범한 접근 대비 2.54배,
+`benchmarks/typedpython/results/2026-10-17-buffer-atomic-vs-plain.txt`). 해결은 두 레인 벡터 접근 도우미다.
+
+- **적용 모양**: `for i in range(a, b)` (step 1), 본문이 문장 하나, 증명된(`proven`) `x[i] = f(...)` 로 F64 버퍼에 저장하고,
+  `f` 는 같은 `i` 에서 증명된 F64 버퍼 읽기, F64 지역, 상수, `+ - *`, 단항 `-` 만 쓴다. 나눗셈(예외 가능), I64(검사된 오버플로),
+  증명되지 않은 읽기는 스칼라 경로를 유지한다. 벡터 루프는 두 원소씩 처리하고 나머지는 기존 스칼라 루프(폴링 청크 안)가 맡는다.
+- **도우미**: `tp_buf_ld2_f64` / `tp_buf_st2_f64` (`tp_runtime.h`), `TP_BUF_VEC2` 가 정의된 빌드에서만 존재한다.
+  `-DTP_NO_BUF_VEC` 는 도우미를 끄고(스칼라 relaxed atomic 으로 되돌림), `-DTP_BUFFER_PLAIN` 빌드에도 없다.
+- **대상별로 원소 접근이 단일 복사 원자적인 이유 (8바이트 미만으로 찢어지지 않음)**:
+  - aarch64: Arm ARM (DDI 0487), B2.2.1 "Requirements for single-copy atomicity": a read/write of a 128-bit SIMD&FP register
+    that is 64-bit aligned is treated as a pair of single-copy atomic 64-bit accesses. The helper is inline asm (`ldr q` / `str q`
+    with precise "Q" memory operands, no clobber), so the compiler sees no non-atomic C access and cannot split, merge or reorder it as a plain one;
+    the per-element guarantee comes from the architecture. It needs 8-byte alignment, which `tp_buf_enter` already checks.
+    **Check the exact wording against the manual revision**: this environment had no copy of the manual, the summary is from memory.
+  - x86-64: SDM Vol. 3A 8.1.1 이 보장하는 것은 64비트 경계에 정렬된 쿼드워드 읽기/쓰기까지다. 128비트 SSE/AVX 접근의
+    원소별 원자성은 SDM 이 보장하지 않으므로(벤더 별도 문서에만 있음) **인용할 보장이 없어** 구현하지 않았다. x86-64 는
+    스칼라 relaxed atomic 경로를 쓴다.
+  - 그 밖의 대상(wasm 등): 스칼라 relaxed atomic.
+- **C 메모리 모델**: 벡터 접근은 인라인 어셈블리라 C 컴파일러가 보는 비원자 C 접근이 없다. 보장은 C 표준이 아니라 위 아키텍처
+  규칙에 있고, 도우미 한 곳에 모아 두어 대상 판별이 한 줄(`#if`)로 드러난다. N-9 의 의미(찢어진 값 없음)는 유지된다.
+- **측정** (`benchmarks/typedpython/results/2026-10-25-buffer-vec-scale.txt`): scale 커널, 인라인 어셈블리 도우미로 로드 약 9 에서 평범한 접근 대비
+  atomic 스칼라 1.25배, 새 경로(정밀 메모리 피연산자, 4원소 언롤) 0.94배.
+
 ## 5. 등급 2, 별도 언어 없는 커널
 
 > **결정 (2026-10-04, 메인테이너):** GPU 커널과 CPU SIMD 는 탐색이 아니라 **약속**이다 (INTENT §1.5).

python-multiplatform-ksp/src/main/python/typedpython/cgen.py

@@ -1903,6 +1903,83 @@ def st_While(self, s: ir.While) -> None:
         self.depth -= 1
         self.emit("}")
 
+    def _vec_plan(self, s: ir.ForRange):
+        """The element-wise buffer store loop `for i in range(a, b): x[i] = f(y[i], ..., k)` that
+        may run two F64 lanes at a time (TP_BUF_VEC2, tp_runtime.h; design/typedpython.md
+        "Vectorised buffer loops"). Returns the store statement, or None. Only step 1, one
+        statement, a proven store at the loop variable into an F64 buffer array, and a value built
+        from proven loads at the same variable, F64 locals and constants, + - * and negation.
+        No division (it raises), no I64 (checked overflow), no other statement."""
+        if s.step.value != 1 or len(s.body) != 1 or not isinstance(s.body[0], ir.StoreIndex):
+            return None
+        st = s.body[0]
+
+        def at_var(ix, proven):
+            return proven and isinstance(ix, ir.Local) and ix.name == s.var
+
+        def ok(e):
+            if e.type is not Type.F64:
+                return False
+            if isinstance(e, ir.Const):
+                return True
+            if isinstance(e, ir.Local):
+                return e.name != s.var and self.types.get(e.name) is Type.F64
+            if isinstance(e, ir.Index):
+                a = self.arrays.get(e.array)
+                return (a is not None and a.buffer and a.type is Type.F64_ARRAY
+                        and at_var(e.index, e.proven))
+            if isinstance(e, ir.BinOp):
+                return (e.op in (ir.BinOpKind.ADD, ir.BinOpKind.SUB, ir.BinOpKind.MUL)
+                        and ok(e.left) and ok(e.right))
+            if isinstance(e, ir.UnaryOp):
+                return e.op is ir.UnaryOpKind.NEG and ok(e.operand)
+            return False
+
+        a = self.arrays.get(st.array)
+        if a is None or not a.buffer or a.type is not Type.F64_ARRAY or not a.stored:
+            return None
+        if not at_var(st.index, st.proven) or not ok(st.value):
+            return None
+        return st
+
+    def _vec_expr(self, e: ir.Expr, base: str, hoist: list[str]) -> str:
+        if isinstance(e, (ir.Const, ir.Local)):
+            v = self.tmp("double")
+            hoist.append(f"{v} = {self.ev(e)};")
+            return f"vdupq_n_f64({v})"
+        if isinstance(e, ir.Index):
+            _, ptr = self._array(e.array)
+            return f"tp_buf_ld2_f64({ptr}, {base})"
+        if isinstance(e, ir.UnaryOp):
+            return f"vnegq_f64({self._vec_expr(e.operand, base, hoist)})"
+        fn = {ir.BinOpKind.ADD: "vaddq_f64", ir.BinOpKind.SUB: "vsubq_f64",
+              ir.BinOpKind.MUL: "vmulq_f64"}[e.op]
+        return f"{fn}({self._vec_expr(e.left, base, hoist)}, {self._vec_expr(e.right, base, hoist)})"
+
+    def _emit_vec_loop(self, st: ir.StoreIndex, ts: str, tk: str, limit: str) -> None:
+        """Four lanes per step (two independent 2-lane loads, then two stores) while tk + 4 <=
+        limit, then two lanes while tk + 2 <= limit; the scalar loop that follows takes the rest."""
+        hoist: list[str] = []
+        self.emit("#ifdef TP_BUF_VEC2")
+        base = f"(Py_ssize_t)((uint64_t){ts} + {tk})"
+        base2 = f"(Py_ssize_t)((uint64_t){ts} + {tk} + 2u)"
+        e1 = self._vec_expr(st.value, base, hoist)
+        e2 = self._vec_expr(st.value, base2, hoist)
+        for h in hoist:
+            self.emit(h)
+        _, ptr = self._array(st.array)
+        v1, v2 = "tp_vv1", "tp_vv2"
+        self.emit(f"for (; {tk} + 4u <= {limit}; {tk} += 4u) {{")
+        self.emit(f"    float64x2_t {v1} = {e1};")
+        self.emit(f"    float64x2_t {v2} = {e2};")
+        self.emit(f"    tp_buf_st2_f64({ptr}, {base}, {v1});")
+        self.emit(f"    tp_buf_st2_f64({ptr}, {base2}, {v2});")
+        self.emit("}")
+        self.emit(f"for (; {tk} + 2u <= {limit}; {tk} += 2u) {{")
+        self.emit(f"    tp_buf_st2_f64({ptr}, {base}, {e1});")
+        self.emit("}")
+        self.emit("#endif")
+
     def st_ForRange(self, s: ir.ForRange) -> None:
         if self.types.get(s.var) is not Type.I64:
             raise CGenError(f"{self.f.name}: range variable {s.var} must be an I64 local")
@@ -1931,9 +2008,13 @@ def st_ForRange(self, s: ir.ForRange) -> None:
             self.emit(f"    {tn} = ((uint64_t){ts} - (uint64_t){te} - 1u) / {mag} + 1u;")
             sign = "-"
         self.depth += 1
+        vec = self._vec_plan(s) if step == 1 else None
         if self.snapshot:
             # no polling in a snapshot holder: a plain loop
-            self.emit(f"for ({tk} = 0; {tk} < {tn}; {tk}++) {{")
+            self.emit(f"{tk} = 0;")
+            if vec is not None:
+                self._emit_vec_loop(vec, ts, tk, tn)
+            self.emit(f"for (; {tk} < {tn}; {tk}++) {{")
             self.depth += 1
             self.brk.append(None)
         else:
@@ -1945,6 +2026,8 @@ def st_ForRange(self, s: ir.ForRange) -> None:
             self.depth += 1
             self.poll()
             self.emit(f"{te2} = {tn} - {tk} > (uint64_t)TP_POLL_INTERVAL ? {tk} + (uint64_t)TP_POLL_INTERVAL : {tn};")
+            if vec is not None:
+                self._emit_vec_loop(vec, ts, tk, te2)
             self.emit(f"for (; {tk} < {te2}; {tk}++) {{")
             self.depth += 1
             self.brk.append(label)

python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h

@@ -1240,6 +1240,33 @@ static inline void tp_buf_st_i64(const tp_buf *b, Py_ssize_t i, int64_t v)
     memcpy(&u, &v, sizeof u);
     __atomic_store_n((uint64_t *)b->data + i, u, __ATOMIC_RELAXED);
 }
+
+/* Two-lane F64 access for the vectorised element-wise store loop (cgen.py _vec_plan; issue #226;
+ * design/typedpython.md "Vectorised buffer loops"). aarch64 only: LDR/STR Qt on an address
+ * aligned to the 8-byte element (tp_buf_enter guarantees 8-byte alignment of the data) is two
+ * element accesses, each single-copy atomic (Arm ARM DDI 0487 B2.2.1: a 64-bit aligned read/write of
+ * a 128-bit SIMD&FP register is a pair of single-copy atomic 64-bit accesses; check the wording
+ * against the manual revision). Written as inline asm (ldr/str q), not NEON intrinsics, so the
+ * compiler sees no plain C access it could split, merge or reorder. Neither lane is
+ * ever split below 8 bytes. Other targets (x86-64 included: the Intel SDM Vol. 3A 8.1.1 lists
+ * 128-bit SSE accesses as NOT guaranteed atomic) get no TP_BUF_VEC2 and keep the scalar relaxed
+ * atomic loop above. */
+#if defined(__aarch64__) && defined(__ARM_NEON) && !defined(TP_NO_BUF_VEC)
+#include <arm_neon.h>
+#define TP_BUF_VEC2 1
+static inline float64x2_t tp_buf_ld2_f64(const tp_buf *b, Py_ssize_t i)
+{
+    float64x2_t v;                  /* inline asm: no C-level non-atomic access exists to split or merge */
+    const double *p = (const double *)b->data + i;
+    __asm__ volatile("ldr %q0, %1" : "=w"(v) : "Q"(*(const double (*)[2])p));
+    return v;
+}
+static inline void tp_buf_st2_f64(const tp_buf *b, Py_ssize_t i, float64x2_t v)
+{
+    double *p = (double *)b->data + i;
+    __asm__ volatile("str %q1, %0" : "=Q"(*(double (*)[2])p) : "w"(v));
+}
+#endif
 #endif
 
 /* ------------------------------------------------------------------------------------------

python-multiplatform-ksp/src/test/python/test_buffer_arrays.py

@@ -496,3 +496,142 @@ def many(n, buffer):
 
     assert_accepted(many(5, "array"))
     assert_rejected("verify/structure", many(3, None))
+
+
+# --- vectorised element-wise store loops (#226) ---------------------------------------------------
+
+VEC_SOURCE = textwrap.dedent('''
+    # typedpython: compiled
+    import array
+
+
+    def vscale(a: array.array[float], k: float) -> None:
+        for i in range(len(a)):
+            a[i] = a[i] * k
+
+
+    def vmix(a: array.array[float], k: float) -> None:
+        for i in range(len(a)):
+            a[i] = -(a[i] - k) * 2.5 + k
+
+
+    def vtwo(a: array.array[float], b: array.array[float]) -> None:
+        for i in range(len(a)):
+            a[i] = a[i] + b[i]
+
+
+    def vfill(a: array.array[float], k: float) -> None:
+        for i in range(len(a)):
+            a[i] = k
+
+
+    def vdiv(a: array.array[float], k: float) -> None:
+        for i in range(len(a)):
+            a[i] = a[i] / k
+''')
+
+
+@pytest.fixture(scope="module")
+def vec_built(tmp_path_factory):
+    d = tmp_path_factory.mktemp("vecmod")
+    path, proved, diags, c = build_source(d, VEC_SOURCE, "vecmod")
+    return Built(proved, diags, c, load_so(cbuild.build(c, path.stem, d / "out"))), d
+
+
+def _body(c, name):
+    start = c.index(f"static int tp_impl_{name}(", c.index("static const char tp_source"))
+    while not c[c.index(")", start) + 1:].startswith("\n{"):
+        start = c.index(f"static int tp_impl_{name}(", start + 1)
+    return c[start:c.index("\n}\n", start)]
+
+
+def test_vector_helper_only_where_the_architecture_guarantees_it():
+    h = (cbuild.RUNTIME_DIR / "tp_runtime.h").read_text()
+    guard = h[h.index("#define TP_BUF_VEC2") - 120:h.index("#define TP_BUF_VEC2")]
+    assert "__aarch64__" in guard and "x86" not in guard and "__SSE" not in guard
+    ld = h[h.index("tp_buf_ld2_f64("):h.index("tp_buf_st2_f64(")]
+    st = h[h.index("tp_buf_st2_f64("):h.index("#endif", h.index("tp_buf_st2_f64("))]
+    assert '"ldr %q0, %1"' in ld and '"str %q1, %0"' in st and '"Q"' in ld and '"=Q"' in st
+    assert '"memory"' not in ld + st
+    assert "vld1q" not in ld + st and "vst1q" not in ld + st      # no plain C access the compiler may split
+    plain = h.index("#ifdef TP_BUFFER_PLAIN")
+    assert h.index("#else", plain) < h.index("#define TP_BUF_VEC2")     # absent from the plain build
+
+
+def test_vectorised_loops_are_emitted_for_the_pattern_only(vec_built):
+    b, _ = vec_built
+    for n in ("vscale", "vmix", "vfill"):
+        body = _body(b.c, n)
+        assert "#ifdef TP_BUF_VEC2" in body and "tp_buf_st2_f64(" in body, n
+        assert "tp_buf_st_f64(" in body, n          # the scalar tail stays
+    for n in ("vdiv", "vtwo"):      # division can raise; b[i] is not proven in range: scalar only
+        assert "tp_buf_st2_f64" not in _body(b.c, n), n
+
+
+@pytest.mark.parametrize("n", [0, 1, 2, 3, 7, 1000, 4097])
+def test_vectorised_results_equal_the_interpreter(vec_built, n):
+    m = vec_built[0].module
+    mk = lambda: array.array("d", [i * 0.37 - 5 for i in range(n)])  # noqa: E731
+    differential(m, "vscale", lambda: (mk(), 1.7))
+    differential(m, "vmix", lambda: (mk(), 0.5))
+    differential(m, "vfill", lambda: (mk(), 3.25))
+    differential(m, "vdiv", lambda: (mk(), 4.0))
+    differential(m, "vtwo", lambda: (mk(), array.array("d", [i * 0.11 for i in range(n)])))
+
+
+def test_fallback_build_without_the_vector_helper_agrees(vec_built, tmp_path):
+    b, _ = vec_built
+    so = cbuild.build(b.c, "vecmod", tmp_path / "novec",
+                      flags=list(cbuild.DEFAULT_FLAGS) + ["-DTP_NO_BUF_VEC"])
+    m = load_so(so)
+    a = array.array("d", [i * 1.5 for i in range(33)])
+    m.vscale(a, 2.0)
+    assert a.tolist() == [i * 3.0 for i in range(33)]
+
+
+def test_concurrent_reader_sees_no_torn_element(vec_built):
+    import struct
+    import threading
+    m = vec_built[0].module
+    x = struct.unpack("<d", struct.pack("<Q", 0x1111111122222222))[0]
+    y = struct.unpack("<d", struct.pack("<Q", 0x3333333344444444))[0]
+    n = 4096
+    a = array.array("d", [x] * n)
+    ok = {struct.pack("<d", x), struct.pack("<d", y)}
+    stop = threading.Event()
+
+    def writer():
+        k = 0
+        while not stop.is_set():
+            m.vfill(a, y if k & 1 else x)
+            k += 1
+
+    t = threading.Thread(target=writer)
+    t.start()
+    bad = 0
+    try:
+        for _ in range(300):
+            raw = a.tobytes()
+            for off in range(0, len(raw), 8):
+                if raw[off:off + 8] not in ok:
+                    bad += 1
+    finally:
+        stop.set()
+        t.join()
+    assert bad == 0
+
+
+def test_no_vec_build_does_not_use_the_helper(vec_built, tmp_path):
+    import subprocess
+    c = tmp_path / "v.c"
+    c.write_text(vec_built[0].c)
+    import sysconfig
+    base = ["cc", "-E", "-P", "-I", str(cbuild.RUNTIME_DIR), "-I", sysconfig.get_paths()["include"]]
+    def body(*flags):
+        out = subprocess.run(base + list(flags) + [str(c)], capture_output=True, text=True, check=True).stdout
+        return _body(out, "vscale")
+    import platform
+    if platform.machine() != "arm64" and platform.machine() != "aarch64":
+        pytest.skip("the vector helper exists on aarch64 only")
+    assert "tp_buf_st2_f64(" in body()
+    assert "tp_buf_st2_f64(" not in body("-DTP_NO_BUF_VEC")