diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f1e0cef..7819b8c 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/BUILD.bazel b/BUILD.bazel index 11e6b8b..064a56c 100644 --- a/BUILD.bazel +++ b/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", +) diff --git a/README.md b/README.md index a2eb1d6..98e5f3d 100644 --- a/README.md +++ b/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. diff --git a/include/turboquant/types.h b/include/turboquant/types.h index d2c58ac..9cef5a5 100644 --- a/include/turboquant/types.h +++ b/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; diff --git a/src/rotation.cpp b/src/rotation.cpp index 2a09947..0dab8fa 100644 --- a/src/rotation.cpp +++ b/src/rotation.cpp @@ -3,33 +3,93 @@ #include "turboquant/rotation.h" +#include +#include +#include + 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(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(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((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& 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 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 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 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; } diff --git a/tests/determinism_test.cpp b/tests/determinism_test.cpp new file mode 100644 index 0000000..815dc8c --- /dev/null +++ b/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 + +#include +#include +#include +#include +#include +#include + +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(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(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(static_cast(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(gaussianHash(g.dim)), + static_cast(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(gh), + static_cast(g.gaussian)); + ++failures; + } + if (ch != g.codes) { + std::printf("FAIL d=%d codes: got 0x%016llx, want 0x%016llx\n", g.dim, + static_cast(ch), + static_cast(g.codes)); + ++failures; + } + } + if (failures == 0) std::printf("PASS: %zu dims match goldens\n", std::size(kGoldens)); + return failures == 0 ? 0 : 1; +}