Adding additional finite field and carryless multiplication ops

[johnplatts]

There are some additional finite field and carryless multiplication ops that have actual use cases, including cryptography (including SM4, Camellia, Grain128-AEADv2, and post-quantum crytography), hashing (including Groestl hash, Toeplitz hash, and CRC checksums), and Reed-Solomon error correction.

Here are the proposed additional finite field and carryless multiplication ops:

• VU8 CLMul(VU8 a, VU8 b) - does carryless multiplication of a[i] and b[i], returning the lower 8 bits without any modular reduction, equivalent to NEON vmul_p8 or SVE2 svpmul_u8
• VU16 WidenCLMul(DU16, VU8 a, VU8 b) - does carryless multiplication of a[i] and b[i], returned as a U16, equivalent to NEON vmull_p8
• VU8 CLMulHigh(VU8 a, VU8 b) - does carryless multiplication of a[i] and b[i], returning the upper 8 bits, equivalent to DemoteTo(du8, ShiftRight<8>(WidenCLMulLower(du16, a, b)))
• VU8 GF2P8Mul(VU8 a, VU8 b) - does carryless multiplication of a[i] and b[i] followed by finite field modular reduction by x⁸ + x⁴+ x³ + x + 1, equivalent to GFNI _mm_gf2p8mul_epi8
• VU16 CLMulEven(DU16, VU8 a, VU8 b) - does carryless multiplication of a[2*i] and b[2*i] - equivalent to SVE2 svpmullb_u16
• VU16 CLMulOdd(DU16, VU8 a, VU8 b) - does carryless multiplication of a[2*i+1] and b[2*i+1] - equivalent to SVE2 svpmullt_u16
• VU16 CLWidenMulPairwiseAdd(DU16, VU8 a, VU8 b) - equivalent to Xor(CLMulEven(du16, a, b), CLMulOdd(du16, a, b)) or PPC8 vec_pmsum_be(a, b)
• VU8 GF2P8MulInverse(VU8 v) - computes the GF(2^8) multiplicative inverse of v[i] (modulo x⁸ + x⁴+ x³ + x + 1) - equivalent to GFNI _mm_gf2p8affineinv_epi64_epi8(v, _mm_set1_epi64x(0x0102040810204080), 0x00)
• VU8 GaloisAffine<uint64_k kMatrix, uint8_t kXorMask>(VU8 v) - applies affine transform kMatrix to v[i] - equivalent to GFNI _mm_gf2p8affine_epi64_epi8(v, _mm_set1_epi64x(static_cast<int64_t>(kMatrix)), kXorMask)
• VU8 GaloisAffineInv<uint64_k kMatrix, uint8_t kXorMask>(VU8 v) - equivalent to GaloisAffine<kMatrix, kXorMask>(GF2P8MulInverse(v)) or GFNI _mm_gf2p8affineinv_epi64_epi8(v, _mm_set1_epi64x(static_cast<int64_t>(kMatrix)), kXorMask)

[johnplatts]

Here is how GF2P8MulInverse can be implemented on targets other than HWY_SCALAR:

template <class V, HWY_IF_U8(TFromV<V>)>
HWY_API V GF2P8MulInverse(V v) {
  using D = DFromV<decltype(v)>;

  const D d;

#if HWY_TARGET == HWY_RVV
  // Need at least LMUL=1 on RVV to ensure that Lanes(d_mul_inv) is at least 16
  const ScalableTag<uint8_t, HWY_MAX(HWY_POW2_D(D), 0)> d_mul_inv;
#else
  const FixedTag<uint8_t, HWY_MAX(HWY_MAX_LANES_D(D), 16)> d_mul_inv;
#endif

  const auto aes_enc_result =
      AESLastRound(ResizeBitCast(d_mul_inv, v), Zero(d_mul_inv));
  const auto inv_affine_result =
      Xor(Xor3(RotateRight<7>(aes_enc_result), RotateRight<5>(aes_enc_result),
               RotateRight<2>(aes_enc_result)),
          Set(d_mul_inv, 0x05));
  const auto mul_inverse_result = TableLookupBytes(
      inv_affine_result, Dup128VecFromValues(d_mul_inv, 0, 13, 10, 7, 4, 1, 14,
                                             11, 8, 5, 2, 15, 12, 9, 6, 3));

  return ResizeBitCast(d, mul_inverse_result);
}

[johnplatts]

Here is how GaloisAffine can be implemented:

