XrdOucCRC32C.cc

Go to the documentation of this file.

/* crc32c.c -- compute CRC-32C using the Intel crc32 instruction
 * Copyright (C) 2013, 2015 Mark Adler
 * Version 1.3  31 Dec 2015  Mark Adler
 */
 
/*
  This software is provided 'as-is', without any express or implied
  warranty.  In no event will the author be held liable for any damages
  arising from the use of this software.
 
  Permission is granted to anyone to use this software for any purpose,
  including commercial applications, and to alter it and redistribute it
  freely, subject to the following restrictions:
 
  1. The origin of this software must not be misrepresented; you must not
     claim that you wrote the original software. If you use this software
     in a product, an acknowledgment in the product documentation would be
     appreciated but is not required.
  2. Altered source versions must be plainly marked as such, and must not be
     misrepresented as being the original software.
  3. This notice may not be removed or altered from any source distribution.
 
  Mark Adler
  madler@alumni.caltech.edu
 */
 
/* Use hardware CRC instruction on Intel SSE 4.2 processors.  This computes a
   CRC-32C, *not* the CRC-32 used by Ethernet and zip, gzip, etc.  A software
   version is provided as a fall-back, as well as for speed comparisons. */
 
/* Version history:
   1.0  10 Feb 2013  First version
   1.1   1 Aug 2013  Correct comments on why three crc instructions in parallel
   1.2   1 Nov 2015  Add const qualifier to avoid compiler warning
                     Load entire input into memory (test code)
                     Argument gives number of times to repeat (test code)
                     Argument < 0 forces software implementation (test code)
   1.3  31 Dec 2015  Check for Intel architecture using compiler macro
                     Support big-endian processors in software calculation
                     Add header for external use
 */
 
/* Modification history:
        22 Nov 2019  Rename orignal crc32.c to XrdOucCRC32.cc and crc32c.h to
                     XrdOucCRC32C.hh with corresponding change to include
                     statement herein. Add required casts to allow C++
                     compilation.
        28 Oct 2026  Add a hardware version for ARMv8 processors.
 */
 
#include <pthread.h>
#include "XrdOuc/XrdOucCRC32C.hh"
 
/* Little-endian 64-bit ARM processors may provide the ARMv8 crc32c
   instructions, mandatory since ARMv8.1 and optional in ARMv8.0. */
#if defined(__aarch64__) && !defined(__AARCH64EB__)
#define CRC32C_HW_ARM
#include <arm_acle.h>
#include <cstring>
#if !defined(__ARM_FEATURE_CRC32) && defined(__linux__)
#include <sys/auxv.h>
#include <asm/hwcap.h>
#endif
#endif
 
/* CRC-32C (iSCSI) polynomial in reversed bit order. */
#define POLY 0x82f63b78
 
#if defined(__x86_64__) || defined(CRC32C_HW_ARM)
 
/* Shift tables shared by the hardware CRC-32C implementations. */
 
/* Multiply a matrix times a vector over the Galois field of two elements,
   GF(2).  Each element is a bit in an unsigned integer.  mat must have at
   least as many entries as the power of two for most significant one bit in
   vec. */
static inline uint32_t gf2_matrix_times(uint32_t *mat, uint32_t vec) {
    uint32_t sum = 0;
    while (vec) {
        if (vec & 1)
            sum ^= *mat;
        vec >>= 1;
        mat++;
    }
    return sum;
}
 
/* Multiply a matrix by itself over GF(2).  Both mat and square must have 32
   rows. */
static inline void gf2_matrix_square(uint32_t *square, uint32_t *mat) {
    for (unsigned n = 0; n < 32; n++)
        square[n] = gf2_matrix_times(mat, mat[n]);
}
 
/* Construct an operator to apply len zeros to a crc.  len must be a power of
   two.  If len is not a power of two, then the result is the same as for the
   largest power of two less than len.  The result for len == 0 is the same as
   for len == 1.  A version of this routine could be easily written for any
   len, but that is not needed for this application. */
