fix: portable rotation sampler; exact reflection sign; pin output in CI

.github/workflows/ci.yml

@@ -27,5 +27,9 @@ jobs:
       - name: Build
         run: bazel build --disk_cache=~/.cache/bazel-disk //...
 
+      # determinism_test pins output across libstdc++ (Linux) and libc++ (macOS).
+      - name: Test
+        run: bazel test --disk_cache=~/.cache/bazel-disk --test_output=errors //...
+
       - name: Run example
         run: bazel run --disk_cache=~/.cache/bazel-disk //:turboquant_example

BUILD.bazel

@@ -3,6 +3,7 @@
 
 load("@rules_cc//cc:cc_binary.bzl", "cc_binary")
 load("@rules_cc//cc:cc_library.bzl", "cc_library")
+load("@rules_cc//cc:cc_test.bzl", "cc_test")
 
 cc_library(
     name = "turboquant",
@@ -38,3 +39,11 @@ cc_binary(
     srcs = ["benchmarks/benchmark.cpp"],
     deps = [":turboquant"],
 )
+
+# Pins exact output for fixed seeds across platforms and standard libraries.
+cc_test(
+    name = "determinism_test",
+    srcs = ["tests/determinism_test.cpp"],
+    deps = [":turboquant"],
+    size = "small",
+)

README.md

@@ -81,6 +81,12 @@ Single-threaded, Apple M4 Max, `bazel run -c opt //:turboquant_benchmark`. Compr
 
 Throughput is dominated by the dense d x d rotation matvec, so it scales as O(d^2) and is nearly independent of bitwidth (see Limitations). Run `bazel run -c opt //:turboquant_benchmark` for the full sweep (d = 256/768/1536, b = 1-4).
 
+## Determinism
+
+A quantizer built from the same `(dim, bitwidth, seed)` produces the same codes on every supported platform. The rotation and QJL projection are drawn with a portable Box-Muller sampler over `std::mt19937`'s integers rather than `std::normal_distribution`, whose algorithm the C++ standard leaves unspecified (libc++ and libstdc++ return different sequences from the same engine). `tests/determinism_test.cpp` pins the output with golden hashes and runs in CI on Linux and macOS.
+
+`turboquant::kAlgorithmVersion` changes whenever that output changes. If you persist codes, record it alongside `(dim, bitwidth, seed)` and refuse codes from another version.
+
 ## Limitations
 
 - Bitwidths 1-4 only (precomputed codebooks). Extending to higher bitwidths requires solving the Lloyd-Max optimization for the Beta distribution at the desired precision.

include/turboquant/types.h

@@ -10,6 +10,15 @@
 
 namespace turboquant {
 
+// Version of the numeric output: bumped whenever the rotation, projection or
+// codes produced for a given (dim, bitwidth, seed) change. Anything that
+// persists codes should record it, so codes from another version are refused.
+//   1: rotation and projection drawn with std::normal_distribution, whose
+//      output differs between standard libraries (libc++ vs libstdc++).
+//   2: portable Box-Muller sampler; reflection sign from the QR coefficients.
+//      Pinned by tests/determinism_test.cpp.
+inline constexpr int kAlgorithmVersion = 2;
+
 using Vec = Eigen::VectorXf;
 using Mat = Eigen::MatrixXf;
 

src/rotation.cpp

@@ -3,33 +3,93 @@
 
 #include "turboquant/rotation.h"
 
+#include <cmath>
+#include <cstdint>
+#include <numbers>
+
 namespace turboquant {
 
+namespace {
+
+// N(0, 1) samples defined only in terms of mt19937's integer output.
+//
+// std::normal_distribution is not portable: the standard leaves its algorithm
+// unspecified, and libc++ and libstdc++ return different sequences from the
+// same engine, so a rotation built from one seed differed between macOS and
+// Linux. This is Box-Muller in double over 53-bit uniforms. The integer
+// handling is exact; log, sqrt, sin and cos were measured bit-identical on
+// libc++/Apple clang and libstdc++/GCC (see tests/determinism_test.cpp).
+class PortableNormal {
+public:
+    explicit PortableNormal(std::mt19937& rng) : rng_(rng) {}
+
+    float operator()() {
+        if (has_spare_) {
+            has_spare_ = false;
+            return static_cast<float>(spare_);
+        }
+        const double u1 = open01();
+        const double u2 = open01();
+        const double r = std::sqrt(-2.0 * std::log(u1));
+        const double theta = 2.0 * std::numbers::pi * u2;
+        spare_ = r * std::sin(theta);
+        has_spare_ = true;
+        return static_cast<float>(r * std::cos(theta));
+    }
+
+private:
+    // Uniform on the open interval (0, 1) from two 32-bit draws; never 0 or 1.
+    double open01() {
+        const std::uint64_t a = rng_() >> 5;  // 27 bits
+        const std::uint64_t b = rng_() >> 6;  // 26 bits
+        const double k = static_cast<double>((a << 26) | b);  // < 2^53, exact
+        return (k + 0.5) / 9007199254740992.0;                // / 2^53, exact
+    }
+
+    std::mt19937& rng_;
+    bool has_spare_ = false;
+    double spare_ = 0.0;
+};
+
+// Q from a Householder QR is the product of one reflector per nonzero tau, so
+// det(Q) = (-1)^(number of nonzero taus): exact, O(dim). Q.determinant() cost
+// another O(dim^3) and, in float, collapsed to -0 from dim ~384 up, so the
+// sign fix below silently never ran at practical dimensions.
+bool is_reflection(const Eigen::HouseholderQR<Mat>& qr) {
+    int reflectors = 0;
+    for (Eigen::Index i = 0; i < qr.hCoeffs().size(); ++i)
+        if (qr.hCoeffs()(i) != 0.0f)
+            ++reflectors;
+    return (reflectors % 2) == 1;
+}
+
+}  // namespace
+
 Mat make_rotation_matrix(int dim, std::mt19937& rng) {
-    std::normal_distribution<float> normal(0.0f, 1.0f);
+    PortableNormal normal(rng);
 
     Mat G(dim, dim);
     for (int i = 0; i < dim; ++i)
         for (int j = 0; j < dim; ++j)
-            G(i, j) = normal(rng);
+            G(i, j) = normal();
 
     Eigen::HouseholderQR<Mat> qr(G);
     Mat Q = qr.householderQ() * Mat::Identity(dim, dim);
 
     // Ensure det = +1 (proper rotation, not reflection)
-    if (Q.determinant() < 0.0f)
+    if (is_reflection(qr))
         Q.col(0) *= -1.0f;
 
     return Q;
 }
 
 Mat make_gaussian_matrix(int rows, int cols, std::mt19937& rng) {
-    std::normal_distribution<float> normal(0.0f, 1.0f);
+    PortableNormal normal(rng);
 
     Mat S(rows, cols);
     for (int i = 0; i < rows; ++i)
         for (int j = 0; j < cols; ++j)
-            S(i, j) = normal(rng);
+            S(i, j) = normal();
 
     return S;
 }

tests/determinism_test.cpp

@@ -0,0 +1,107 @@
+// Copyright (c) 2026 Edge AI
+// SPDX-License-Identifier: MIT
+
+// Pins the exact output for fixed seeds, so a standard library, compiler or
+// platform that changes it fails here instead of silently producing codes that
+// another build — same seed, same parameters — cannot read.
+//
+// CI runs this on both libstdc++ (Linux) and libc++ (macOS). If a change to the
+// output is intentional, bump kAlgorithmVersion and regenerate the goldens:
+//   bazel run //:determinism_test -- --print
+
+#include <turboquant/turboquant.h>
+
+#include <cstdint>
+#include <cstdio>
+#include <cstring>
+#include <iterator>
+#include <random>
+#include <vector>
+
+namespace {
+
+using namespace turboquant;
+
+std::uint64_t fnv1a(const void* data, std::size_t n,
+                    std::uint64_t h = 1469598103934665603ull) {
+    const auto* p = static_cast<const unsigned char*>(data);
+    for (std::size_t i = 0; i < n; ++i) {
+        h ^= p[i];
+        h *= 1099511628211ull;
+    }
+    return h;
+}
+
+// Gaussian matrix bits: exact by construction, independent of Eigen's QR.
+std::uint64_t gaussianHash(int dim) {
+    std::mt19937 rng(42);
+    const Mat s = make_gaussian_matrix(dim, dim, rng);
+    return fnv1a(s.data(), static_cast<std::size_t>(s.size()) * sizeof(float));
+}
+
+// Codes (MSE indices + QJL signs) for 8 test vectors built from integers only.
+std::uint64_t codesHash(int dim) {
+    std::mt19937 rng(42);
+    const QuantizerProd qp(dim, 3, rng);
+    std::mt19937 vec_rng(7);
+    std::uint64_t h = 1469598103934665603ull;
+    for (int k = 0; k < 8; ++k) {
+        Vec x(dim);
+        for (int i = 0; i < dim; ++i)
+            x(i) = static_cast<float>(static_cast<std::int64_t>(vec_rng()) - 2147483648LL) /
+                   2147483648.0f;
+        x.normalize();
+        const QuantizedProd q = qp.quantize(x);
+        h = fnv1a(q.mse_part.indices.data(), q.mse_part.indices.size(), h);
+        h = fnv1a(q.qjl_signs.data(), q.qjl_signs.size(), h);
+    }
+    return h;
+}
+
+struct Golden {
+    int dim;
+    std::uint64_t gaussian;
+    std::uint64_t codes;
+};
+
+// Generated with --print, algorithm version 2. Verified identical on macOS arm64
+// (libc++, Apple clang 21) and Linux arm64 (libstdc++, g++ 11.4), each in
+// fastbuild, -c opt, and -c opt with -std=c++20.
+constexpr Golden kGoldens[] = {
+    {64, 0x4a335dee61b821bcull, 0xcfda2dd49e33c8bbull},
+    {256, 0x04e52ec727f12707ull, 0xae17ac1b26b87814ull},
+};
+
+}  // namespace
+
+int main(int argc, char** argv) {
+    const bool print = argc > 1 && std::strcmp(argv[1], "--print") == 0;
+    if (print) {
+        std::printf("// kAlgorithmVersion = %d\n", kAlgorithmVersion);
+        for (const auto& g : kGoldens)
+            std::printf("    {%d, 0x%016llxull, 0x%016llxull},\n", g.dim,
+                        static_cast<unsigned long long>(gaussianHash(g.dim)),
+                        static_cast<unsigned long long>(codesHash(g.dim)));
+        return 0;
+    }
+
+    int failures = 0;
+    for (const auto& g : kGoldens) {
+        const std::uint64_t gh = gaussianHash(g.dim);
+        const std::uint64_t ch = codesHash(g.dim);
+        if (gh != g.gaussian) {
+            std::printf("FAIL d=%d gaussian matrix: got 0x%016llx, want 0x%016llx\n", g.dim,
+                        static_cast<unsigned long long>(gh),
+                        static_cast<unsigned long long>(g.gaussian));
+            ++failures;
+        }
+        if (ch != g.codes) {
+            std::printf("FAIL d=%d codes: got 0x%016llx, want 0x%016llx\n", g.dim,
+                        static_cast<unsigned long long>(ch),
+                        static_cast<unsigned long long>(g.codes));
+            ++failures;
+        }
+    }
+    if (failures == 0) std::printf("PASS: %zu dims match goldens\n", std::size(kGoldens));
+    return failures == 0 ? 0 : 1;
+}