folly/detail/ChecksumDetail.cpp

   1 /*
   2  * crc32_impl.h
   3  *
   4  * Copyright 2016 Eric Biggers
   5  *
   6  * Permission is hereby granted, free of charge, to any person
   7  * obtaining a copy of this software and associated documentation
   8  * files (the "Software"), to deal in the Software without
   9  * restriction, including without limitation the rights to use,
  10  * copy, modify, merge, publish, distribute, sublicense, and/or sell
  11  * copies of the Software, and to permit persons to whom the
  12  * Software is furnished to do so, subject to the following
  13  * conditions:
  14  *
  15  * The above copyright notice and this permission notice shall be
  16  * included in all copies or substantial portions of the Software.
  17  *
  18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  19  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  20  * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  21  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  22  * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  23  * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  24  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  25  * OTHER DEALINGS IN THE SOFTWARE.
  26  */
  27
  28 /*
  29  * CRC-32 folding with PCLMULQDQ.
  30  *
  31  * The basic idea is to repeatedly "fold" each 512 bits into the next
  32  * 512 bits, producing an abbreviated message which is congruent the
  33  * original message modulo the generator polynomial G(x).
  34  *
  35  * Folding each 512 bits is implemented as eight 64-bit folds, each of
  36  * which uses one carryless multiplication instruction.  It's expected
  37  * that CPUs may be able to execute some of these multiplications in
  38  * parallel.
  39  *
  40  * Explanation of "folding": let A(x) be 64 bits from the message, and
  41  * let B(x) be 95 bits from a constant distance D later in the
  42  * message.  The relevant portion of the message can be written as:
  43  *
  44  *      M(x) = A(x)*x^D + B(x)
  45  *
  46  * ... where + and * represent addition and multiplication,
  47  * respectively, of polynomials over GF(2).  Note that when
  48  * implemented on a computer, these operations are equivalent to XOR
  49  * and carryless multiplication, respectively.
  50  *
  51  * For the purpose of CRC calculation, only the remainder modulo the
  52  * generator polynomial G(x) matters:
  53  *
  54  * M(x) mod G(x) = (A(x)*x^D + B(x)) mod G(x)
  55  *
  56  * Since the modulo operation can be applied anywhere in a sequence of
  57  * additions and multiplications without affecting the result, this is
  58  * equivalent to:
  59  *
  60  * M(x) mod G(x) = (A(x)*(x^D mod G(x)) + B(x)) mod G(x)
  61  *
  62  * For any D, 'x^D mod G(x)' will be a polynomial with maximum degree
  63  * 31, i.e.  a 32-bit quantity.  So 'A(x) * (x^D mod G(x))' is
  64  * equivalent to a carryless multiplication of a 64-bit quantity by a
  65  * 32-bit quantity, producing a 95-bit product.  Then, adding
  66  * (XOR-ing) the product to B(x) produces a polynomial with the same
  67  * length as B(x) but with the same remainder as 'A(x)*x^D + B(x)'.
  68  * This is the basic fold operation with 64 bits.
  69  *
  70  * Note that the carryless multiplication instruction PCLMULQDQ
  71  * actually takes two 64-bit inputs and produces a 127-bit product in
  72  * the low-order bits of a 128-bit XMM register.  This works fine, but
  73  * care must be taken to account for "bit endianness".  With the CRC
  74  * version implemented here, bits are always ordered such that the
  75  * lowest-order bit represents the coefficient of highest power of x
  76  * and the highest-order bit represents the coefficient of the lowest
  77  * power of x.  This is backwards from the more intuitive order.
  78  * Still, carryless multiplication works essentially the same either
  79  * way.  It just must be accounted for that when we XOR the 95-bit
  80  * product in the low-order 95 bits of a 128-bit XMM register into
  81  * 128-bits of later data held in another XMM register, we'll really
  82  * be XOR-ing the product into the mathematically higher degree end of
  83  * those later bits, not the lower degree end as may be expected.
  84  *
  85  * So given that caveat and the fact that we process 512 bits per
  86  * iteration, the 'D' values we need for the two 64-bit halves of each
  87  * 128 bits of data are:
  88  *
  89  * D = (512 + 95) - 64 for the higher-degree half of each 128
  90  *                 bits, i.e. the lower order bits in
  91  *                 the XMM register
  92  *
  93  *    D = (512 + 95) - 128 for the lower-degree half of each 128
  94  *                 bits, i.e. the higher order bits in
  95  *                 the XMM register
  96  *
  97  * The required 'x^D mod G(x)' values were precomputed.
  98  *
  99  * When <= 512 bits remain in the message, we finish up by folding
 100  * across smaller distances.  This works similarly; the distance D is
 101  * just different, so different constant multipliers must be used.
 102  * Finally, once the remaining message is just 64 bits, it is is
 103  * reduced to the CRC-32 using Barrett reduction (explained later).
 104  *
 105  * For more information see the original paper from Intel: "Fast CRC
 106  *    Computation for Generic Polynomials Using PCLMULQDQ
 107  *    Instruction" December 2009
 108  *    http://www.intel.com/content/dam/www/public/us/en/documents/
 109  *    white-papers/
 110  *    fast-crc-computation-generic-polynomials-pclmulqdq-paper.pdf
 111  */
 112
 113 #include <folly/detail/ChecksumDetail.h>
 114
 115 namespace folly {
 116 namespace detail {
 117
 118 uint32_t
 119 crc32_hw_aligned(uint32_t remainder, const __m128i* p, size_t vec_count) {
 120   /* Constants precomputed by gen_crc32_multipliers.c.  Do not edit! */
 121   const __m128i multipliers_4 = _mm_set_epi32(0, 0x1D9513D7, 0, 0x8F352D95);
 122   const __m128i multipliers_2 = _mm_set_epi32(0, 0x81256527, 0, 0xF1DA05AA);
 123   const __m128i multipliers_1 = _mm_set_epi32(0, 0xCCAA009E, 0, 0xAE689191);
 124   const __m128i final_multiplier = _mm_set_epi32(0, 0, 0, 0xB8BC6765);
 125   const __m128i mask32 = _mm_set_epi32(0, 0, 0, 0xFFFFFFFF);
 126   const __m128i barrett_reduction_constants =
 127       _mm_set_epi32(0x1, 0xDB710641, 0x1, 0xF7011641);
 128
 129   const __m128i* const end = p + vec_count;
 130   const __m128i* const end512 = p + (vec_count & ~3);
 131   __m128i x0, x1, x2, x3;
 132
 133   /*
 134    * Account for the current 'remainder', i.e. the CRC of the part of
 135    * the message already processed.  Explanation: rewrite the message
 136    * polynomial M(x) in terms of the first part A(x), the second part
 137    * B(x), and the length of the second part in bits |B(x)| >= 32:
 138    *
 139    *    M(x) = A(x)*x^|B(x)| + B(x)
 140    *
 141    * Then the CRC of M(x) is:
 142    *
 143    *    CRC(M(x)) = CRC(A(x)*x^|B(x)| + B(x))
 144    *              = CRC(A(x)*x^32*x^(|B(x)| - 32) + B(x))
 145    *              = CRC(CRC(A(x))*x^(|B(x)| - 32) + B(x))
 146    *
 147    * Note: all arithmetic is modulo G(x), the generator polynomial; that's
 148    * why A(x)*x^32 can be replaced with CRC(A(x)) = A(x)*x^32 mod G(x).
 149    *
 150    * So the CRC of the full message is the CRC of the second part of the
 151    * message where the first 32 bits of the second part of the message
 152    * have been XOR'ed with the CRC of the first part of the message.
 153    */
 154   x0 = *p++;
 155   x0 ^= _mm_set_epi32(0, 0, 0, remainder);
 156
 157   if (p > end512) /* only 128, 256, or 384 bits of input? */
 158     goto _128_bits_at_a_time;
 159   x1 = *p++;
 160   x2 = *p++;
 161   x3 = *p++;
 162
 163   /* Fold 512 bits at a time */
 164   for (; p != end512; p += 4) {
 165     __m128i y0, y1, y2, y3;
 166
 167     y0 = p[0];
 168     y1 = p[1];
 169     y2 = p[2];
 170     y3 = p[3];
 171
 172     /*
 173      * Note: the immediate constant for PCLMULQDQ specifies which
 174      * 64-bit halves of the 128-bit vectors to multiply:
 175      *
 176      * 0x00 means low halves (higher degree polynomial terms for us)
 177      * 0x11 means high halves (lower degree polynomial terms for us)
 178      */
 179     y0 ^= _mm_clmulepi64_si128(x0, multipliers_4, 0x00);
 180     y1 ^= _mm_clmulepi64_si128(x1, multipliers_4, 0x00);
 181     y2 ^= _mm_clmulepi64_si128(x2, multipliers_4, 0x00);
 182     y3 ^= _mm_clmulepi64_si128(x3, multipliers_4, 0x00);
 183     y0 ^= _mm_clmulepi64_si128(x0, multipliers_4, 0x11);
 184     y1 ^= _mm_clmulepi64_si128(x1, multipliers_4, 0x11);
 185     y2 ^= _mm_clmulepi64_si128(x2, multipliers_4, 0x11);
 186     y3 ^= _mm_clmulepi64_si128(x3, multipliers_4, 0x11);
 187
 188     x0 = y0;
 189     x1 = y1;
 190     x2 = y2;
 191     x3 = y3;
 192   }
 193
 194   /* Fold 512 bits => 128 bits */
 195   x2 ^= _mm_clmulepi64_si128(x0, multipliers_2, 0x00);
 196   x3 ^= _mm_clmulepi64_si128(x1, multipliers_2, 0x00);
 197   x2 ^= _mm_clmulepi64_si128(x0, multipliers_2, 0x11);
 198   x3 ^= _mm_clmulepi64_si128(x1, multipliers_2, 0x11);
 199   x3 ^= _mm_clmulepi64_si128(x2, multipliers_1, 0x00);
 200   x3 ^= _mm_clmulepi64_si128(x2, multipliers_1, 0x11);
 201   x0 = x3;
 202
 203 _128_bits_at_a_time:
 204   while (p != end) {
 205     /* Fold 128 bits into next 128 bits */
 206     x1 = *p++;
 207     x1 ^= _mm_clmulepi64_si128(x0, multipliers_1, 0x00);
 208     x1 ^= _mm_clmulepi64_si128(x0, multipliers_1, 0x11);
 209     x0 = x1;
 210   }
 211
 212   /* Now there are just 128 bits left, stored in 'x0'. */
 213
 214   /*
 215    * Fold 128 => 96 bits.  This also implicitly appends 32 zero bits,
 216    * which is equivalent to multiplying by x^32.  This is needed because
 217    * the CRC is defined as M(x)*x^32 mod G(x), not just M(x) mod G(x).
 218    */
 219   x0 = _mm_srli_si128(x0, 8) ^ _mm_clmulepi64_si128(x0, multipliers_1, 0x10);
 220
 221   /* Fold 96 => 64 bits */
 222   x0 = _mm_srli_si128(x0, 4) ^
 223       _mm_clmulepi64_si128(x0 & mask32, final_multiplier, 0x00);
 224
 225   /*
 226    * Finally, reduce 64 => 32 bits using Barrett reduction.
 227    *
 228    * Let M(x) = A(x)*x^32 + B(x) be the remaining message.  The goal is to
 229    * compute R(x) = M(x) mod G(x).  Since degree(B(x)) < degree(G(x)):
 230    *
 231    *    R(x) = (A(x)*x^32 + B(x)) mod G(x)
 232    *         = (A(x)*x^32) mod G(x) + B(x)
 233    *
 234    * Then, by the Division Algorithm there exists a unique q(x) such that:
 235    *
 236    *    A(x)*x^32 mod G(x) = A(x)*x^32 - q(x)*G(x)
 237    *
 238    * Since the left-hand side is of maximum degree 31, the right-hand side
 239    * must be too.  This implies that we can apply 'mod x^32' to the
 240    * right-hand side without changing its value:
 241    *
 242    *    (A(x)*x^32 - q(x)*G(x)) mod x^32 = q(x)*G(x) mod x^32
 243    *
 244    * Note that '+' is equivalent to '-' in polynomials over GF(2).
 245    *
 246    * We also know that:
 247    *
 248    *                  / A(x)*x^32 \
 249    *    q(x) = floor (  ---------  )
 250    *                  \    G(x)   /
 251    *
 252    * To compute this efficiently, we can multiply the top and bottom by
 253    * x^32 and move the division by G(x) to the top:
 254    *
 255    *                  / A(x) * floor(x^64 / G(x)) \
 256    *    q(x) = floor (  -------------------------  )
 257    *                  \           x^32            /
 258    *
 259    * Note that floor(x^64 / G(x)) is a constant.
 260    *
 261    * So finally we have:
 262    *
 263    *                              / A(x) * floor(x^64 / G(x)) \
 264    *    R(x) = B(x) + G(x)*floor (  -------------------------  )
 265    *                              \           x^32            /
 266    */
 267   x1 = x0;
 268   x0 = _mm_clmulepi64_si128(x0 & mask32, barrett_reduction_constants, 0x00);
 269   x0 = _mm_clmulepi64_si128(x0 & mask32, barrett_reduction_constants, 0x10);
 270   return _mm_cvtsi128_si32(_mm_srli_si128(x0 ^ x1, 4));
 271 }
 272 }
 273 } // namespace