static void crc32c_zeros_op(uint32_t *even, size_t len) {
    uint32_t odd[32];       /* odd-power-of-two zeros operator */
 
    /* put operator for one zero bit in odd */
    odd[0] = POLY;              /* CRC-32C polynomial */
    uint32_t row = 1;
    for (unsigned n = 1; n < 32; n++) {
        odd[n] = row;
        row <<= 1;
    }
 
    /* put operator for two zero bits in even */
    gf2_matrix_square(even, odd);
 
    /* put operator for four zero bits in odd */
    gf2_matrix_square(odd, even);
 
    /* first square will put the operator for one zero byte (eight zero bits),
       in even -- next square puts operator for two zero bytes in odd, and so
       on, until len has been rotated down to zero */
    do {
        gf2_matrix_square(even, odd);
        len >>= 1;
        if (len == 0)
            return;
        gf2_matrix_square(odd, even);
        len >>= 1;
    } while (len);
 
    /* answer ended up in odd -- copy to even */
    for (unsigned n = 0; n < 32; n++)
        even[n] = odd[n];
}
 
/* Take a length and build four lookup tables for applying the zeros operator
   for that length, byte-by-byte on the operand. */
static void crc32c_zeros(uint32_t zeros[][256], size_t len) {
    uint32_t op[32];
 
    crc32c_zeros_op(op, len);
    for (unsigned n = 0; n < 256; n++) {
        zeros[0][n] = gf2_matrix_times(op, n);
        zeros[1][n] = gf2_matrix_times(op, n << 8);
        zeros[2][n] = gf2_matrix_times(op, n << 16);
        zeros[3][n] = gf2_matrix_times(op, n << 24);
    }
}
 
/* Apply the zeros operator table to crc. */
static inline uint32_t crc32c_shift(uint32_t zeros[][256], uint32_t crc) {
    return zeros[0][crc & 0xff] ^ zeros[1][(crc >> 8) & 0xff] ^
           zeros[2][(crc >> 16) & 0xff] ^ zeros[3][crc >> 24];
}
 
/* Block sizes for three-way parallel crc computation.  LONG and SHORT must
   both be powers of two.  The associated string constants must be set
   accordingly, for use in constructing the assembler instructions. */
#define LONG 8192
#define LONGx1 "8192"
#define LONGx2 "16384"
#define SHORT 256
#define SHORTx1 "256"
#define SHORTx2 "512"
 
/* Tables for hardware crc that shift a crc by LONG and SHORT zeros. */
static pthread_once_t crc32c_once_hw = PTHREAD_ONCE_INIT;
static uint32_t crc32c_long[4][256];
static uint32_t crc32c_short[4][256];
 
/* Initialize tables for shifting crcs. */
static void crc32c_init_hw(void) {
    crc32c_zeros(crc32c_long, LONG);
    crc32c_zeros(crc32c_short, SHORT);
}
 
#endif /* __x86_64__ || CRC32C_HW_ARM */
 
#ifdef __x86_64__
 
