npu: a launch's result is dequantised with one multiply per element, not one ldexp(): prefill 24.4-25.1 -> 29.8-30.6 tok/s on the served 27B pack (Refs #1990)

docs/npu.md

@@ -518,6 +518,23 @@ checks keep the contract honest without a device: `tests/npu_t2s_layout_test.cpp
 against the numbers the tool measured on a real image (gate/up: 50,135,040 B of pb16i → 23,674,880 B
 of stream), and `tests/q4nx_t2s_test.py` pins the code order byte for byte.
 
+**That single rounding costs one multiply per element, not one `ldexp()` (#1990).** A launch's result is
+`float32(S * 2^e_ref) * ascale` for every element, and `2^e_ref` depends only on the column. It was computed
+with a libm `ldexp()` call per element, after readback, where no launch phase column counted it: a 16-row
+prefill block on `bonsai2-27b-t2s2` is 416 launches of up to 16 x 17,408 elements, and that call was the
+block's largest host term. `q4nx_pb16i::dequant_rows()` computes the power once per column and multiplies.
+The value is the same, not close: a double times a normal power of two is exact while the product is a
+normal double, which holds for every exponent from -1022 to 960, and any column outside that range takes
+the per-element form as before. Measured in one `npu-hold` ticket, one binary, `NPU_DEQ_LDEXP=1` (the old
+form) against the default, 103 blocks per arm: the dequant term **137.4 -> 32.3 ms per block**, the block
+median **630.4 -> 515.0 ms** (25.4 -> 31.1 rows/s), and three 510-536-token prompts reached their first
+token at 24.4-25.1 tok/s before and **29.8-30.6 tok/s** after. That is the first pair of arms; the box
+slowed for both forms in the second pair (the old form measured 357 ms of dequant per block there), so the
+second pair is not quoted as a speedup. `NPU_DEQ_LDEXP=2` runs both forms and counts differing elements:
+0 of 568,483,840 on the device, and the greedy text is identical in all five arms.
+`tests/npu_pb16i_layout_test.cpp` holds the two forms to the same bits for one and two planes, the edges
+of the exponent range and exponents outside it. The `[bl]` line names the term as `deq`.
+
 **A t2s container wants a sealed set for the same reason a pb16i one does.** The layout is not what makes
 a phase change expensive: the **hw_context** is. Measured on this box in-loop on 2026-11-02, a t2s decode
 family shipped as eight per-role `_t2s2_m16` ELFs paid **+2.0-2.3 ms on every launch, in every role** — more

npu/forward/engine-src/npu_engine_i8ctx_inc.h

@@ -462,17 +462,35 @@ struct I8Ctx {
     // row is dequantised with its own activation scale.
     void dequant_pb16i(float* out, int am, int an, const float* ascales,
                        const std::vector<int16_t>& eref, int planes) {
+        // engine#1990: one 2^e_ref per column instead of one ldexp() per element; same bits (q4nx_pb16i.h).
+        // NPU_DEQ_LDEXP=1 keeps the per-element form for an A/B; =2 runs both and counts differing elements.
+        static const int legacy = [] {
+            const char* e = getenv("NPU_DEQ_LDEXP");
+            return e ? atoi(e) : 0;
+        }();
+        if (legacy != 1)
+            q4nx_pb16i::dequant_rows(out, Cm, ND, am, an, ascales, eref.data(), planes, deq_pow2_);
+        if (legacy == 0) return;
         const int row = planes > 1 ? 2 * ND : ND;
+        static long calls = 0, elems = 0, differing = 0;
         for (int m = 0; m < am; m++) {
             const int32_t* lo = Cm + (size_t)m * (size_t)row;
             const float asc = ascales[m];
             for (int n = 0; n < an; n++) {
                 const int64_t s = q4nx_pb16i::sum_of(lo + n, lo + ND + n, planes > 1 ? 1 : 0);
-                out[(size_t)m * (size_t)an + n] =
-                    q4nx_pb16i::scale_i64(s, eref[(size_t)n]) * asc;
+                const float v = q4nx_pb16i::scale_i64(s, eref[(size_t)n]) * asc;
+                float& o = out[(size_t)m * (size_t)an + n];
+                if (legacy == 2 && memcmp(&o, &v, sizeof v) != 0) differing++;
+                o = v;
             }
         }
+        if (legacy == 2) {
+            elems += (long)am * an;
+            if (++calls % 2000 == 0)
+                fprintf(stderr, "[deqcheck] calls=%ld elements=%ld differing=%ld\n", calls, elems, differing);
+        }
     }
+    std::vector<double> deq_pow2_;   // dequant_pb16i's per-column scratch
 
     // Cache-friendly per-column int8 quantization into a [KD][ND] padded buffer.
     // The original column-at-a-time loop reads w[i*N + j] with a stride of N

npu/forward/include/q4nx_lh.h

@@ -55,6 +55,7 @@ enum LhSlot {
     LH_RQ_RESP,         // response serialise + write
     LH_PACK_UPLOAD,     // weight BO create + pack + upload (wcache_ miss only)
     LH_HEAD_ROWS,       // run_lm_head_rows(): the head for a batched verify pass
+    LH_DEQ_OUT,         // a layer launch's result -> f32 (after readback, so in no launch phase column)
     LH_N
 };
 inline bool lh_on() {
@@ -96,7 +97,7 @@ inline const char* lh_name(int s) {
         "decode (step loop)", "  of-which: sample+grammar",
         "  of-which: detokenise token", "response write",
         "  of-which: weight pack+upload (wcache miss)",
-        "lm head (batch rows)"};
+        "lm head (batch rows)", "launch output dequant"};
     return (s >= 0 && s < LH_N) ? kn[s] : "?";
 }
 }  // namespace q4nx_lh

npu/forward/include/q4nx_pb16i.h

@@ -152,6 +152,49 @@ inline int64_t sum_of(const int32_t* lo, const int32_t* hi, int mode) {
     return ((int64_t)*hi << 21) + (int64_t)*lo;
 }
 
+// A launch's whole result, dequantised: out[m][n] = scale_i64(S[m][n], e_ref[n]) * ascales[m], bit for bit.
+//
+// scale_i64() is one libm call per element, and a prefill block is 416 launches of up to 16 x 17408 elements, so
+// that call was a block's largest host term outside the device wait (engine#1990: 0.15-0.50 ms per 16-row launch,
+// measured on the host). 2^e_ref depends only on the column, so it is computed once per column here and each element
+// is a multiply. That is the same value, not an approximation: for a normal power of two p, (double)s * p is exact
+// whenever the product is a normal double, which holds for every nonzero s when e_ref is in [kExpLo, kExpHi]
+// (|s| < 2^63, so the product stays under 2^1023, and it is at least 2^-1022). ldexp() is exact under the same
+// condition, so both forms hand the same double to the one float rounding. Outside that range -- no container
+// has such a column, but the header does not get to assume it -- every element takes scale_i64() as before.
+//
+// `Cm` is the result as the kernel wrote it: per row, one int32 plane of ND values, or [lo | hi] when planes > 1.
+// `pow2` is caller-owned scratch so a launch does not allocate.
+constexpr int kExpLo = -1022, kExpHi = 960;
+inline void dequant_rows(float* out, const int32_t* Cm, int ND, int am, int an, const float* ascales,
+                         const int16_t* eref, int planes, std::vector<double>& pow2) {
+    const int row = planes > 1 ? 2 * ND : ND;
+    const int mode = planes > 1 ? 1 : 0;
+    bool table = true;
+    pow2.resize((size_t)(an > 0 ? an : 0));
+    for (int n = 0; n < an; n++) {
+        const int e = (int)eref[n];
+        if (e < kExpLo || e > kExpHi) { table = false; break; }
+        pow2[(size_t)n] = std::ldexp(1.0, e);
+    }
+    for (int m = 0; m < am; m++) {
+        const int32_t* lo = Cm + (size_t)m * (size_t)row;
+        const int32_t* hi = lo + ND;
+        const float asc = ascales[m];
+        float* o = out + (size_t)m * (size_t)an;
+        if (!table) {
+            for (int n = 0; n < an; n++) o[n] = scale_i64(sum_of(lo + n, hi + n, mode), eref[n]) * asc;
+        } else if (mode) {
+            const double* p = pow2.data();
+            for (int n = 0; n < an; n++)
+                o[n] = (float)((double)(((int64_t)hi[n] << 21) + (int64_t)lo[n]) * p[n]) * asc;
+        } else {
+            const double* p = pow2.data();
+            for (int n = 0; n < an; n++) o[n] = (float)((double)lo[n] * p[n]) * asc;
+        }
+    }
+}
+
 // The sidecar files a container's pb16i set is keyed on: blk_<layer>_<tag>_weight.
 inline std::string sidecar_stem(int layer, const char* tag) {
     return std::string("blk_") + std::to_string(layer) + "_" + tag + "_weight";

npu/forward/src/q4nx_forward.cpp

@@ -2483,10 +2483,10 @@ struct BatchPhases {
         const double row = M > 0 ? (double)M : 1.0;
         std::fprintf(stderr,
                      "[bl] M=%2d gdn %7.1f ab %7.1f core %7.1f attn %6.1f mlp %6.1f rms %5.1f "
-                     "res %5.1f gemm %7.1f head %6.1f | per row: gdn %6.1f attn %5.1f "
+                     "res %5.1f gemm %7.1f head %6.1f deq %6.1f | per row: gdn %6.1f attn %5.1f "
                      "gemm %6.1f head %5.1f\n",
                      M, ms(LH_GDN), ms(LH_GDN_AB), ms(LH_GDN_CORE), ms(LH_ATTN), ms(LH_MLP),
-                     ms(LH_RMS), ms(LH_RES), gemm, ms(LH_HEAD_ROWS),
+                     ms(LH_RMS), ms(LH_RES), gemm, ms(LH_HEAD_ROWS), ms(LH_DEQ_OUT),
                      ms(LH_GDN) / row, ms(LH_ATTN) / row, gemm / row, ms(LH_HEAD_ROWS) / row);
     }
 };

npu/target/xdna/xdna_target.cpp

@@ -1080,8 +1080,13 @@ void XdnaTarget::npu_gemm_rows(const GemvArgs& ar, std::vector<float>& out, std:
         const auto tw = prof ? std::chrono::steady_clock::now() : std::chrono::steady_clock::time_point{};
         kctx.readback();
         const auto tr = prof ? std::chrono::steady_clock::now() : std::chrono::steady_clock::time_point{};
-        out.assign((size_t)M * (size_t)N, 0.0f);
-        kctx.dequant_pb16i(out.data(), M, N, ascales.data(), slot->eref, kctx.pb16i_planes);
+        {
+            // `out` is already M x N (sized before the refusal paths above) and every element is written here, so it
+            // is not zero-filled a second time. Timed on its own: this runs after `tr`, so no launch phase column
+            // ever counted it, and the batched [lh] table had it under UNATTRIBUTED (engine#1990).
+            q4nx_lh::LhScope _t_deq(LH_DEQ_OUT);
+            kctx.dequant_pb16i(out.data(), M, N, ascales.data(), slot->eref, kctx.pb16i_planes);
+        }
         // The batched (prefill) twin of the trace above: a short prompt runs step_batch, so this
         // is the path a 1-token request actually takes.
         if (const char* tl = getenv("NPU_PB16I_TRACE")) {

tests/npu_pb16i_layout_test.cpp

@@ -206,6 +206,56 @@ int main() {
     for (const geom& g : bad)
         check(!q4nx_pb16i::fits(g.K, g.N), "a non-whole-tile geometry must be refused", g.K, g.N);
 
+    // 7: dequant_rows() (engine#1990) is scale_i64() per element, bit for bit -- one plane and two, the exponents a
+    // container uses, the edges of the table's range, and exponents outside it (which must take the per-element form
+    // and still agree). The sums include 0, +-1, the int32 extremes and two-plane values past 2^53, where
+    // (double)s itself rounds: both forms round it the same way first.
+    {
+        struct Case { int planes, M, N, e_lo, e_hi; const char* what; };
+        const Case cases[] = {
+            {1, 16, 1024, -40, -8, "one plane, container-like exponents"},
+            {2, 16, 1024, -40, -8, "two planes, container-like exponents"},
+            {1, 3, 257, q4nx_pb16i::kExpLo, q4nx_pb16i::kExpLo + 3, "one plane, lowest table exponents"},
+            {2, 3, 257, q4nx_pb16i::kExpHi - 3, q4nx_pb16i::kExpHi, "two planes, highest table exponents"},
+            {2, 2, 64, -1100, -1000, "exponents below the table range (per-element path)"},
+            {1, 2, 64, 950, 1100, "exponents above the table range (per-element path)"},
+            {1, 1, 1, 0, 0, "a single element"},
+        };
+        uint32_t x = 0x9e3779b9u;
+        auto next = [&]() { x ^= x << 13; x ^= x >> 17; x ^= x << 5; return x; };
+        size_t compared = 0, nonzero = 0;
+        for (const Case& c : cases) {
+            const int ND = c.N + 5;                      // the kernel's row is wider than the columns read
+            const int row = c.planes > 1 ? 2 * ND : ND;
+            std::vector<int32_t> Cm((size_t)c.M * (size_t)row);
+            for (int32_t& v : Cm) v = (int32_t)next();
+            const int32_t edge[] = {0, 1, -1, INT32_MAX, INT32_MIN};
+            for (size_t i = 0; i < Cm.size() && i < 40; i++) Cm[i] = edge[i % 5];
+            std::vector<int16_t> eref((size_t)c.N);
+            for (int n = 0; n < c.N; n++) eref[(size_t)n] = (int16_t)(c.e_lo + (int)(next() % (uint32_t)(c.e_hi - c.e_lo + 1)));
+            std::vector<float> asc((size_t)c.M);
+            for (float& a : asc) a = 1e-4f * (float)(1 + next() % 5000);
+            std::vector<float> want((size_t)c.M * (size_t)c.N), got(want.size(), -7.0f);
+            for (int m = 0; m < c.M; m++) {
+                const int32_t* lo = Cm.data() + (size_t)m * (size_t)row;
+                for (int n = 0; n < c.N; n++)
+                    want[(size_t)m * (size_t)c.N + n] =
+                        q4nx_pb16i::scale_i64(q4nx_pb16i::sum_of(lo + n, lo + ND + n, c.planes > 1 ? 1 : 0),
+                                              eref[(size_t)n]) * asc[(size_t)m];
+            }
+            std::vector<double> scratch;
+            q4nx_pb16i::dequant_rows(got.data(), Cm.data(), ND, c.M, c.N, asc.data(), eref.data(), c.planes, scratch);
+            size_t diff = 0;
+            for (size_t i = 0; i < want.size(); i++) {
+                if (std::memcmp(&want[i], &got[i], sizeof(float)) != 0) diff++;
+                if (want[i] != 0.0f) nonzero++;
+            }
+            compared += want.size();
+            if (diff) { std::printf("FAIL dequant_rows != scale_i64 per element: %s (%zu of %zu differ)\n", c.what, diff, want.size()); failures++; }
+        }
+        if (nonzero < compared / 4) { std::printf("FAIL dequant_rows arm compares mostly zeros (%zu of %zu nonzero)\n", nonzero, compared); failures++; }
+    }
+
     if (failures) {
         std::printf("npu_pb16i_layout_test: %d failure(s)\n", failures);
         return 1;