fix(deltanet): align decay formula with llama.cpp (#95)

Changed DeltaNet recurrent update from:
  OLD: S=decay*S; sk=S@K; d=beta*(V-sk); S=S+K*d
  NEW: sk=S@K; d=(V-decay*sk)*beta; S=decay*S+K*d

Matches llama.cpp gated_delta_net.cu:
  sk computed on ORIGINAL state before decay
  decay applied to sk in delta AND to S in update

Impact: fixes numerical difference with reference, but does NOT fix
the short-prompt instability alone. Document QA still works (2/2).
Short prompts still fail (0/5). The issue requires additional fixes
(L2 norm removal or Q scaling timing change per reference analysis).

Refs #95

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: unamedkr

quant.h

@@ -13922,50 +13922,56 @@ static void deltanet_forward(tq_model_t* model, tq_state_t* s, int l) {
         float decay = decay_vals[h]; /* precomputed exp(gate) */
 
 #ifdef __ARM_NEON
-        /* NEON-optimized: fused decay + sk computation.
-         * For each row i of state: decay state, accumulate sk.
-         * sk[j] = sum_i(S[i,j] * K[i]) after decay */
+        /* NEON-optimized: llama.cpp-aligned delta rule.
+         * Formula (matches gated_delta_net.cu):
+         *   sk = S @ K           (BEFORE decay)
+         *   d  = (V - g*sk) * beta
+         *   S  = g*S + K * d
+         *   o  = S @ Q
+         * The key difference from the previous impl: sk is computed
+         * on the ORIGINAL state, then decay is applied to both sk
+         * (in the delta) and S (in the update). This prevents
+         * short-prompt instability where early tokens have near-zero
+         * state and the decay-first approach loses information. */
         float* sk = s->delta_sk;
         memset(sk, 0, (size_t)dv * sizeof(float));
 
-        float32x4_t vdecay = vdupq_n_f32(decay);
+        /* Step A: sk = S @ K (on original state, BEFORE decay) */
         for (int i = 0; i < dk; i++) {
             float* sp = sh + i * dv;
             float ki = kh[i];
             float32x4_t vki = vdupq_n_f32(ki);
             int j = 0;
             for (; j + 3 < dv; j += 4) {
                 float32x4_t vs = vld1q_f32(sp + j);
-                vs = vmulq_f32(vs, vdecay);  /* decay */
-                vst1q_f32(sp + j, vs);        /* store decayed state */
                 float32x4_t vsk = vld1q_f32(sk + j);
-                vsk = vfmaq_f32(vsk, vs, vki); /* accumulate sk */
+                vsk = vfmaq_f32(vsk, vs, vki);
                 vst1q_f32(sk + j, vsk);
             }
             for (; j < dv; j++) {
-                sp[j] *= decay;
                 sk[j] += sp[j] * ki;
             }
         }
 
-        /* Delta: d = beta * (V - sk) */
+        /* Step B: d = (V - g*sk) * beta */
         float* d_vec = s->delta_dvec;
         float32x4_t vbeta = vdupq_n_f32(beta_h);
+        float32x4_t vdecay = vdupq_n_f32(decay);
         {
             int j = 0;
             for (; j + 3 < dv; j += 4) {
                 float32x4_t vv = vld1q_f32(vh + j);
-                float32x4_t vs = vld1q_f32(sk + j);
-                float32x4_t vd = vmulq_f32(vbeta, vsubq_f32(vv, vs));
+                float32x4_t vsk = vld1q_f32(sk + j);
+                float32x4_t vd = vmulq_f32(vbeta, vsubq_f32(vv, vmulq_f32(vdecay, vsk)));
                 vst1q_f32(d_vec + j, vd);
             }
             for (; j < dv; j++) {
-                d_vec[j] = beta_h * (vh[j] - sk[j]);
+                d_vec[j] = beta_h * (vh[j] - decay * sk[j]);
             }
         }
 
-        /* State update: S[i][j] += K[i] * d[j] (rank-1 outer product)
-         * + Output: o[j] = sum_i(S[i,j] * Q[i]) (simultaneously) */
+        /* Step C: S = g*S + K*d (state update)
+         * + Output: o = S @ Q (simultaneously) */
         float* oh = s->delta_out + h * dv;
         memset(oh, 0, (size_t)dv * sizeof(float));
 
@@ -13978,26 +13984,24 @@ static void deltanet_forward(tq_model_t* model, tq_state_t* s, int l) {
             int j = 0;
             for (; j + 3 < dv; j += 4) {
                 float32x4_t vs = vld1q_f32(sp + j);
+                vs = vmulq_f32(vs, vdecay);           /* S = g*S */
                 float32x4_t vd = vld1q_f32(d_vec + j);
-                vs = vfmaq_f32(vs, vki, vd);  /* S += K[i] * d */
+                vs = vfmaq_f32(vs, vki, vd);           /* S += K[i] * d */
                 vst1q_f32(sp + j, vs);
                 float32x4_t vo = vld1q_f32(oh + j);
-                vo = vfmaq_f32(vo, vs, vqi);   /* o += S * Q[i] */
+                vo = vfmaq_f32(vo, vs, vqi);           /* o += S * Q[i] */
                 vst1q_f32(oh + j, vo);
             }
             for (; j < dv; j++) {
-                sp[j] += ki * d_vec[j];
+                sp[j] = decay * sp[j] + ki * d_vec[j];
                 oh[j] += sp[j] * qi;
             }
         }
 #else
-        /* Scalar fallback */
-        /* Decay: S = S * exp(gate) */
-        for (int i = 0; i < dk * dv; i++) {
-            sh[i] *= decay;
-        }
+        /* Scalar fallback — llama.cpp-aligned formula:
+         * sk = S @ K, d = (V - g*sk) * beta, S = g*S + K*d, o = S @ Q */
 
-        /* Compute sk */
+        /* Compute sk = S @ K (original state, before decay) */
         float* sk = s->delta_sk;
         for (int j = 0; j < dv; j++) {
             float sum = 0.0f;
@@ -14007,20 +14011,20 @@ static void deltanet_forward(tq_model_t* model, tq_state_t* s, int l) {
             sk[j] = sum;
         }
 
-        /* Delta */
+        /* Delta: d = (V - g*sk) * beta */
         float* d_vec = s->delta_dvec;
         for (int j = 0; j < dv; j++) {
-            d_vec[j] = beta_h * (vh[j] - sk[j]);
+            d_vec[j] = beta_h * (vh[j] - decay * sk[j]);
         }
 
-        /* State update */
+        /* State update: S = g*S + K*d */
         for (int i = 0; i < dk; i++) {
             for (int j = 0; j < dv; j++) {
-                sh[i * dv + j] += kh[i] * d_vec[j];
+                sh[i * dv + j] = decay * sh[i * dv + j] + kh[i] * d_vec[j];
             }
         }
 
-        /* Output */
+        /* Output: o = S @ Q */
         float* oh = s->delta_out + h * dv;
         for (int j = 0; j < dv; j++) {
             float sum = 0.0f;