/* Compute CRC-32C using the Intel hardware instruction. */
static uint32_t crc32c_hw(uint32_t crc, void const *buf, size_t len) {
    /* populate shift tables the first time through */
    pthread_once(&crc32c_once_hw, crc32c_init_hw);
 
    /* pre-process the crc */
    crc = ~crc;
    uint64_t crc0 = crc;            /* 64-bits for crc32q instruction */
 
    /* the crc32 instructions accept unaligned operands, so the data pointer is
       not first brought to an eight-byte boundary: doing that a byte at a time
       costs more than it saves */
    unsigned char const *next = (unsigned char const *)buf;
 
    /* compute the crc on sets of LONG*3 bytes, executing three independent crc
       instructions, each on LONG bytes -- this is optimized for the Nehalem,
       Westmere, Sandy Bridge, and Ivy Bridge architectures, which have a
       throughput of one crc per cycle, but a latency of three cycles */
    while (len >= LONG*3) {
        uint64_t crc1 = 0;
        uint64_t crc2 = 0;
        unsigned char const * const end = next + LONG;
        do {
            __asm__("crc32q\t" "(%3), %0\n\t"
                    "crc32q\t" LONGx1 "(%3), %1\n\t"
                    "crc32q\t" LONGx2 "(%3), %2"
                    : "=r"(crc0), "=r"(crc1), "=r"(crc2)
                    : "r"(next), "0"(crc0), "1"(crc1), "2"(crc2));
            next += 8;
        } while (next < end);
        crc0 = crc32c_shift(crc32c_long, crc0) ^ crc1;
        crc0 = crc32c_shift(crc32c_long, crc0) ^ crc2;
        next += LONG*2;
        len -= LONG*3;
    }
 
    /* do the same thing, but now on SHORT*3 blocks for the remaining data less
       than a LONG*3 block */
    while (len >= SHORT*3) {
        uint64_t crc1 = 0;
        uint64_t crc2 = 0;
        unsigned char const * const end = next + SHORT;
        do {
            __asm__("crc32q\t" "(%3), %0\n\t"
                    "crc32q\t" SHORTx1 "(%3), %1\n\t"
                    "crc32q\t" SHORTx2 "(%3), %2"
                    : "=r"(crc0), "=r"(crc1), "=r"(crc2)
                    : "r"(next), "0"(crc0), "1"(crc1), "2"(crc2));
            next += 8;
        } while (next < end);
        crc0 = crc32c_shift(crc32c_short, crc0) ^ crc1;
        crc0 = crc32c_shift(crc32c_short, crc0) ^ crc2;
        next += SHORT*2;
        len -= SHORT*3;
    }
 
    /* compute the crc on the remaining eight-byte units less than a SHORT*3
       block */
    {
        unsigned char const * const end = next + (len - (len & 7));
        while (next < end) {
            __asm__("crc32q\t" "(%1), %0"
                    : "=r"(crc0)
                    : "r"(next), "0"(crc0));
            next += 8;
        }
        len &= 7;
    }
 
    /* compute the crc for up to seven trailing bytes, four, two and one byte at
       a time */
    crc = crc0;
    if (len & 4) {
        __asm__("crc32l\t" "(%1), %0"
                : "=r"(crc)
                : "r"(next), "0"(crc));
        next += 4;
    }
    if (len & 2) {
        __asm__("crc32w\t" "(%1), %0"
                : "=r"(crc)
                : "r"(next), "0"(crc));
        next += 2;
    }
    if (len & 1) {
        __asm__("crc32b\t" "(%1), %0"
                : "=r"(crc)
                : "r"(next), "0"(crc));
    }
 
    /* return a post-processed crc */
    return ~crc;
}
 
/* Check for SSE 4.2.  SSE 4.2 was first supported in Nehalem processors
   introduced in November, 2008.  This does not check for the existence of the
   cpuid instruction itself, which was introduced on the 486SL in 1992, so this
   will fail on earlier x86 processors.  cpuid works on all Pentium and later
   processors. */
#define SSE42(have) \
    do { \
        uint32_t eax, ecx; \
        eax = 1; \
        __asm__("cpuid" \
                : "=c"(ecx) \
                : "a"(eax) \
                : "%ebx", "%edx"); \
        (have) = (ecx >> 20) & 1; \
    } while (0)
 
static int crc32c_have_sse42(void) {
    int sse42;
 
    SSE42(sse42);
    return sse42;
}
 
/* Compute a CRC-32C.  If the crc32 instruction is available, use the hardware
   version.  Otherwise, use the software version.  The cpuid check is done only
   once, since cpuid is slow and traps to the hypervisor in virtual machines. */
uint32_t crc32c(uint32_t crc, void const *buf, size_t len) {
    static const int sse42 = crc32c_have_sse42();
 
    return sse42 ? crc32c_hw(crc, buf, len) : crc32c_sw(crc, buf, len);
}
 
#elif defined(CRC32C_HW_ARM)
 
