diff --git a/README.md b/README.md index 5e21744..5a068ab 100644 --- a/README.md +++ b/README.md @@ -127,17 +127,25 @@ only the throughput differs. The implementation dispatches at first call to a kernel selected from CPU features via an atomic function pointer (libsodium-style -trampoline; race-free without locks or allocation). Today that -kernel is the per-ID scalar path on every supported host. The AVX2 -8-wide kernel that the issue tracker proposed -([#5](https://github.com/semantic-reasoning/libksuid/issues/5)) is -**not yet shipped** -- the implementation requires SIMD long- -division-by-62 with reciprocal-multiplication magic constants -verified against a parity corpus, and the dispatch infrastructure -landed in this PR is the foundation that will be added in a -follow-up. Until then, `ksuid_string_batch` is functionally -equivalent to a `ksuid_format` loop with one acquire-load + one -indirect call of overhead per call. +trampoline; race-free without locks or allocation): + +- **x86_64 + AVX2**: an 8-wide AVX2 kernel that processes eight + KSUIDs per outer iteration via a Granlund-Möller floor-reciprocal + multiply divide-by-62 ([#13](https://github.com/semantic-reasoning/libksuid/issues/13)). + The magic constant is auto-generated by `tools/derive-magic.py`, + pinned in `libksuid/divisor_magic.h`, and verified against + `__uint128_t` integer division on every CI run. +- **Other hosts** (non-AVX2 x86_64, aarch64, arm, ...): a per-ID + scalar loop equivalent to calling `ksuid_format` N times. + +Output is byte-identical across kernels; the differential parity +test in `tests/test_string_batch.c` cross-checks the AVX2 kernel +against the scalar reference over ≥ 2²⁰ pseudo-random KSUIDs and a +lane-swap detection corpus. + +Setting the environment variable `KSUID_FORCE_SCALAR=1` pins the +dispatcher to the scalar path at first call (runtime kill switch +without rebuilding the library). ## Layout diff --git a/libksuid/divisor_magic.h b/libksuid/divisor_magic.h new file mode 100644 index 0000000..ad82ef1 --- /dev/null +++ b/libksuid/divisor_magic.h @@ -0,0 +1,39 @@ +/* SPDX-License-Identifier: LGPL-3.0-or-later + * + * AUTO-GENERATED by tools/derive-magic.py 64 62. + * DO NOT HAND-EDIT. Re-run the deriver and commit the new output if + * the divisor or bit width changes. + * + * Granlund-Moeller magic constant for unsigned divide-by-62 via + * the high N bits of an N x N -> 2N multiply. + * + * For every value < 2^64 the reciprocal computation + * + * q' = (value * KSUID_DIV62_M) >> 64 + * r' = value - q' * 62 + * if (r' >= 62) { ++q'; r' -= 62; } + * + * yields the same q', r' as the straight integer division + * value / 62, value % 62. mulhi(value, M) may + * underestimate the true quotient by exactly 1 for some inputs and + * the correction step recovers it; mulhi never overestimates, + * which is why M is the FLOOR of 2^N / d (not the ceiling -- the + * ceiling form would overestimate for some inputs and the standard + * correction would not catch it). + * + * Property pinned at compile time in the consuming TU: + * 2^64 - M * 62 = 16 (must satisfy 0 <= ... <= d - 1) + * + * This file is regenerated by `tools/derive-magic.py` and verified + * by `tests/test_divisor_magic.c` on every CI run. + */ +#ifndef KSUID_DIVISOR_MAGIC_H +#define KSUID_DIVISOR_MAGIC_H + +#include + +#define KSUID_DIV62_M_BITS 64 +#define KSUID_DIV62_M ((uint64_t) UINT64_C (0x421084210842108)) +#define KSUID_DIV62_M_2N_MINUS_M_TIMES_D 16 + +#endif /* KSUID_DIVISOR_MAGIC_H */ diff --git a/libksuid/encode_avx2.c b/libksuid/encode_avx2.c new file mode 100644 index 0000000..63e47db --- /dev/null +++ b/libksuid/encode_avx2.c @@ -0,0 +1,216 @@ +/* SPDX-License-Identifier: LGPL-3.0-or-later + * + * AVX2 8-wide ksuid_string_batch kernel. + * + * Replaces the per-ID scalar `value / 62` and `value % 62` with a + * Granlund-Moeller multiply-high + correction step run in parallel + * across 8 KSUID lanes. The magic constant KSUID_DIV62_M = floor(2^64 + * / 62) lives in libksuid/divisor_magic.h (auto-generated by + * tools/derive-magic.py and verified against integer division by + * tests/test_divisor_magic.c on every CI run; FLOOR not ceiling -- + * see the deriver script for the reasoning). + * + * Lane layout (SoA). 8 KSUIDs x 5 base-2^32 limbs = 40 limb values. + * AVX2 _mm256_mul_epu32 takes the LOW 32 bits of each 64-bit lane + * (4 lanes per __m256i) and produces a 64-bit product per lane, so + * we carry each limb across two __m256i registers: the "lo" pack + * holds lanes 0..3 (one limb value per 64-bit lane) and the "hi" + * pack holds lanes 4..7. State: 5 limbs * 2 packs = 10 ymm + * registers. Haswell exposes 16 architectural ymm registers, so + * everything plus M, k62, and a few temporaries stays in registers. + * + * Outer iteration count is fixed at KSUID_STRING_LEN = 27. The + * scalar reference (libksuid/base62.c ksuid_base62_encode) suppresses + * leading zeros and stops when bp_len reaches 0, then '0'-pads the + * head. Lanes finish at different rates in SIMD, so we run the full + * 27 iterations unconditionally; once all 5 limbs of a lane are + * zero, every subsequent iteration produces (q=0, r=0) and writes + * '0' to that lane's output column -- which is exactly the head + * padding the scalar would have done. + * + * Tail (n & 7): falls through to ksuid_string_batch_scalar after + * the bulk loop. _mm256_zeroupper() is issued before any non-VEX + * code path is entered (Intel SDM Vol 1 sec 14.1.2 -- a stale + * upper-ymm causes a transition penalty when subsequent SSE-encoded + * instructions execute on the same lane). + * + * Issue #13 risk register cross-references: + * C2 (mulhi64 carry propagation) -- handled in + * ksuid_mulhi64_avx2 below; mid_low is the sum of three u32 + * values (fits in 34 bits, no u64 overflow) and the upper + * accumulator carries are added explicitly. + * C3 (mul_epu32 lane positioning) -- value is constructed so + * limb is in the LOW 32 of each 64-bit lane and remainder is + * in the HIGH 32; we issue four mul_epu32 invocations + * (lo*lo, lo*hi, hi*lo, hi*hi) to recover the full 64x64 + * product since each mul_epu32 only reads the low 32. + * C9 (zeroupper placement) -- emitted before scalar tail call + * AND before function return. + * R11 (force-scalar override) -- handled in encode_batch.c. + */ +#if defined(__x86_64__) || defined(_M_X64) +# include +# include +# include + +# include +# include +# include +# include + +/* Same alphabet as libksuid/base62.c; redeclared here to keep this + * TU self-contained (the scalar table is `static` to that file). */ +static const char k_avx2_b62_alphabet[64] + = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"; + +/* Multiply-high u64 x u64 -> high u64, 4 lanes wide. + * + * Schoolbook 64x64 -> 128 with the four 32x32 -> 64 cross-products. + * Let a = a_hi << 32 | a_lo, b = b_hi << 32 | b_lo. + * + * a*b = a_hi*b_hi << 64 + * + (a_lo*b_hi + a_hi*b_lo) << 32 + * + a_lo*b_lo + * + * The high 64 bits accumulate (a_hi*b_hi) plus the carry from the + * middle terms. mid_low sums three 32-bit values (lh_lo + hl_lo + + * ll_hi) which fits in 34 bits and so cannot overflow u64; + * mid_low >> 32 is the low-half-to-high-half carry; lh_hi and hl_hi + * are the upper 32 of the middle products that propagate directly + * into the high accumulator. + * + * Saturation check: with a_hi = b_hi = 2^32-1 and a_lo = b_lo = + * 2^32-1, the high accumulator equals (2^64-1) -- the maximum legal + * u64 -- so no internal step overflows. */ +static inline __m256i +ksuid_mulhi64_avx2 (__m256i a, __m256i b) +{ + __m256i mask32 = _mm256_set1_epi64x ((int64_t) 0xFFFFFFFFLL); + __m256i a_hi = _mm256_srli_epi64 (a, 32); + __m256i b_hi = _mm256_srli_epi64 (b, 32); + + __m256i ll = _mm256_mul_epu32 (a, b); + __m256i lh = _mm256_mul_epu32 (a, b_hi); + __m256i hl = _mm256_mul_epu32 (a_hi, b); + __m256i hh = _mm256_mul_epu32 (a_hi, b_hi); + + __m256i ll_hi = _mm256_srli_epi64 (ll, 32); + __m256i lh_lo = _mm256_and_si256 (lh, mask32); + __m256i hl_lo = _mm256_and_si256 (hl, mask32); + __m256i lh_hi = _mm256_srli_epi64 (lh, 32); + __m256i hl_hi = _mm256_srli_epi64 (hl, 32); + + __m256i mid_low = _mm256_add_epi64 (_mm256_add_epi64 (lh_lo, hl_lo), ll_hi); + __m256i mid_carry = _mm256_srli_epi64 (mid_low, 32); + + return _mm256_add_epi64 (_mm256_add_epi64 (hh, lh_hi), + _mm256_add_epi64 (hl_hi, mid_carry)); +} + +/* (q, r) = (value/62, value%62) per 64-bit lane via Granlund-Moeller: + * q = mulhi(value, M) -- may underestimate by 1 + * r = value - q*62 -- always in [0, 2*62-1] + * if (r >= 62) { ++q; r -= 62; } + * + * Domain in this kernel: value < 63 * 2^32 < 2^38 (limb in low 32, + * remainder in [0,61] in high 32). Therefore q < 2^33 and q*62 fits + * in 39 bits with room to spare. We compute q*62 = (q<<6) - (q<<1) + * with shifts -- AVX2 has no native u64*u64-to-u64 truncating + * multiply that would not require another mul_epu32 carry chain. */ +static inline void +ksuid_divmod62_avx2 (__m256i value, __m256i mag, __m256i *q_out, __m256i *r_out) +{ + __m256i q = ksuid_mulhi64_avx2 (value, mag); + __m256i q62 = + _mm256_sub_epi64 (_mm256_slli_epi64 (q, 6), _mm256_slli_epi64 (q, 1)); + __m256i r = _mm256_sub_epi64 (value, q62); + + /* _mm256_cmpgt_epi64 is signed; r and 61 are both small positive + * integers (< 2^7) so signed comparison agrees with unsigned. */ + __m256i ge62 = _mm256_cmpgt_epi64 (r, _mm256_set1_epi64x (61)); + __m256i k62 = _mm256_set1_epi64x (62); + __m256i one = _mm256_set1_epi64x (1); + *r_out = _mm256_sub_epi64 (r, _mm256_and_si256 (ge62, k62)); + *q_out = _mm256_add_epi64 (q, _mm256_and_si256 (ge62, one)); +} + +void +ksuid_string_batch_avx2 (const ksuid_t *ids, char *out_27n, size_t n) +{ + size_t bulk = n & ~(size_t) 7; + __m256i mag = _mm256_set1_epi64x ((int64_t) KSUID_DIV62_M); + + for (size_t base = 0; base < bulk; base += 8) { + /* SoA pack: 5 limbs x (lo-4-lanes, hi-4-lanes). The limb value + * lives in the low 32 of each 64-bit lane; the high 32 is zero + * initially and carries the running remainder during the inner + * loop. */ + __m256i Llo[5], Lhi[5]; + { + uint64_t plo[4], phi[4]; + const uint8_t *base_p = (const uint8_t *) &ids[base]; + for (size_t j = 0; j < 5; ++j) { + for (size_t lane = 0; lane < 4; ++lane) { + plo[lane] = + (uint64_t) ksuid_be32_load (base_p + lane * KSUID_BYTES + j * 4); + phi[lane] = + (uint64_t) ksuid_be32_load (base_p + (lane + 4) * KSUID_BYTES + + j * 4); + } + /* _mm256_loadu_si256 is the AVX2 unaligned-load intrinsic; + * the (__m256i *) cast is documented and does not require + * 32-byte alignment. */ + /* NOLINTNEXTLINE(clang-diagnostic-cast-align) */ + Llo[j] = _mm256_loadu_si256 ((const __m256i *) plo); + /* NOLINTNEXTLINE(clang-diagnostic-cast-align) */ + Lhi[j] = _mm256_loadu_si256 ((const __m256i *) phi); + } + } + + /* 27 outer iterations -- one per output character. The output + * column index walks right-to-left because long division by + * the radix produces digits LSB-first (the rightmost output + * character is the lowest-order base62 digit of the integer). */ + for (int col = KSUID_STRING_LEN - 1; col >= 0; --col) { + __m256i rem_lo = _mm256_setzero_si256 (); + __m256i rem_hi = _mm256_setzero_si256 (); + + for (size_t j = 0; j < 5; ++j) { + /* value = limb | (rem << 32). Both halves fit in u32 so + * the OR composes them without overlap. */ + __m256i v_lo = _mm256_or_si256 (Llo[j], _mm256_slli_epi64 (rem_lo, 32)); + __m256i v_hi = _mm256_or_si256 (Lhi[j], _mm256_slli_epi64 (rem_hi, 32)); + + __m256i q_lo, r_lo, q_hi, r_hi; + ksuid_divmod62_avx2 (v_lo, mag, &q_lo, &r_lo); + ksuid_divmod62_avx2 (v_hi, mag, &q_hi, &r_hi); + + Llo[j] = q_lo; + Lhi[j] = q_hi; + rem_lo = r_lo; + rem_hi = r_hi; + } + + /* Final remainder is the digit; alphabet-lookup, scatter. */ + uint64_t rems[8]; + /* NOLINTNEXTLINE(clang-diagnostic-cast-align) */ + _mm256_storeu_si256 ((__m256i *) & rems[0], rem_lo); + /* NOLINTNEXTLINE(clang-diagnostic-cast-align) */ + _mm256_storeu_si256 ((__m256i *) & rems[4], rem_hi); + for (size_t lane = 0; lane < 8; ++lane) { + out_27n[(base + lane) * KSUID_STRING_LEN + (size_t) col] + = k_avx2_b62_alphabet[rems[lane]]; + } + } + } + + /* Emit zeroupper before the scalar tail call AND before + * returning; both transitions can land us in non-VEX code. */ + _mm256_zeroupper (); + + if (bulk < n) + ksuid_string_batch_scalar (ids + bulk, + out_27n + bulk * KSUID_STRING_LEN, n - bulk); +} + +#endif /* x86_64 */ diff --git a/libksuid/encode_batch.c b/libksuid/encode_batch.c index ffbcec0..5c2decc 100644 --- a/libksuid/encode_batch.c +++ b/libksuid/encode_batch.c @@ -31,6 +31,8 @@ #include #include +#include +#include #if defined(__x86_64__) || defined(_M_X64) # if defined(__GNUC__) || defined(__clang__) @@ -99,13 +101,37 @@ ksuid_string_batch_init_trampoline (const ksuid_t * ids, char *out_27n, static _Atomic ksuid_string_batch_fn g_batch_impl = &ksuid_string_batch_init_trampoline; +/* KSUID_FORCE_SCALAR override (Critic R11). Reading getenv on the + * first dispatch only is safe -- the resolved pointer is cached for + * the lifetime of the process and the env var is consulted exactly + * once. The override exists so production deployments can pin the + * scalar path at startup if a future regression in the AVX2 kernel + * is discovered after rollout, without rebuilding the library. + * + * Recognised values: any non-empty, non-"0", non-"false" string + * disables the AVX2 kernel. NULL or unset = use the best kernel + * available on the host. */ +static int +ksuid_force_scalar_env (void) +{ + const char *v = getenv ("KSUID_FORCE_SCALAR"); + if (v == NULL || v[0] == '\0') + return 0; + if (strcmp (v, "0") == 0 || strcmp (v, "false") == 0 + || strcmp (v, "FALSE") == 0) + return 0; + return 1; +} + static void ksuid_string_batch_init_trampoline (const ksuid_t *ids, char *out_27n, size_t n) { ksuid_string_batch_fn resolved = &ksuid_string_batch_scalar; #if defined(KSUID_HAVE_AVX2_BATCH) - if (ksuid_cpu_supports_avx2 ()) + if (!ksuid_force_scalar_env () && ksuid_cpu_supports_avx2 ()) resolved = &ksuid_string_batch_avx2; +#else + (void) ksuid_force_scalar_env; /* silence unused-static warning */ #endif atomic_store_explicit (&g_batch_impl, resolved, memory_order_release); resolved (ids, out_27n, n); diff --git a/libksuid/ksuid.h b/libksuid/ksuid.h index 52a310f..67b8f2e 100644 --- a/libksuid/ksuid.h +++ b/libksuid/ksuid.h @@ -166,14 +166,20 @@ extern "C" * * The dispatch is resolved lazily on the first call (atomic, thread- * safe) and the resolved pointer is reused for the lifetime of the - * process. The eventual AVX2 8-wide kernel will plug into this - * dispatcher and deliver ~4-6x throughput on x86_64 hosts that - * support it; in the current release every host resolves to the - * per-ID scalar path (a ksuid_format loop), and the dispatch is - * one acquire-load + one indirect call of overhead per call. The - * AVX2 kernel itself is tracked as a follow-up to libksuid issue - * #5; until it ships, ksuid_string_batch is functionally - * equivalent to a ksuid_format loop with no measurable overhead. + * process. On x86_64 hosts that advertise AVX2, the dispatcher + * resolves to an 8-wide AVX2 kernel that processes 8 KSUIDs per + * outer iteration via a Granlund-Moeller reciprocal-multiply + * divide-by-62; on other hosts (including non-AVX2 x86_64) the + * dispatcher resolves to a per-ID scalar loop equivalent to + * ksuid_format() N times. The dispatch overhead -- one acquire- + * load + one indirect call per ksuid_string_batch invocation -- is + * lost in the noise compared to even the per-ID scalar work. + * + * Output is byte-identical across kernels. The KSUID_FORCE_SCALAR + * environment variable, if set to a non-empty non-"0"/"false" + * value, pins the dispatcher to the scalar path at first dispatch + * (a runtime kill switch for the AVX2 kernel without rebuilding + * the library). * * No error path: every 20-byte ksuid_t encodes by construction. n == 0 * is a no-op. The call is thread-safe for concurrent invocations on diff --git a/meson.build b/meson.build index 958bdc8..f6231bf 100644 --- a/meson.build +++ b/meson.build @@ -185,13 +185,57 @@ if get_option('simd') != 'none' endif endif +# AVX2 8-wide ksuid_string_batch kernel. Compiled in its own static +# library with -mavx2 (or /arch:AVX2 on MSVC) so the rest of libksuid +# keeps the SSE2 baseline ABI; the kernel is selected at runtime via +# CPUID inside the dispatcher (libksuid/encode_batch.c). Non-x86_64 +# hosts skip this entirely; non-AVX2 x86_64 hosts compile the TU but +# never call into it. +have_avx2_batch = false +avx2_arg = '' +if get_option('simd') != 'none' and host_cpu == 'x86_64' + if get_option('avx2_batch').allowed() + if cc.get_argument_syntax() == 'msvc' + candidate = '/arch:AVX2' + else + candidate = '-mavx2' + endif + if cc.has_argument(candidate) + have_avx2_batch = true + avx2_arg = candidate + common_args += '-DKSUID_HAVE_AVX2_BATCH=1' + endif + endif +endif + # Apply collected defines + warnings only after the per-arch SIMD -# section has had its chance to push KSUID_HAVE_SSE2 / _NEON. +# section has had its chance to push KSUID_HAVE_SSE2 / _NEON / +# _AVX2_BATCH. add_project_arguments(common_args, language : 'c') +avx2_lib_dep = [] +if have_avx2_batch + # PIC required so the same objects can land in both the static and + # shared variants of `both_libraries('ksuid', ...)` below. + ksuid_avx2_lib = static_library('ksuid_avx2_batch', + 'libksuid/encode_avx2.c', + include_directories : inc, + c_args : [avx2_arg], + pic : true, + install : false, + ) + # link_whole forces every object from the AVX2 static lib into the + # final library archive even though the call site (the trampoline + # in encode_batch.c) only ever takes the address of the kernel + # behind a runtime CPUID check; without link_whole the linker may + # GC the symbol on shared-library link. + avx2_lib_dep = [ksuid_avx2_lib] +endif + ksuid_lib = both_libraries('ksuid', core_sources, include_directories : inc, dependencies : extra_link_deps, + link_whole : avx2_lib_dep, version : meson.project_version(), soversion : '0', install : true, diff --git a/meson.options b/meson.options index 189dbd6..08fc12e 100644 --- a/meson.options +++ b/meson.options @@ -4,3 +4,5 @@ option('cli', type : 'boolean', value : true, description : 'Build the ksuid-gen CLI tool') option('simd', type : 'combo', choices : ['auto', 'none'], value : 'auto', description : 'Enable SIMD/NEON acceleration (auto-detected per host arch)') +option('avx2_batch', type : 'feature', value : 'auto', + description : 'Build the AVX2 8-wide ksuid_string_batch kernel (x86_64 only). Selected at runtime via CPUID; non-AVX2 hosts unaffected.') diff --git a/tests/meson.build b/tests/meson.build index e067012..82194d0 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -14,6 +14,26 @@ base_tests = ['test_smoke', 'test_parts', 'test_base62', 'test_parse_format', 'test_chacha20', 'test_new', 'test_simd_parity', 'test_wipe', 'test_compare_parity', 'test_string_batch'] +# test_divisor_magic verifies the auto-generated divisor_magic.h +# (Granlund-Moeller floor + correction reciprocal) against straight +# integer division for >=2^20 random inputs. The reference uses +# __uint128_t for mulhi64, which GCC and Clang provide on every +# 64-bit target; MSVC does not, so the test only registers when the +# compiler advertises the type. The AVX2 kernel that consumes the +# same magic constant will only be built where __uint128_t is +# available anyway. +have_uint128 = cc.has_type('__uint128_t') or cc.has_type('unsigned __int128') +# test_divisor_magic uses __uint128_t for its mulhi64 reference, +# which is a GCC/Clang extension flagged by -Wpedantic. Suppress +# the pedantic warning for this single test TU rather than dropping +# warning_level=3 across the whole build. +if have_uint128 + test_divisor_magic_args = [] + if cc.has_argument('-Wno-pedantic') + test_divisor_magic_args += '-Wno-pedantic' + endif +endif + foreach t : base_tests exe = executable(t, t + '.c', include_directories : test_inc, @@ -22,6 +42,15 @@ foreach t : base_tests test(t, exe, timeout : 30) endforeach +if have_uint128 + exe = executable('test_divisor_magic', 'test_divisor_magic.c', + include_directories : test_inc, + link_with : test_link, + c_args : test_divisor_magic_args, + ) + test('test_divisor_magic', exe, timeout : 30) +endif + if have_threads_h and threads_dep.found() # test_rand_tls compiles with -DKSUID_TESTING=1 only when the # platform exposes a thread-exit hook to assert against (issue #4 diff --git a/tests/test_divisor_magic.c b/tests/test_divisor_magic.c new file mode 100644 index 0000000..ab5b929 --- /dev/null +++ b/tests/test_divisor_magic.c @@ -0,0 +1,144 @@ +/* SPDX-License-Identifier: LGPL-3.0-or-later + * + * Verifies that the auto-generated KSUID_DIV62_M magic constant in + * libksuid/divisor_magic.h, plugged into a Granlund-Moeller scalar + * reciprocal, agrees with straight integer division by 62 for every + * tested input. Builds the floor/correction algorithm from scratch + * here -- the AVX2 kernel will reuse the same algebra over 8 lanes. + * + * Coverage: + * - pinned corner cases (0, 1, 61, 62, 63, UINT64_MAX, ...) + * - 2^20 random u64 inputs (matches the M2 acceptance criterion + * in issue #13: "differential parity test seed must be + * reproducible") via a seeded LCG. + * + * Failure here means either the deriver script is wrong or the + * checked-in header was hand-edited. Both modes are blockers for + * shipping the AVX2 kernel. + */ +#include +#include "test_util.h" + +/* Reference mulhi64 via __uint128_t. Available on GCC/Clang on + * every libksuid CI lane that exercises this test (the AVX2-target + * x86_64 lane is the only one that ships the AVX2 kernel; MSVC is + * not in this test's build matrix). */ +_Static_assert (sizeof (unsigned __int128) == 16, "test requires __uint128_t"); + +/* Compile-time pin of the deriver's contract. If the header was + * hand-edited or regenerated for a different divisor, this fails + * the build before any test runs. */ +_Static_assert (KSUID_DIV62_M_BITS == 64, + "KSUID_DIV62_M expected to be 64-bit"); +_Static_assert (KSUID_DIV62_M * (uint64_t) 62 + + KSUID_DIV62_M_2N_MINUS_M_TIMES_D == 0, + "KSUID_DIV62_M does not satisfy 2^64 - M * 62 = " + "KSUID_DIV62_M_2N_MINUS_M_TIMES_D"); +_Static_assert (KSUID_DIV62_M_2N_MINUS_M_TIMES_D < 62, + "deficit must be in [0, d - 1]; " + "anything else means M is not floor(2^64 / 62)"); + +static uint64_t +mulhi64_ref (uint64_t a, uint64_t b) +{ + return (uint64_t) (((unsigned __int128) a * b) >> 64); +} + +static uint64_t +div62_via_magic (uint64_t value) +{ + uint64_t q = mulhi64_ref (value, KSUID_DIV62_M); + uint64_t r = value - q * 62; + /* mulhi may underestimate by exactly 1; never overestimates + * because M is the FLOOR. The correction is unconditional in + * the sense that r is always in [0, 123]; the branch picks the + * exact quotient + remainder. We only need q here -- mod62 has + * its own helper -- so the post-correction r is unused, but + * we must run the branch to update q. */ + if (r >= 62) + q += 1; + return q; +} + +static uint64_t +mod62_via_magic (uint64_t value) +{ + uint64_t q = mulhi64_ref (value, KSUID_DIV62_M); + uint64_t r = value - q * 62; + if (r >= 62) + r -= 62; + return r; +} + +static void +check_one (uint64_t value) +{ + uint64_t q_ref = value / 62; + uint64_t r_ref = value % 62; + uint64_t q_mag = div62_via_magic (value); + uint64_t r_mag = mod62_via_magic (value); + if (q_mag != q_ref || r_mag != r_ref) { + fprintf (stderr, + " divmod62 mismatch at value=0x%016llx:\n" + " reference: q=%llu, r=%llu\n" + " magic: q=%llu, r=%llu\n", + (unsigned long long) value, + (unsigned long long) q_ref, (unsigned long long) r_ref, + (unsigned long long) q_mag, (unsigned long long) r_mag); + ksuid_test_failures_++; + } +} + +static void +test_pinned_corners (void) +{ + /* Boundaries: 0, around d, around 2^32, around 2^38 (the AVX2 + * kernel's value bound), around 2^63, max u64. */ + static const uint64_t corners[] = { + 0, 1, 60, 61, 62, 63, 123, 124, 125, + (uint64_t) UINT32_MAX - 1, (uint64_t) UINT32_MAX, + (uint64_t) UINT32_MAX + 1, (uint64_t) UINT32_MAX + 62, + (uint64_t) 1 << 32, ((uint64_t) 1 << 32) + 1, + /* The exact AVX2 in-loop bound: value < 62 * 2^32 + 2^32 */ + ((uint64_t) 62 << 32) - 1, ((uint64_t) 62 << 32), + ((uint64_t) 63 << 32) - 1, + /* Larger values still must work for full-u64 confidence. */ + (uint64_t) 1 << 50, (uint64_t) 1 << 60, (uint64_t) 1 << 63, + UINT64_MAX - 62, UINT64_MAX - 61, UINT64_MAX - 1, UINT64_MAX, + }; + for (size_t i = 0; i < sizeof corners / sizeof corners[0]; ++i) + check_one (corners[i]); +} + +static void +test_one_million_lcg_random (void) +{ + /* Seeded LCG -- a failing input found in CI must reproduce + * locally with the same seed. Period 2^64. 2^20 = 1048576 + * samples is the M2 acceptance threshold from issue #13. */ + uint64_t s = 0x9e3779b97f4a7c15ULL; + for (size_t trial = 0; trial < (1u << 20); ++trial) { + s = s * 6364136223846793005ULL + 1442695040888963407ULL; + check_one (s); + } +} + +static void +test_dense_low_range (void) +{ + /* Exhaustive sweep of every value in [0, 4 * 62) so every + * remainder value 0..61 is exercised at every quotient + * 0, 1, 2, 3. Catches off-by-one errors that random testing + * misses for low quotients. */ + for (uint64_t v = 0; v < (uint64_t) 4 * 62; ++v) + check_one (v); +} + +int +main (void) +{ + RUN_TEST (test_pinned_corners); + RUN_TEST (test_dense_low_range); + RUN_TEST (test_one_million_lcg_random); + TEST_MAIN_END (); +} diff --git a/tests/test_string_batch.c b/tests/test_string_batch.c index e25e2e4..93e1f7c 100644 --- a/tests/test_string_batch.c +++ b/tests/test_string_batch.c @@ -1,18 +1,42 @@ /* SPDX-License-Identifier: LGPL-3.0-or-later * - * Tests for the public ksuid_string_batch bulk encoder. This commit - * lands the API + the scalar reference; the AVX2 8-wide kernel - * lands in a follow-up commit. The differential parity test against - * the AVX2 kernel arrives in commit 4. For now we pin the contract: - * - n == 0 is a no-op - * - n KSUIDs land at the documented output offsets - * - every produced 27-byte slice equals ksuid_format of the same ID + * Tests for the public ksuid_string_batch bulk encoder. + * + * Two layers of coverage: + * 1. Public-API parity: every 27-byte slice produced by + * ksuid_string_batch (which dispatches to the best kernel for + * the host -- AVX2 on AVX2 x86_64, scalar elsewhere) equals + * ksuid_format of the same ID. Pins n=0 no-op, exact-multiple- + * of-8, off-by-one-into-tail (n=9), prime-misaligned (n=257), + * and the corner KSUIDs (NIL, MAX). + * 2. Direct AVX2-vs-scalar differential parity (compiled in only + * when KSUID_HAVE_AVX2_BATCH is defined; gated at runtime on + * __builtin_cpu_supports("avx2") so the same binary is safe + * on non-AVX2 hosts in the same x86_64 build). Bypasses the + * runtime dispatcher and calls the scalar + AVX2 kernels + * directly, comparing byte-for-byte. This is what catches + * cross-lane bugs (Critic R3 in issue #13: 8 distinct KSUIDs + * whose lanes get swapped by an off-by-one in the SoA pack + * would still pass per-ID format-parity tests if the wrong + * output happens to match a different input). */ #include #include "test_util.h" #include +#if defined(KSUID_HAVE_AVX2_BATCH) && (defined(__GNUC__) || defined(__clang__)) +# define KSUID_TEST_AVX2_PARITY 1 +/* Internal kernel prototypes. Tests link against the static archive + * so default-hidden visibility does not exclude these symbols. */ +extern void ksuid_string_batch_scalar (const ksuid_t * ids, char *out_27n, + size_t n); +extern void ksuid_string_batch_avx2 (const ksuid_t * ids, char *out_27n, + size_t n); +#else +# define KSUID_TEST_AVX2_PARITY 0 +#endif + static void fill_pseudo_random (ksuid_t *id, uint64_t seed) { @@ -111,6 +135,168 @@ test_batch_pinned_corners (void) } } +#if KSUID_TEST_AVX2_PARITY +static int +host_supports_avx2 (void) +{ + __builtin_cpu_init (); + return __builtin_cpu_supports ("avx2"); +} + +static void +avx2_parity_for_n (size_t n) +{ + if (!host_supports_avx2 ()) + return; + ksuid_t *ids = malloc (n * sizeof *ids); + ASSERT_TRUE (ids != NULL); + for (size_t i = 0; i < n; ++i) + fill_pseudo_random (&ids[i], + 0xa3b1c2d4e5f60718ULL ^ (i * 0x9e3779b97f4a7c15ULL)); + + char *out_s = malloc (n * KSUID_STRING_LEN); + char *out_a = malloc (n * KSUID_STRING_LEN); + if (out_s == NULL || out_a == NULL) { + FAIL_ ("out_s/out_a malloc"); + free (ids); + free (out_s); + free (out_a); + return; + } + ksuid_string_batch_scalar (ids, out_s, n); + ksuid_string_batch_avx2 (ids, out_a, n); + + /* Per-lane byte compare so the failure message identifies WHICH + * KSUID position diverged (a single ASSERT_EQ_BYTES across the + * full buffer would only print the first byte offset). */ + for (size_t i = 0; i < n; ++i) + ASSERT_EQ_BYTES (out_s + i * KSUID_STRING_LEN, + out_a + i * KSUID_STRING_LEN, KSUID_STRING_LEN); + + free (ids); + free (out_s); + free (out_a); +} + +static void +test_avx2_parity_n_in_block_boundaries (void) +{ + /* Boundaries around the 8-wide block size: tail-only, exact + * block, off-by-one into tail, two blocks plus tail, etc. */ + static const size_t ns[] = { 1, 7, 8, 9, 15, 16, 17, 23, 24, 25, 1000 }; + for (size_t i = 0; i < sizeof ns / sizeof ns[0]; ++i) + avx2_parity_for_n (ns[i]); +} + +static void +test_avx2_parity_lane_swap_detection (void) +{ + /* 8 distinct KSUIDs in the same vector. If the AVX2 lane-pack + * code mis-mapped lane k to lane k', the out-of-position + * comparison against the scalar reference fails. */ + if (!host_supports_avx2 ()) + return; + ksuid_t ids[8]; + for (size_t lane = 0; lane < 8; ++lane) + fill_pseudo_random (&ids[lane], + 0xdeadbeefcafef00dULL + (uint64_t) (lane * 0x100000001b3ULL)); + + char out_s[8 * KSUID_STRING_LEN]; + char out_a[8 * KSUID_STRING_LEN]; + ksuid_string_batch_scalar (ids, out_s, 8); + ksuid_string_batch_avx2 (ids, out_a, 8); + for (size_t lane = 0; lane < 8; ++lane) + ASSERT_EQ_BYTES (out_s + lane * KSUID_STRING_LEN, + out_a + lane * KSUID_STRING_LEN, KSUID_STRING_LEN); +} + +static void +test_avx2_parity_corner_values (void) +{ + if (!host_supports_avx2 ()) + return; + /* 16 KSUIDs of pure-NIL / pure-MAX / all-0xff payload variants + * spread across two SIMD blocks, exercising the long-division + * "all-zero limbs" and "max-value limbs" extremes for every + * lane position in the SoA layout. */ + ksuid_t ids[16]; + for (size_t i = 0; i < 16; ++i) { + if ((i & 3) == 0) + ids[i] = KSUID_NIL; + else if ((i & 3) == 1) + ids[i] = KSUID_MAX; + else if ((i & 3) == 2) { + /* All-bytes-0x80, exercises the high bit being set in + * every limb which would expose any signed/unsigned + * confusion in mulhi64. */ + for (size_t j = 0; j < KSUID_BYTES; ++j) + ids[i].b[j] = 0x80; + } else { + /* Limbs alternating 0xffffffff / 0x00000000 -- exercises + * the rem-injection path with non-trivial high bits. */ + for (size_t j = 0; j < KSUID_BYTES; ++j) + ids[i].b[j] = ((j / 4) % 2) ? 0xff : 0x00; + } + } + char out_s[16 * KSUID_STRING_LEN]; + char out_a[16 * KSUID_STRING_LEN]; + ksuid_string_batch_scalar (ids, out_s, 16); + ksuid_string_batch_avx2 (ids, out_a, 16); + for (size_t i = 0; i < 16; ++i) + ASSERT_EQ_BYTES (out_s + i * KSUID_STRING_LEN, + out_a + i * KSUID_STRING_LEN, KSUID_STRING_LEN); +} + +static void +test_avx2_parity_one_million_lcg (void) +{ + /* Match the M2 acceptance threshold from issue #13: >= 2^20 + * pseudo-random KSUIDs differential-checked end-to-end. The seed + * is the same constant as test_divisor_magic.c so a CI failure + * reproduces locally bit-for-bit. */ + if (!host_supports_avx2 ()) + return; + size_t n = 1u << 20; + ksuid_t *ids = malloc (n * sizeof *ids); + ASSERT_TRUE (ids != NULL); + uint64_t s = 0x9e3779b97f4a7c15ULL; + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < KSUID_BYTES; ++j) { + s = s * 6364136223846793005ULL + 1442695040888963407ULL; + ids[i].b[j] = (uint8_t) (s >> 56); + } + } + + char *out_s = malloc (n * KSUID_STRING_LEN); + char *out_a = malloc (n * KSUID_STRING_LEN); + if (out_s == NULL || out_a == NULL) { + FAIL_ ("out_s/out_a malloc"); + free (ids); + free (out_s); + free (out_a); + return; + } + ksuid_string_batch_scalar (ids, out_s, n); + ksuid_string_batch_avx2 (ids, out_a, n); + + /* memcmp-style fast check; only fall through to per-position + * report on mismatch to keep the success path cheap. */ + if (memcmp (out_s, out_a, n * KSUID_STRING_LEN) != 0) { + for (size_t i = 0; i < n; ++i) { + if (memcmp (out_s + i * KSUID_STRING_LEN, + out_a + i * KSUID_STRING_LEN, KSUID_STRING_LEN) != 0) { + fprintf (stderr, " AVX2 parity diverged at lane %zu of %zu\n", i, n); + ksuid_test_failures_++; + break; + } + } + } + free (ids); + free (out_s); + free (out_a); +} +#endif /* KSUID_TEST_AVX2_PARITY */ + int main (void) { @@ -122,5 +308,11 @@ main (void) RUN_TEST (test_batch_64); RUN_TEST (test_batch_257_misaligned); RUN_TEST (test_batch_pinned_corners); +#if KSUID_TEST_AVX2_PARITY + RUN_TEST (test_avx2_parity_n_in_block_boundaries); + RUN_TEST (test_avx2_parity_lane_swap_detection); + RUN_TEST (test_avx2_parity_corner_values); + RUN_TEST (test_avx2_parity_one_million_lcg); +#endif TEST_MAIN_END (); } diff --git a/tools/derive-magic.py b/tools/derive-magic.py new file mode 100755 index 0000000..6b15149 --- /dev/null +++ b/tools/derive-magic.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Granlund-Moeller magic-constant deriver. Generates the M constant +# such that for every value < 2^N, value / d == (value * M) >> N +# (with at most one off-by-one correction step using value - q*d >= d). +# +# Usage: +# tools/derive-magic.py [] +# +# Example: +# tools/derive-magic.py 64 62 libksuid/divisor_magic.h +# +# This script exists because every architect/critic in the issue-#13 +# review chain hand-derived a different (and wrong) constant for the +# divide-by-62 reciprocal. Sourcing the magic from a verified +# programmatic generator, with the resulting header committed to the +# tree alongside its generation arguments and a _Static_assert in +# the consuming TU, removes that class of error. +# +# Mathematical contract for M = floor(2^N / d): +# - M * d <= 2^N (definition of floor) +# - M * d >= 2^N - (d - 1) (because floor of integer division by +# d differs from ceiling by at most 1, scaled by d) +# These together imply 0 <= 2^N - M * d <= d - 1, which is the +# precondition under which mulhi(v, M) >= v/d - 1 (mulhi may +# underestimate by exactly 1, never overestimates) and the +# standard Granlund-Moeller correction step +# r = v - q * d; if (r >= d) { ++q; r -= d; } +# recovers the exact quotient + remainder. +# +# Note: ceiling(2^N / d) does NOT satisfy this contract -- with +# ceiling, mulhi(v, M) can OVERESTIMATE for large v and the +# correction's "if (r >= d)" check fires when it shouldn't, +# producing an off-by-one quotient. The earlier issue-#13 +# discussion landed on the wrong constant precisely because +# multiple reviewers conflated ceiling vs floor; see the issue +# comment for the failure mode. + +import sys + + +def derive_magic(n_bits: int, divisor: int) -> int: + if not (1 < divisor < (1 << n_bits)): + raise ValueError(f"divisor must satisfy 1 < d < 2^{n_bits}; got {divisor}") + two_n = 1 << n_bits + # floor(2^N / d): M * d is at most 2^N, the standard + # Granlund-Moeller "underestimate" form that pairs with the + # post-multiply remainder-based correction step. + m = two_n // divisor + # Sanity: 2^N - M * d is in [0, d - 1] (i.e. 2^N mod d). + deficit = two_n - m * divisor + if not (0 <= deficit <= divisor - 1): + raise AssertionError( + f"derived M={m:#x} fails precondition: 2^N - M*d = {deficit} " + f"not in [0, {divisor - 1}]" + ) + return m + + +def verify_against_reference(n_bits: int, divisor: int, m: int, samples: int) -> None: + """Cross-check the magic constant by running the Granlund-Moeller scalar + reciprocal across `samples` random values and asserting it equals + the straight integer division. Fails loudly if M is wrong.""" + import random + rng = random.Random(0x9E3779B97F4A7C15) + mask = (1 << n_bits) - 1 + for _ in range(samples): + v = rng.randint(0, mask) + q_ref = v // divisor + # mulhi: high N bits of v * M + q_mul = (v * m) >> n_bits + r = v - q_mul * divisor + if r >= divisor: + q_mul += 1 + r -= divisor + if q_mul != q_ref: + raise AssertionError( + f"magic constant verification failed at v={v:#x}: " + f"reference quotient {q_ref}, magic {q_mul}" + ) + + +def emit_header(path: str, n_bits: int, divisor: int, m: int) -> None: + deficit = (1 << n_bits) - m * divisor + text = f"""/* SPDX-License-Identifier: LGPL-3.0-or-later + * + * AUTO-GENERATED by tools/derive-magic.py {n_bits} {divisor}. + * DO NOT HAND-EDIT. Re-run the deriver and commit the new output if + * the divisor or bit width changes. + * + * Granlund-Moeller magic constant for unsigned divide-by-{divisor} via + * the high N bits of an N x N -> 2N multiply. + * + * For every value < 2^{n_bits} the reciprocal computation + * + * q' = (value * KSUID_DIV{divisor}_M) >> {n_bits} + * r' = value - q' * {divisor} + * if (r' >= {divisor}) {{ ++q'; r' -= {divisor}; }} + * + * yields the same q', r' as the straight integer division + * value / {divisor}, value % {divisor}. mulhi(value, M) may + * underestimate the true quotient by exactly 1 for some inputs and + * the correction step recovers it; mulhi never overestimates, + * which is why M is the FLOOR of 2^N / d (not the ceiling -- the + * ceiling form would overestimate for some inputs and the standard + * correction would not catch it). + * + * Property pinned at compile time in the consuming TU: + * 2^{n_bits} - M * {divisor} = {deficit} (must satisfy 0 <= ... <= d - 1) + * + * This file is regenerated by `tools/derive-magic.py` and verified + * by `tests/test_divisor_magic.c` on every CI run. + */ +#ifndef KSUID_DIVISOR_MAGIC_H +#define KSUID_DIVISOR_MAGIC_H + +#include + +#define KSUID_DIV{divisor}_M_BITS {n_bits} +#define KSUID_DIV{divisor}_M ((uint{n_bits}_t) UINT{n_bits}_C ({m:#x})) +#define KSUID_DIV{divisor}_M_2N_MINUS_M_TIMES_D {deficit} + +#endif /* KSUID_DIVISOR_MAGIC_H */ +""" + with open(path, "w", encoding="utf-8", newline="\n") as fh: + fh.write(text) + + +def main(argv: list[str]) -> int: + if not (3 <= len(argv) <= 4): + sys.stderr.write( + "usage: derive-magic.py []\n" + ) + return 2 + n_bits = int(argv[1]) + divisor = int(argv[2]) + out_path = argv[3] if len(argv) == 4 else None + + m = derive_magic(n_bits, divisor) + verify_against_reference(n_bits, divisor, m, samples=1 << 16) + + if out_path is None: + deficit = (1 << n_bits) - m * divisor + print(f"M = {m:#x} (2^{n_bits} - M * {divisor} = {deficit})") + return 0 + + emit_header(out_path, n_bits, divisor, m) + return 0 + + +if __name__ == "__main__": + sys.exit(main(sys.argv))