bong-water-water-bong · Nov 11, 2026
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -207,6 +207,11 @@ target_include_directories(npu_decode_note_test PRIVATE app)
 add_test(NAME npu_decode_note COMMAND npu_decode_note_test)
 add_dependencies(onebit npu_decode_note_test)
 add_test(NAME bump_hrx_report COMMAND bash ${CMAKE_SOURCE_DIR}/tests/bump_hrx_report_test.sh ${CMAKE_SOURCE_DIR}/scripts/bump-hrx-report.sh)
+# engine#1931 (first slice): the Whisper front-end's four f32 host kernels -- conv1d with stride/pad, layer norm with the
+# biased variance, the exact erf-GELU, and a row softmax. Header-only and device-free, so it runs wherever CI builds.
+add_executable(audio_ops_test tests/audio_ops_test.cpp)
+target_include_directories(audio_ops_test PRIVATE npu/forward/include)
+add_test(NAME audio_ops COMMAND audio_ops_test)
 add_test(NAME metrics_serve COMMAND ${CMAKE_SOURCE_DIR}/tests/metrics_serve.sh $<TARGET_FILE:onebit>)
 # serve --print-backend: the llama-server Lemonade starts itself, one hop (docs/serve.md): no GPU.
 add_test(NAME print_backend COMMAND ${CMAKE_SOURCE_DIR}/tests/print_backend.sh $<TARGET_FILE:onebit>)