/* The crc32c intrinsics need the crc extension: enable it for the functions
   using them when it is not enabled for the whole compilation unit. */
#ifdef __ARM_FEATURE_CRC32
#define CRC32C_TARGET_CRC
#else
#define CRC32C_TARGET_CRC __attribute__((target("+crc")))
#endif
 
/* Load eight bytes from a possibly unaligned address. */
static inline uint64_t crc32c_load64(unsigned char const *next) {
    uint64_t word;
    memcpy(&word, next, sizeof(word));
    return word;
}
 
/* Compute CRC-32C using the ARMv8 hardware instructions, with the same
   three-way parallel scheme as on x86 to hide the latency of the crc32cx
   instruction. */
CRC32C_TARGET_CRC
static uint32_t crc32c_hw(uint32_t crc, void const *buf, size_t len) {
    /* populate shift tables the first time through */
    pthread_once(&crc32c_once_hw, crc32c_init_hw);
 
    uint32_t crc0 = ~crc;
    unsigned char const *next = (unsigned char const *)buf;
 
    /* compute the crc on sets of LONG*3 bytes, three independent crc
       instructions at a time, each on LONG bytes */
    while (len >= LONG*3) {
        uint32_t crc1 = 0;
        uint32_t crc2 = 0;
        unsigned char const * const end = next + LONG;
        do {
            crc0 = __crc32cd(crc0, crc32c_load64(next));
            crc1 = __crc32cd(crc1, crc32c_load64(next + LONG));
            crc2 = __crc32cd(crc2, crc32c_load64(next + LONG*2));
            next += 8;
        } while (next < end);
        crc0 = crc32c_shift(crc32c_long, crc0) ^ crc1;
        crc0 = crc32c_shift(crc32c_long, crc0) ^ crc2;
        next += LONG*2;
        len -= LONG*3;
    }
 
    /* do the same thing, but now on SHORT*3 blocks for the remaining data less
       than a LONG*3 block */
    while (len >= SHORT*3) {
        uint32_t crc1 = 0;
        uint32_t crc2 = 0;
        unsigned char const * const end = next + SHORT;
        do {
            crc0 = __crc32cd(crc0, crc32c_load64(next));
            crc1 = __crc32cd(crc1, crc32c_load64(next + SHORT));
            crc2 = __crc32cd(crc2, crc32c_load64(next + SHORT*2));
            next += 8;
        } while (next < end);
        crc0 = crc32c_shift(crc32c_short, crc0) ^ crc1;
        crc0 = crc32c_shift(crc32c_short, crc0) ^ crc2;
        next += SHORT*2;
        len -= SHORT*3;
    }
 
    /* compute the crc on the remaining eight-byte units less than a SHORT*3
       block */
    while (len >= 8) {
        crc0 = __crc32cd(crc0, crc32c_load64(next));
        next += 8;
        len -= 8;
    }
 
    /* compute the crc for up to seven trailing bytes, four, two and one byte at
       a time */
    if (len & 4) {
        uint32_t word;
        memcpy(&word, next, sizeof(word));
        crc0 = __crc32cw(crc0, word);
        next += 4;
    }
    if (len & 2) {
        uint16_t half;
        memcpy(&half, next, sizeof(half));
        crc0 = __crc32ch(crc0, half);
        next += 2;
    }
    if (len & 1)
        crc0 = __crc32cb(crc0, *next);
 
    /* return a post-processed crc */
    return ~crc0;
}
 
/* Check for the crc32c instructions: known at compile time when the crc
   extension is enabled, otherwise asked from the kernel on Linux. */
static int crc32c_have_crc(void) {
#if defined(__ARM_FEATURE_CRC32)
    return 1;
#elif defined(__linux__) && defined(HWCAP_CRC32)
    return (getauxval(AT_HWCAP) & HWCAP_CRC32) != 0;
#else
    return 0;
#endif
}
 
/* Compute a CRC-32C.  If the crc32c instructions are available, use the
   hardware version.  Otherwise, use the software version.  The check is done
   only once. */