template <uint64_t kMatrix, uint8_t kXorMask, class V, HWY_IF_U8(TFromV<V>)>
HWY_API V GaloisAffine(V v) {
  constexpr uint8_t kMask0 =
      static_cast<uint8_t>(((kMatrix >> 56) & 0x01) | ((kMatrix >> 47) & 0x02) |
                           ((kMatrix >> 38) & 0x04) | ((kMatrix >> 29) & 0x08) |
                           ((kMatrix >> 20) & 0x10) | ((kMatrix >> 11) & 0x20) |
                           ((kMatrix >> 2) & 0x40) | ((kMatrix << 7) & 0x80));
  constexpr uint8_t kMask1 =
      static_cast<uint8_t>(((kMatrix >> 57) & 0x01) | ((kMatrix >> 48) & 0x02) |
                           ((kMatrix >> 39) & 0x04) | ((kMatrix >> 30) & 0x08) |
                           ((kMatrix >> 21) & 0x10) | ((kMatrix >> 12) & 0x20) |
                           ((kMatrix >> 3) & 0x40) | ((kMatrix << 6) & 0x80));
  constexpr uint8_t kMask2 =
      static_cast<uint8_t>(((kMatrix >> 58) & 0x01) | ((kMatrix >> 49) & 0x02) |
                           ((kMatrix >> 40) & 0x04) | ((kMatrix >> 31) & 0x08) |
                           ((kMatrix >> 22) & 0x10) | ((kMatrix >> 13) & 0x20) |
                           ((kMatrix >> 4) & 0x40) | ((kMatrix << 5) & 0x80));
  constexpr uint8_t kMask3 =
      static_cast<uint8_t>(((kMatrix >> 59) & 0x01) | ((kMatrix >> 50) & 0x02) |
                           ((kMatrix >> 41) & 0x04) | ((kMatrix >> 32) & 0x08) |
                           ((kMatrix >> 23) & 0x10) | ((kMatrix >> 14) & 0x20) |
                           ((kMatrix >> 5) & 0x40) | ((kMatrix << 4) & 0x80));
  constexpr uint8_t kMask4 =
      static_cast<uint8_t>(((kMatrix >> 60) & 0x01) | ((kMatrix >> 51) & 0x02) |
                           ((kMatrix >> 42) & 0x04) | ((kMatrix >> 33) & 0x08) |
                           ((kMatrix >> 24) & 0x10) | ((kMatrix >> 15) & 0x20) |
                           ((kMatrix >> 6) & 0x40) | ((kMatrix << 3) & 0x80));
  constexpr uint8_t kMask5 =
      static_cast<uint8_t>(((kMatrix >> 61) & 0x01) | ((kMatrix >> 52) & 0x02) |
                           ((kMatrix >> 43) & 0x04) | ((kMatrix >> 34) & 0x08) |
                           ((kMatrix >> 25) & 0x10) | ((kMatrix >> 16) & 0x20) |
                           ((kMatrix >> 7) & 0x40) | ((kMatrix << 2) & 0x80));
  constexpr uint8_t kMask6 =
      static_cast<uint8_t>(((kMatrix >> 62) & 0x01) | ((kMatrix >> 53) & 0x02) |
                           ((kMatrix >> 44) & 0x04) | ((kMatrix >> 35) & 0x08) |
                           ((kMatrix >> 26) & 0x10) | ((kMatrix >> 17) & 0x20) |
                           ((kMatrix >> 8) & 0x40) | ((kMatrix << 1) & 0x80));
  constexpr uint8_t kMask7 =
      static_cast<uint8_t>(((kMatrix >> 63) & 0x01) | ((kMatrix >> 54) & 0x02) |
                           ((kMatrix >> 45) & 0x04) | ((kMatrix >> 36) & 0x08) |
                           ((kMatrix >> 27) & 0x10) | ((kMatrix >> 18) & 0x20) |
                           ((kMatrix >> 9) & 0x40) | (kMatrix & 0x80));

  const DFromV<decltype(v)> d;
  const RebindToSigned<decltype(d)> di;

  auto result0 = Xor(IfThenElseZero(TestBit(v, Set(d, 0x01)), Set(d, kMask0)),
                     Set(d, kXorMask));
  auto result1 = IfThenElseZero(TestBit(v, Set(d, 0x02)), Set(d, kMask1));
  auto result2 = IfThenElseZero(TestBit(v, Set(d, 0x04)), Set(d, kMask2));
  auto result3 = IfThenElseZero(TestBit(v, Set(d, 0x08)), Set(d, kMask3));
  auto result4 = IfThenElseZero(TestBit(v, Set(d, 0x10)), Set(d, kMask4));
  auto result5 = IfThenElseZero(TestBit(v, Set(d, 0x20)), Set(d, kMask5));
  auto result6 = IfThenElseZero(TestBit(v, Set(d, 0x40)), Set(d, kMask6));
  auto result7 =
      BitCast(d, IfNegativeThenElseZero(BitCast(di, v),
                                        Set(di, static_cast<int8_t>(kMask7))));

  result0 = Xor(result0, result1);
  result2 = Xor(result2, result3);
  result4 = Xor(result4, result5);
  result6 = Xor(result6, result7);

  result0 = Xor(result0, result2);
  result4 = Xor(result4, result6);

  result0 = Xor(result0, result4);

  return result0;
}

[jan-wassenberg]

8-bit CLMul, WidenCLMul, CLMulHigh seem useful.
Hoisting detail::GaloisAffine to a supported op sounds fine to me.
Do you know of someone who soon wants to use GaloisAffineInv, GF2P8MulInverse (also Inv for consistency?), GF2P8Mul?
I'm a bit concerned about the Even/Odd/Pairwise add - how would those be implemented efficiently on non-SVE? Is that something we should encourage the use of?