diff --git a/npu/forward/include/audio_ops.h b/npu/forward/include/audio_ops.h
@@ -0,0 +1,86 @@
+// Copyright 2026 bong-water-water-bong
+// SPDX-License-Identifier: Apache-2.0
+//
+// audio_ops.h -- the four f32 host kernels a Whisper encoder front-end is made of: a temporal conv with stride and zero
+// padding, layer norm, erf-GELU and a row softmax. Header-only, no dependencies.
+//
+// Why these and not a graph executor: the encoder's arithmetic that is easy to get subtly wrong is here -- the mel
+// front-end is two strided 3-tap convs ([1280,128,3] stride 1 then [1280,1280,3] stride 2, both pad 1), and the GELU is
+// the *exact* erf form, which numpy does not have (the reference builds it with A&S 7.1.26, |err| < 1.5e-7). Every number
+// in tests/audio_ops_test.cpp was computed in double from the same formulas the reference uses, so a change here that
+// moves a value is visible before it reaches the device.
+
+#pragma once
+
+#include <cmath>
+#include <cstddef>
+
+namespace onebit::audio {
+
+// (T + 2*pad - k) / stride + 1, floor -- the ONNX Conv default for a valid (non-SAME) convolution.
+inline int conv1d_out_len(int T, int k, int stride, int pad) {
+    const int span = T + 2 * pad - k;
+    return span < 0 ? 0 : span / stride + 1;
+}
+
+// x: (cin, T) row-major. w: (cout, cin, k) row-major -- the layout the ONNX initializer carries, so a caller reading
+// weights through the pack's own loader needs no transpose. b: (cout,) or null. out: (cout, tout) row-major.
+inline void conv1d(const float* x, int cin, int T, const float* w, const float* b, int cout, int k, int stride, int pad,
+                   float* out) {
+    const int tout = conv1d_out_len(T, k, stride, pad);
+    for (int co = 0; co < cout; ++co) {
+        for (int t = 0; t < tout; ++t) {
+            float acc = (b != nullptr) ? b[co] : 0.0f;
+            for (int ci = 0; ci < cin; ++ci) {
+                const float* xr = x + static_cast<std::size_t>(ci) * T;
+                const float* wr = w + (static_cast<std::size_t>(co) * cin + ci) * k;
+                for (int j = 0; j < k; ++j) {
+                    const int src = t * stride + j - pad;
+                    if (src >= 0 && src < T) acc += wr[j] * xr[src];
+                }
+            }
+            out[static_cast<std::size_t>(co) * tout + t] = acc;
+        }
+    }
+}
+
+// y = (x - mean) / sqrt(var + eps) * gamma + beta, biased variance (divide by n), one row.
+inline void layer_norm(const float* x, const float* gamma, const float* beta, int n, float eps, float* out) {
+    double sum = 0.0;
+    for (int i = 0; i < n; ++i) sum += x[i];
+    const double mean = sum / n;
+    double var = 0.0;
+    for (int i = 0; i < n; ++i) {
+        const double d = x[i] - mean;
+        var += d * d;
+    }
+    var /= n;
+    const double inv = 1.0 / std::sqrt(var + eps);
+    for (int i = 0; i < n; ++i) {
+        const float g = (gamma != nullptr) ? gamma[i] : 1.0f;
+        const float b = (beta != nullptr) ? beta[i] : 0.0f;
+        out[i] = static_cast<float>((x[i] - mean) * inv) * g + b;
+    }
+}
+
+// The exact form: 0.5*x*(1+erf(x/sqrt(2))). `approximate=none` in the export, so tanh is not a drop-in.
+inline float gelu_erf(float x) {
+    return 0.5f * x * (1.0f + std::erf(x * 0.70710678118654752f));
+}
+
+// In place, one row: numerically stable (subtract the row max) so a logit of 100 does not overflow.
+inline void softmax_row(float* x, int n) {
+    if (n <= 0) return;
+    float mx = x[0];
+    for (int i = 1; i < n; ++i)
+        if (x[i] > mx) mx = x[i];
+    double sum = 0.0;
+    for (int i = 0; i < n; ++i) {
+        x[i] = std::exp(x[i] - mx);
+        sum += x[i];
+    }
+    const double inv = 1.0 / sum;
+    for (int i = 0; i < n; ++i) x[i] = static_cast<float>(x[i] * inv);
+}
+
+}  // namespace onebit::audio
diff --git a/tests/audio_ops_test.cpp b/tests/audio_ops_test.cpp
@@ -0,0 +1,137 @@
+// Copyright 2026 bong-water-water-bong
+// SPDX-License-Identifier: Apache-2.0
+//
+// audio_ops_test.cpp -- known-answer checks for npu/forward/include/audio_ops.h. No container, no device, no pack: every
+// expected value was computed in double in Python from the same formulas the Whisper encoder reference uses
+// (tools/whisper_encoder_ref.py), so this binary runs wherever CI builds, which is the point.
+
+#include <cmath>
+#include <cstdio>
+#include <string>
+#include <vector>
+
+#include "audio_ops.h"
+
+namespace {
+
+int g_failures = 0;
+
+void check(bool cond, const std::string& what) {
+    if (!cond) {
+        std::printf("FAIL %s\n", what.c_str());
+        ++g_failures;
+    }
+}
+
+void expect_eq(int got, int want, const std::string& what) {
+    if (got != want) {
+        std::printf("FAIL %s: got %d want %d\n", what.c_str(), got, want);
+        ++g_failures;
+    }
+}
+
+void expect_close(float got, double want, double tol, const std::string& what) {
+    const double d = std::fabs(static_cast<double>(got) - want);
+    if (!(d <= tol)) {
+        std::printf("FAIL %s: got %.10g want %.10g (|d| %.3g > %.3g)\n", what.c_str(), static_cast<double>(got), want, d,
+                    tol);
+        ++g_failures;
+    }
+}
+
+// 1. The front-end conv, at the shapes the encoder actually uses: two strided 3-tap convs, pad 1 on both ends.
+void conv1d_arms() {
+    // cin 2, T 5, cout 3, k 3, stride 2, pad 1 -> tout 3
+    const std::vector<float> x = {0.5f, -1.0f, 0.25f, 2.0f, -0.75f,   // channel 0
+                                  1.5f, 0.0f, -0.5f, 1.0f, 0.25f};   // channel 1
+    std::vector<float> w(3 * 2 * 3);
+    for (int i = 0; i < 18; ++i) w[i] = static_cast<float>((i + 1) * 0.0625);
+    const std::vector<float> b = {0.125f, -0.25f, 0.5f};
+
+    expect_eq(onebit::audio::conv1d_out_len(5, 3, 2, 1), 3, "conv1d_out_len T5 k3 s2 p1");
+    expect_eq(onebit::audio::conv1d_out_len(3000, 3, 2, 1), 1500, "conv1d_out_len 3000 -> 1500");
+
+    std::vector<float> out(3 * 3);
+    onebit::audio::conv1d(x.data(), 2, 5, w.data(), b.data(), 3, 3, 2, 1, out.data());
+
+    // double, same formulas: [0.46875, 0.6875, 0.484375, 0.46875, 0.96875, 1.046875, 1.59375, 2.375, 2.734375]
+    const double want[9] = {0.46875, 0.6875, 0.484375, 0.46875, 0.96875, 1.046875, 1.59375, 2.375, 2.734375};
+    for (int i = 0; i < 9; ++i) expect_close(out[i], want[i], 1e-6, "conv1d value " + std::to_string(i));
+
+    // edge case: stride 1, pad 0, k 3 -- the second conv's arithmetic without the padding, and a 1-tap full-window
+    // where tout * stride + k - pad == T exactly (no tail row may be emitted past the input)
+    std::vector<float> out2(2);
+    onebit::audio::conv1d(x.data(), 1, 5, w.data(), nullptr, 1, 3, 1, 0, out2.data());
+    expect_close(out2[0], 0.03076171875 * 0.0 + 0.0625 * 0.5 + 0.125 * -1.0 + 0.1875 * 0.25, 1e-6, "conv1d no bias");
+    expect_eq(onebit::audio::conv1d_out_len(5, 3, 1, 0), 3, "conv1d_out_len T5 k3 s1 p0");
+}
+
+// 2. Layer norm: biased variance, eps added under the root, gamma/beta applied after the normalise.
+void layernorm_arm() {
+    const std::vector<float> v = {0.25f, -1.5f, 2.0f, 0.75f};
+    const std::vector<float> g = {1.0f, 0.5f, -1.0f, 2.0f};
+    const std::vector<float> be = {0.0f, 1.0f, 0.25f, -0.5f};
+    std::vector<float> out(4);
+    onebit::audio::layer_norm(v.data(), g.data(), be.data(), 4, 1e-5f, out.data());
+    // double: [-0.0995034038, 0.2537244718, -1.0435442489, 0.0970204226]
+    const double want[4] = {-0.0995034038, 0.2537244718, -1.0435442489, 0.0970204226};
+    for (int i = 0; i < 4; ++i) expect_close(out[i], want[i], 1e-6, "layer_norm value " + std::to_string(i));
+
+    // A constant row normalises to exactly beta whatever eps is: the variance is 0, so the numerator is 0 too.
+    const std::vector<float> flat = {2.0f, 2.0f, 2.0f, 2.0f};
+    onebit::audio::layer_norm(flat.data(), nullptr, nullptr, 4, 1e-5f, out.data());
+    for (int i = 0; i < 4; ++i) expect_close(out[i], 0.0, 1e-7, "layer_norm flat " + std::to_string(i));
+    // ... and with no gamma/beta the default is identity scale, zero shift.
+    std::vector<float> out2(4);
+    const std::vector<float> v2 = {1.0f, 2.0f, 3.0f, 4.0f};
+    onebit::audio::layer_norm(v2.data(), nullptr, nullptr, 4, 1e-5f, out2.data());
+    expect_close(out2[0], -1.34163542, 1e-6, "layer_norm null gamma/beta");
+}
+
+// 3. erf-GELU -- the export is approximate=none, so tanh must not pass this.
+void gelu_arm() {
+    // double: -0.0455002639, 0.3457312306, 2.9959503059
+    expect_close(onebit::audio::gelu_erf(-2.0f), -0.0455002639, 1e-6, "gelu(-2)");
+    expect_close(onebit::audio::gelu_erf(0.5f), 0.3457312306, 1e-6, "gelu(0.5)");
+    expect_close(onebit::audio::gelu_erf(3.0f), 2.9959503059, 1e-6, "gelu(3)");
+
+    // The two properties that separate erf-GELU from tanh-GELU: GELU(x) - GELU(-x) == x exactly, and the tanh form is
+    // 1.7e-4 away at x = 2 -- far outside the 1e-6 this test runs at, so a tanh substitute cannot pass silently.
+    expect_close(onebit::audio::gelu_erf(2.0f) - onebit::audio::gelu_erf(-2.0f), 2.0, 1e-6, "gelu odd part");
+    const float tanh_form = 0.5f * 2.0f * (1.0f + std::tanh(0.7978845608f * (2.0f + 0.044715f * 8.0f)));
+    check(std::fabs(tanh_form - onebit::audio::gelu_erf(2.0f)) > 1e-5, "tanh GELU is distinguishable from erf GELU");
+}
+
+// 4. Softmax: normalised, row-max stable, and order-preserving.
+void softmax_arm() {
+    std::vector<float> s = {1.0f, -2.0f, 0.5f, 0.25f};
+    onebit::audio::softmax_row(s.data(), 4);
+    // double: [0.4697737513, 0.0233886579, 0.2849321833, 0.2219054075]
+    const double want[4] = {0.4697737513, 0.0233886579, 0.2849321833, 0.2219054075};
+    for (int i = 0; i < 4; ++i) expect_close(s[i], want[i], 1e-6, "softmax value " + std::to_string(i));
+    double sum = 0.0;
+    for (float v : s) sum += v;
+    expect_close(static_cast<float>(sum), 1.0, 1e-6, "softmax sums to 1");
+
+    // +1000 on every logit must not change the result: the encoder's attention logits reach this range after the
+    // 1/sqrt(d) scale is dropped on a long window.
+    std::vector<float> big = {1001.0f, 998.0f, 1000.5f, 1000.25f};
+    onebit::audio::softmax_row(big.data(), 4);
+    for (int i = 0; i < 4; ++i) expect_close(big[i], want[i], 1e-6, "softmax shift-invariance " + std::to_string(i));
+}
+
+}  // namespace
+
+int main() {
+    conv1d_arms();
+    layernorm_arm();
+    gelu_arm();
+    softmax_arm();
+    if (g_failures != 0) {
+        std::printf("audio ops: %d failures\n", g_failures);
+        return 1;
+    }
+    std::printf("audio ops: conv1d (strided + padded + edge), layer_norm (biased var, null affine), gelu (erf), softmax "
+                "(stable, shift-invariant) -- all pass\n");
+    return 0;
+}