uint32_t crc32c(uint32_t crc, void const *buf, size_t len) {
    static const int hw = crc32c_have_crc();
 
    return hw ? crc32c_hw(crc, buf, len) : crc32c_sw(crc, buf, len);
}
 
#else /* !__x86_64__ && !CRC32C_HW_ARM */
 
uint32_t crc32c(uint32_t crc, void const *buf, size_t len) {
    return crc32c_sw(crc, buf, len);
}
 
#endif
 
/* Construct table for software CRC-32C little-endian calculation. */
static pthread_once_t crc32c_once_little = PTHREAD_ONCE_INIT;
static uint32_t crc32c_table_little[8][256];
static void crc32c_init_sw_little(void) {
    for (unsigned n = 0; n < 256; n++) {
        uint32_t crc = n;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc32c_table_little[0][n] = crc;
    }
    for (unsigned n = 0; n < 256; n++) {
        uint32_t crc = crc32c_table_little[0][n];
        for (unsigned k = 1; k < 8; k++) {
            crc = crc32c_table_little[0][crc & 0xff] ^ (crc >> 8);
            crc32c_table_little[k][n] = crc;
        }
    }
}
 
/* Compute a CRC-32C in software assuming a little-endian architecture,
   constructing the required table if that hasn't already been done. */
uint32_t crc32c_sw_little(uint32_t crc, void const *buf, size_t len) {
    unsigned char const *next = (unsigned char const *)buf;
 
    pthread_once(&crc32c_once_little, crc32c_init_sw_little);
    crc = ~crc;
    while (len && ((uintptr_t)next & 7) != 0) {
        crc = crc32c_table_little[0][(crc ^ *next++) & 0xff] ^ (crc >> 8);
        len--;
    }
    if (len >= 8) {
        uint64_t crcw = crc;
        do {
            crcw ^= *(uint64_t const *)next;
            crcw = crc32c_table_little[7][crcw & 0xff] ^
                   crc32c_table_little[6][(crcw >> 8) & 0xff] ^
                   crc32c_table_little[5][(crcw >> 16) & 0xff] ^
                   crc32c_table_little[4][(crcw >> 24) & 0xff] ^
                   crc32c_table_little[3][(crcw >> 32) & 0xff] ^
                   crc32c_table_little[2][(crcw >> 40) & 0xff] ^
                   crc32c_table_little[1][(crcw >> 48) & 0xff] ^
                   crc32c_table_little[0][crcw >> 56];
            next += 8;
            len -= 8;
        } while (len >= 8);
        crc = crcw;
    }
    while (len) {
        crc = crc32c_table_little[0][(crc ^ *next++) & 0xff] ^ (crc >> 8);
        len--;
    }
    return ~crc;
}
 
/* Swap the bytes in a uint64_t.  (Only for big-endian.) */
#if defined(__has_builtin) || (defined(__GNUC__) && \
    (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 3)))
#  define swap __builtin_bswap64
#else
static inline uint64_t swap(uint64_t x) {
    x = ((x << 8) & 0xff00ff00ff00ff00) | ((x >> 8) & 0xff00ff00ff00ff);
    x = ((x << 16) & 0xffff0000ffff0000) | ((x >> 16) & 0xffff0000ffff);
    return (x << 32) | (x >> 32);
}
#endif
 
/* Construct tables for software CRC-32C big-endian calculation. */
static pthread_once_t crc32c_once_big = PTHREAD_ONCE_INIT;
static uint32_t crc32c_table_big_byte[256];
static uint64_t crc32c_table_big[8][256];
static void crc32c_init_sw_big(void) {
    for (unsigned n = 0; n < 256; n++) {
        uint32_t crc = n;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc = crc & 1 ? (crc >> 1) ^ POLY : crc >> 1;
        crc32c_table_big_byte[n] = crc;
    }
    for (unsigned n = 0; n < 256; n++) {
        uint32_t crc = crc32c_table_big_byte[n];
        crc32c_table_big[0][n] = swap(crc);
        for (unsigned k = 1; k < 8; k++) {
            crc = crc32c_table_big_byte[crc & 0xff] ^ (crc >> 8);
            crc32c_table_big[k][n] = swap(crc);
        }
    }
}
 
