SSM NEON vectorization + threading, Q6_K drift fix, debug cleanup

Extract SSM forward pass into backend range functions (ssm_neon.c,
ssm_scalar.c) with threadpool dispatch over V-heads. NEON kernels
use 4x-unrolled vmlaq_f32 for delta recurrence state operations,
vectorized L2 norm and RMSNorm+SiLU gate.

Fix Q6_K SDOT accumulation drift by factoring block scale `d` out
of per-sub-block additions — single `d * block_sum` per 256-element
block instead of 8 separate `d * ...` additions.

Remove all 21 #ifdef DEBUG blocks from transformer.c (-332 lines),
delete 5 verify_*.py scripts, update .gitignore for stale artifacts.

Author: farkasmark

.gitignore

@@ -2,13 +2,25 @@
 bitnet
 src/*.o
 src/quant/*.o
+src/transformer/*.o
+*.dSYM/
 test_gguf
 test_quant
 test_tokenizer
 test_transformer
 test_e2e
 test_safety
+test_threadpool
+test_arena
+test_q2k
+test_prefill
+test_kv_f16
 dump_model
+bench_kernels
+bench/*.wasm
+bench/*.js
+bitnet_prof*
+__pycache__/
 
 # WASM build output
 wasm/bitnet.js

Makefile

@@ -46,7 +46,8 @@ ifneq ($(filter arm% aarch%,$(UNAME_M)),)
 
   TRANSFORMER_BACKEND = src/transformer/rmsnorm_neon.c src/transformer/rmsnorm_scalar.c \
     src/transformer/gqa_neon.c src/transformer/gqa_scalar.c \
-    src/transformer/logits_neon.c src/transformer/logits_scalar.c
+    src/transformer/logits_neon.c src/transformer/logits_scalar.c \
+    src/transformer/ssm_neon.c src/transformer/ssm_scalar.c
 else
   # x86: AVX2 + scalar
   QUANT_BACKEND = src/quant/x_quant_avx2.c \
@@ -73,7 +74,8 @@ else
 
   TRANSFORMER_BACKEND = src/transformer/rmsnorm_avx2.c src/transformer/rmsnorm_scalar.c \
     src/transformer/gqa_avx2.c src/transformer/gqa_scalar.c \
-    src/transformer/logits_avx2.c src/transformer/logits_scalar.c
+    src/transformer/logits_avx2.c src/transformer/logits_scalar.c \
+    src/transformer/ssm_scalar.c
 endif
 
 QUANT_SRCS = $(QUANT_COMMON) $(QUANT_BACKEND)
@@ -204,7 +206,8 @@ AVX2_QUANT_SRCS = $(QUANT_COMMON) \
 
 AVX2_TRANSFORMER_BACKEND = src/transformer/rmsnorm_avx2.c src/transformer/rmsnorm_scalar.c \
     src/transformer/gqa_avx2.c src/transformer/gqa_scalar.c \
-    src/transformer/logits_avx2.c src/transformer/logits_scalar.c
+    src/transformer/logits_avx2.c src/transformer/logits_scalar.c \
+    src/transformer/ssm_scalar.c
 
 AVX2_SRCS = src/platform.c src/gguf.c $(AVX2_QUANT_SRCS) src/model.c \
             src/transformer.c $(AVX2_TRANSFORMER_BACKEND) src/tokenizer.c src/sampler.c \

include/transformer_internal.h

@@ -90,4 +90,48 @@ void bn_transformer_logits_f16_wasm_range(void *ctx, int start, int end);
 void bn_transformer_logits_f16_scalar_range(void *ctx, int start, int end);
 void bn_transformer_logits_f32_range(void *ctx, int start, int end);
 
+// --- SSM context structs ---
+
+typedef struct {
+    float *qkv;            // [qkv_dim] input/output
+    float *conv_state;     // [(kern-1) * qkv_dim]
+    const float *conv1d_w; // [qkv_dim * kern]
+    int qkv_dim, kern;
+} BnSSMConvCtx;
+
+typedef struct {
+    float *q, *k;          // [key_dim] each
+    int head_dim;
+} BnSSML2NormCtx;
+
+typedef struct {
+    float *state, *out;
+    const float *q, *k;
+    float *v;              // also temp for sk
+    const float *alpha, *beta;
+    int num_k_heads, head_k_dim, head_v_dim;
+    float q_scale;
+} BnSSMDeltaCtx;
+
+typedef struct {
+    float *out;
+    const float *z, *norm_w;
+    float eps;
+    int head_v_dim;
+} BnSSMGateCtx;
+
+// --- SSM range function declarations ---
+
+void bn_transformer_ssm_conv_silu_neon_range(void *ctx, int start, int end);
+void bn_transformer_ssm_conv_silu_scalar_range(void *ctx, int start, int end);
+
+void bn_transformer_ssm_l2norm_neon_range(void *ctx, int start, int end);
+void bn_transformer_ssm_l2norm_scalar_range(void *ctx, int start, int end);
+
+void bn_transformer_ssm_delta_neon_range(void *ctx, int start, int end);
+void bn_transformer_ssm_delta_scalar_range(void *ctx, int start, int end);
+
+void bn_transformer_ssm_gate_neon_range(void *ctx, int start, int end);
+void bn_transformer_ssm_gate_scalar_range(void *ctx, int start, int end);
+
 #endif // BN_TRANSFORMER_INTERNAL_H

src/quant/q6k_neon_sdot.c

@@ -28,6 +28,8 @@ void bn_quant_q6k_neon_sdot_range(void *ctx, int row_start, int row_end) {
             const int8_t *xb = x_q + (b * BN_QK_K);
             const float *xs = x_scales + (b * 8);
 
+            // Accumulate per-block to reduce float rounding from repeated d*
+            float block_sum = 0.0f;
             for (int chunk = 0; chunk < 2; chunk++) {
                 uint8x16_t ql0 = vld1q_u8(ql);
                 uint8x16_t ql1 = vld1q_u8(ql + 16);
@@ -36,7 +38,7 @@ void bn_quant_q6k_neon_sdot_range(void *ctx, int row_start, int row_end) {
                 uint8x16_t qh0 = vld1q_u8(qh);
                 uint8x16_t qh1 = vld1q_u8(qh + 16);
 
-                // Unpack 8 weight vectors (identical to q6k_neon.c)
+                // Unpack 8 weight vectors
                 int8x16_t w0a = vsubq_s8(vreinterpretq_s8_u8(vorrq_u8(
                     vandq_u8(ql0, mask_lo4),
                     vshlq_n_u8(vandq_u8(qh0, mask_2), 4))), bias32);
@@ -62,40 +64,37 @@ void bn_quant_q6k_neon_sdot_range(void *ctx, int row_start, int row_end) {
                     vshrq_n_u8(ql3, 4),
                     vshlq_n_u8(vshrq_n_u8(qh1, 6), 4))), bias32);
 
-                // 4 pairs, each pair = 32 elements = 2 weight sub-blocks sharing 1 activation scale
-                // Pair 0: w0a(sc[0]) + w0b(sc[1]), dx = xs[chunk*4 + 0]
+                // 4 pairs: d factored out, accumulated into block_sum
                 float dx0 = xs[chunk * 4 + 0];
                 int32x4_t s0a = vdotq_s32(zero, w0a, vld1q_s8(xb));
                 int32x4_t s0b = vdotq_s32(zero, w0b, vld1q_s8(xb + 16));
-                row_sum += d * dx0 * ((float)sc[0] * (float)vaddvq_s32(s0a)
-                                    + (float)sc[1] * (float)vaddvq_s32(s0b));
+                block_sum += dx0 * ((float)sc[0] * (float)vaddvq_s32(s0a)
+                                  + (float)sc[1] * (float)vaddvq_s32(s0b));
 
-                // Pair 1: w1a(sc[2]) + w1b(sc[3]), dx = xs[chunk*4 + 1]
                 float dx1 = xs[chunk * 4 + 1];
                 int32x4_t s1a = vdotq_s32(zero, w1a, vld1q_s8(xb + 32));
                 int32x4_t s1b = vdotq_s32(zero, w1b, vld1q_s8(xb + 48));
-                row_sum += d * dx1 * ((float)sc[2] * (float)vaddvq_s32(s1a)
-                                    + (float)sc[3] * (float)vaddvq_s32(s1b));
+                block_sum += dx1 * ((float)sc[2] * (float)vaddvq_s32(s1a)
+                                  + (float)sc[3] * (float)vaddvq_s32(s1b));
 
-                // Pair 2: w2a(sc[4]) + w2b(sc[5]), dx = xs[chunk*4 + 2]
                 float dx2 = xs[chunk * 4 + 2];
                 int32x4_t s2a = vdotq_s32(zero, w2a, vld1q_s8(xb + 64));
                 int32x4_t s2b = vdotq_s32(zero, w2b, vld1q_s8(xb + 80));
-                row_sum += d * dx2 * ((float)sc[4] * (float)vaddvq_s32(s2a)
-                                    + (float)sc[5] * (float)vaddvq_s32(s2b));
+                block_sum += dx2 * ((float)sc[4] * (float)vaddvq_s32(s2a)
+                                  + (float)sc[5] * (float)vaddvq_s32(s2b));
 
-                // Pair 3: w3a(sc[6]) + w3b(sc[7]), dx = xs[chunk*4 + 3]
                 float dx3 = xs[chunk * 4 + 3];
                 int32x4_t s3a = vdotq_s32(zero, w3a, vld1q_s8(xb + 96));
                 int32x4_t s3b = vdotq_s32(zero, w3b, vld1q_s8(xb + 112));
-                row_sum += d * dx3 * ((float)sc[6] * (float)vaddvq_s32(s3a)
-                                    + (float)sc[7] * (float)vaddvq_s32(s3b));
+                block_sum += dx3 * ((float)sc[6] * (float)vaddvq_s32(s3a)
+                                  + (float)sc[7] * (float)vaddvq_s32(s3b));
 
                 xb += 128;
                 ql += 64;
                 qh += 32;
                 sc += 8;
             }
+            row_sum += d * block_sum;
         }
         c->out[row] = row_sum;
     }