From 9f5fc17e32bf5c7c6cfadf40c29a1dedb4cc03ac Mon Sep 17 00:00:00 2001 From: Lasse Collin Date: Mon, 10 Jun 2024 14:45:44 +0300 Subject: [PATCH] liblzma: Add crc_clmul_consts_gen.c It's a standalone program that prints the required constants. It's won't be a part of the normal build of the package. --- src/liblzma/check/Makefile.inc | 1 + src/liblzma/check/crc_clmul_consts_gen.c | 160 +++++++++++++++++++++++ 2 files changed, 161 insertions(+) create mode 100644 src/liblzma/check/crc_clmul_consts_gen.c diff --git a/src/liblzma/check/Makefile.inc b/src/liblzma/check/Makefile.inc index 04cf1ffc..6c0e3739 100644 --- a/src/liblzma/check/Makefile.inc +++ b/src/liblzma/check/Makefile.inc @@ -5,6 +5,7 @@ ## currently crc32 is always enabled. EXTRA_DIST += \ + check/crc_clmul_consts_gen.c \ check/crc32_tablegen.c \ check/crc64_tablegen.c diff --git a/src/liblzma/check/crc_clmul_consts_gen.c b/src/liblzma/check/crc_clmul_consts_gen.c new file mode 100644 index 00000000..5fe14bd6 --- /dev/null +++ b/src/liblzma/check/crc_clmul_consts_gen.c @@ -0,0 +1,160 @@ +// SPDX-License-Identifier: 0BSD + +/////////////////////////////////////////////////////////////////////////////// +// +/// \file crc_clmul_consts_gen.c +/// \brief Generate constants for CLMUL CRC code +/// +/// Compiling: gcc -std=c99 -o crc_clmul_consts_gen crc_clmul_consts_gen.c +/// +/// This is for CRCs that use reversed bit order (bit reflection). +/// The same CLMUL CRC code can be used with CRC64 and smaller ones like +/// CRC32 apart from one special case: CRC64 needs an extra step in the +/// Barrett reduction to handle the 65th bit; the smaller ones don't. +/// Otherwise it's enough to just change the polynomial and the derived +/// constants and use the same code. +/// +/// See the Intel white paper "Fast CRC Computation for Generic Polynomials +/// Using PCLMULQDQ Instruction" from 2009. +// +// Author: Lasse Collin +// +/////////////////////////////////////////////////////////////////////////////// + +#include +#include + + +/// CRC32 (Ethernet) polynomial in reversed representation +static const uint64_t p32 = 0xedb88320; + +// CRC64 (ECMA-182) polynomial in reversed representation +static const uint64_t p64 = 0xc96c5795d7870f42; + + +/// Calculates floor(x^128 / p) where p is a CRC64 polynomial in +/// reversed representation. The result is in reversed representation too. +static uint64_t +calc_cldiv(uint64_t p) +{ + // Quotient + uint64_t q = 0; + + // Align the x^64 term with the x^128 (the implied high bits of the + // divisor and the dividend) and do the first step of polynomial long + // division, calculating the first remainder. The variable q remains + // zero because the highest bit of the quotient is an implied bit 1 + // (we kind of set q = 1 << -1). + uint64_t r = p; + + // Then process the remaining 64 terms. Note that r has no implied + // high bit, only q and p do. (And remember that a high bit in the + // polynomial is stored at a low bit in the variable due to the + // reversed bit order.) + for (unsigned i = 0; i < 64; ++i) { + q |= (r & 1) << i; + r = (r >> 1) ^ (r & 1 ? p : 0); + } + + return q; +} + + +/// Calculate the remainder of carryless division: +/// +/// x^(bits + n - 1) % p, where n=64 (for CRC64) +/// +/// p must be in reversed representation which omits the bit of +/// the highest term of the polynomial. Instead, it is an implied bit +/// at kind of like "1 << -1" position, as if it had just been shifted out. +/// +/// The return value is in the reversed bit order. (There are no implied bits.) +static uint64_t +calc_clrem(uint64_t p, unsigned bits) +{ + // Do the first step of polynomial long division. + uint64_t r = p; + + // Then process the remaining terms. Start with i = 1 instead of i = 0 + // to account for the -1 in x^(bits + n - 1). This -1 is convenient + // with the reversed bit order. See the "Bit-Reflection" section in + // the Intel white paper. + for (unsigned i = 1; i < bits; ++i) + r = (r >> 1) ^ (r & 1 ? p : 0); + + return r; +} + + +extern int +main(void) +{ + puts("// CRC64"); + + // The order of the two 64-bit constants in a vector don't matter. + // It feels logical to put them in this order as it matches the + // order in which the input bytes are read. + printf("const __m128i fold512 = _mm_set_epi64x(" + "0x%016" PRIx64 ", 0x%016" PRIx64 ");\n", + calc_clrem(p64, 4 * 128 - 64), + calc_clrem(p64, 4 * 128)); + + printf("const __m128i fold128 = _mm_set_epi64x(" + "0x%016" PRIx64 ", 0x%016" PRIx64 ");\n", + calc_clrem(p64, 128 - 64), + calc_clrem(p64, 128)); + + // When we multiply by mu, we care about the high bits of the result + // (in reversed bit order!). It doesn't matter that the low bit gets + // shifted out because the affected output bits will be ignored. + // Below we add the implied high bit with "| 1" after the shifting + // so that the high bits of the multiplication will be correct. + // + // p64 is shifted left by one so that the final multiplication + // in Barrett reduction won't be misaligned by one bit. We could + // use "(p64 << 1) | 1" instead of "p64 << 1" too but it makes + // no difference as that bit won't affect the relevant output bits + // (we only care about the lowest 64 bits of the result, that is, + // lowest in the reversed bit order). + // + // NOTE: The 65rd bit of p64 gets shifted out. It needs to be + // compensated with 64-bit shift and xor in the CRC64 code. + printf("const __m128i mu_p = _mm_set_epi64x(" + "0x%016" PRIx64 ", 0x%016" PRIx64 ");\n", + (calc_cldiv(p64) << 1) | 1, + p64 << 1); + + puts(""); + + puts("// CRC32"); + + printf("const __m128i fold512 = _mm_set_epi64x(" + "0x%08" PRIx64 ", 0x%08" PRIx64 ");\n", + calc_clrem(p32, 4 * 128 - 64), + calc_clrem(p32, 4 * 128)); + + printf("const __m128i fold128 = _mm_set_epi64x(" + "0x%08" PRIx64 ", 0x%08" PRIx64 ");\n", + calc_clrem(p32, 128 - 64), + calc_clrem(p32, 128)); + + // CRC32 calculation is done by modulus scaling it to a CRC64. + // Since the CRC is in reversed representation, only the mu + // constant changes with the modulus scaling. This method avoids + // one additional constant and one additional clmul in the final + // reduction steps, making the code both simpler and faster. + // + // p32 is shifted left by one so that the final multiplication + // in Barrett reduction won't be misaligned by one bit. We could + // use "(p32 << 1) | 1" instead of "p32 << 1" too but it makes + // no difference as that bit won't affect the relevant output bits. + // + // NOTE: The 33-bit value fits in 64 bits so, unlike with CRC64, + // there is no need to compensate for any missing bits in the code. + printf("const __m128i mu_p = _mm_set_epi64x(" + "0x%016" PRIx64 ", 0x%" PRIx64 ");\n", + (calc_cldiv(p32) << 1) | 1, + p32 << 1); + + return 0; +}