/* Compute a CRC-32C in software assuming a big-endian architecture,
   constructing the required tables if that hasn't already been done. */
uint32_t crc32c_sw_big(uint32_t crc, void const *buf, size_t len) {
    unsigned char const *next = (unsigned char const *)buf;
 
    pthread_once(&crc32c_once_big, crc32c_init_sw_big);
    crc = ~crc;
    while (len && ((uintptr_t)next & 7) != 0) {
        crc = crc32c_table_big_byte[(crc ^ *next++) & 0xff] ^ (crc >> 8);
        len--;
    }
    if (len >= 8) {
        uint64_t crcw = swap(crc);
        do {
            crcw ^= *(uint64_t const *)next;
            crcw = crc32c_table_big[0][crcw & 0xff] ^
                   crc32c_table_big[1][(crcw >> 8) & 0xff] ^
                   crc32c_table_big[2][(crcw >> 16) & 0xff] ^
                   crc32c_table_big[3][(crcw >> 24) & 0xff] ^
                   crc32c_table_big[4][(crcw >> 32) & 0xff] ^
                   crc32c_table_big[5][(crcw >> 40) & 0xff] ^
                   crc32c_table_big[6][(crcw >> 48) & 0xff] ^
                   crc32c_table_big[7][(crcw >> 56)];
            next += 8;
            len -= 8;
        } while (len >= 8);
        crc = swap(crcw);
    }
    while (len) {
        crc = crc32c_table_big_byte[(crc ^ *next++) & 0xff] ^ (crc >> 8);
        len--;
    }
    return ~crc;
}
 
/* Table-driven software CRC-32C.  This is about 15 times slower than using the
   hardware instructions.  Determine the endianess of the processor and proceed
   accordingly.  Ideally the endianess will be determined at compile time, in
   which case the unused functions and tables for the other endianess will be
   removed by the optimizer.  If not, then the proper routines and tables will
   be used, even if the endianess is changed mid-stream.  (Yes, there are
   processors that permit that -- go figure.) */
uint32_t crc32c_sw(uint32_t crc, void const *buf, size_t len) {
    static int const little = 1;
    if (*(char const *)&little)
        return crc32c_sw_little(crc, buf, len);
    else
        return crc32c_sw_big(crc, buf, len);
}
 
/******************************************************************************/
/*                        c r c 3 2 c _ c o m b i n e                         */
/******************************************************************************/

uint32_t crc32c_combine(uint32_t crcA, uint32_t crcB, size_t lengthB)
{
  // Original source from Stephan Brumme
  // https://github.com/stbrumme/crc32
  // https://create.stephan-brumme.com/crc32/
  // based on Mark Adler's crc_combine from
  // https://github.com/madler/pigz/blob/master/pigz.c
 
  // main idea:
  // - if you have two equally-sized blocks A and B,
  //   then you can create a block C = A ^ B
  //   which has the property crc(C) = crc(A) ^ crc(B)
  // - if you append length(B) zeros to A and call it A' (think of it as AAAA000)
  //   and   prepend length(A) zeros to B and call it B' (think of it as 0000BBB)
  //   then exists a C' = A' ^ B'
  // - remember: if you XOR someting with zero, it remains unchanged: X ^ 0 = X
  // - that means C' = A concat B so that crc(A concat B) = crc(C') = crc(A') ^ crc(B')
  // - the trick is to compute crc(A') based on crc(A)
  //                       and crc(B') based on crc(B)
  // - since B' starts with many zeros, the crc of those initial zeros is still zero
  // - that means crc(B') = crc(B)
  // - unfortunately the trailing zeros of A' change the crc, so usually crc(A') != crc(A)
  // - the following code is a fast algorithm to compute crc(A')
  // - starting with crc(A) and appending length(B) zeros, needing just log2(length(B)) iterations
  // - the details are explained by the original author at
  //   https://stackoverflow.com/questions/23122312/crc-calculation-of-a-mostly-static-data-stream/23126768
  //
  // notes:
  // - I squeezed everything into one function to keep global namespace clean (original code two helper functions)
  // - most original comments are still in place, I added comments where these helper functions where made inline code
  // - performance-wise there isn't any differenze to the original zlib/pigz code
 
  // degenerated case
  if (lengthB == 0)
    return crcA;

  const uint32_t CrcBits = 32;
 
  uint32_t odd [CrcBits]; // odd-power-of-two  zeros operator
  uint32_t even[CrcBits]; // even-power-of-two zeros operator
 
  // put operator for one zero bit in odd
  odd[0] = POLY;          // CRC-32 polynomial
  for (int i = 1; i < (int)CrcBits; i++)
    odd[i] = 1 << (i - 1);
 
  // put operator for two zero bits in even
  // same as gf2_matrix_square(even, odd);
  for (int i = 0; i < (int)CrcBits; i++)
  {
    uint32_t vec = odd[i];
    even[i] = 0;
    for (int j = 0; vec != 0; j++, vec >>= 1)
      if (vec & 1)
        even[i] ^= odd[j];
  }
  // put operator for four zero bits in odd
  // same as gf2_matrix_square(odd, even);
  for (int i = 0; i < (int)CrcBits; i++)
  {
    uint32_t vec = even[i];
    odd[i] = 0;
    for (int j = 0; vec != 0; j++, vec >>= 1)
      if (vec & 1)
        odd[i] ^= even[j];
  }
 
  // the following loop becomes much shorter if I keep swapping even and odd
  uint32_t* a = even;
  uint32_t* b = odd;
  // apply secondLength zeros to firstCrc32
  for (; lengthB > 0; lengthB >>= 1)
  {
    // same as gf2_matrix_square(a, b);
    for (int i = 0; i < (int)CrcBits; i++)
    {
      uint32_t vec = b[i];
      a[i] = 0;
      for (int j = 0; vec != 0; j++, vec >>= 1)
        if (vec & 1)
          a[i] ^= b[j];
    }
 
    // apply zeros operator for this bit
    if (lengthB & 1)
    {
      // same as firstCrc32 = gf2_matrix_times(a, firstCrc32);
      uint32_t sum = 0;
      for (int i = 0; crcA != 0; i++, crcA >>= 1)
        if (crcA & 1)
          sum ^= a[i];
      crcA = sum;
    }
 
    // switch even and odd
    uint32_t* t = a; a = b; b = t;
  }
 
  // return combined crc
  return crcA ^ crcB;
}
 
/******************************************************************************/
/*                                  m a i n                                   */
/******************************************************************************/
 
#ifdef TEST
 
#include <cstdio>
#include <cstdlib>
#include "load.h"
 
int main(int argc, char **argv) {
    /* interpret argument */
    if (argc > 2) {
        fputs("only one argument permitted\n", stderr);
        return 1;
    }
    long rep = argc == 1 ? 1 : strtol(argv[1], NULL, 10);
    if (rep == 0) {
        fputs("usage: crc32c [[-]nnn] < data\n"
              "  where nnn is the number of times to repeat\n"
              "  negative forces the software implementation\n", stderr);
        return 0;
    }
    int sw = rep < 0;
    if (sw)
        rep = -rep;
 
    /* read input */
    void *dat = NULL;
    size_t size, len;
    int ret = load(stdin, 0, &dat, &size, &len);
    if (ret) {
        fputs("error reading from stdin\n", stderr);
        return 1;
    }
 
    /* compute crc rep times */
    uint32_t crc;
    do {
        crc = sw ? crc32c_sw(0, dat, len) : crc32c(0, dat, len);
    } while (--rep);
 
    /* show crc */
    printf("0x%08x\n", crc);
 
    /* clean up */
    free(dat);
    return 0;
}
 
#endif /* TEST */