convolution_f2_64.hpp

SECTIONMath INCLUDEnoya/convolution_f2_64.hpp

实现 \(\mathrm{GF}(2^{64})\) 元素的异或加法与无进位乘法；用于该有限域上的卷积或线性代数。

\[ \displaystyle c_k=\sum_{i\oplus j=k}a_i b_j \]

Complexity: Time: O((n + m) log(n + m)). Space: O(n + m).

AC 记录：convolution_F_2_64。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O((n + m) log(n + m)).
/// Space: O(n + m).

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <cstring>
#include <random>
#include <span>
#include <vector>

#if defined(__x86_64__)
#include <immintrin.h>
#elif defined(__aarch64__)
#include <arm_neon.h>
#endif

namespace noya {

namespace convolution_f2_64_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

#if defined(__x86_64__)
__attribute__((target("pclmul"))) inline u128 carryless_multiply(u64 a, u64 b) {
  __m128i prd =
      _mm_clmulepi64_si128(_mm_cvtsi64_si128(a), _mm_cvtsi64_si128(b), 0);
  u128 res;
  std::memcpy(&res, &prd, sizeof(res));
  return res;
}
#elif defined(__aarch64__)
__attribute__((target("+crypto"))) inline u128 carryless_multiply(u64 a,
                                                                  u64 b) {
  poly128_t prd = vmull_p64(poly64_t(a), poly64_t(b));
  u128 res;
  std::memcpy(&res, &prd, sizeof(res));
  return res;
}
#else
inline u128 carryless_multiply(u64 a, u64 b) {
  u128 res = 0;
  for (int bit = 0; bit < 64; bit++) {
    if ((a >> bit) & 1) {
      res ^= u128(b) << bit;
    }
  }
  return res;
}
#endif

} // namespace convolution_f2_64_internal

/// @brief Element of GF(2^64) represented modulo
/// x^64+x^4+x^3+x+1. Addition is XOR; multiplication uses a carry-less
/// product and folds its high half twice through x^64=x^4+x^3+x+1.
class gf2_64 {
public:
  using value_type = std::uint64_t;

  gf2_64() = default;
  explicit gf2_64(value_type val) : x_(val) {}

  value_type value() const { return x_; }

  gf2_64 &operator+=(gf2_64 oth) {
    x_ ^= oth.x_;
    return *this;
  }
  gf2_64 &operator-=(gf2_64 oth) { return *this += oth; }
  gf2_64 &operator*=(gf2_64 oth) {
    using namespace convolution_f2_64_internal;
    u128 prd = carryless_multiply(x_, oth.x_);
    u128 onc = u64(prd) ^ carryless_multiply(u64(prd >> 64), 0x1b);
    x_ = u64(onc) ^ u64(carryless_multiply(u64(onc >> 64), 0x1b));
    return *this;
  }

  friend gf2_64 operator+(gf2_64 a, gf2_64 b) { return a += b; }
  friend gf2_64 operator-(gf2_64 a, gf2_64 b) { return a -= b; }
  friend gf2_64 operator*(gf2_64 a, gf2_64 b) { return a *= b; }
  friend bool operator==(gf2_64, gf2_64) = default;

  gf2_64 power(value_type exp) const {
    gf2_64 res(1);
    gf2_64 bs = *this;
    while (exp != 0) {
      if (exp & 1) {
        res *= bs;
      }
      bs *= bs;
      exp >>= 1;
    }
    return res;
  }

  gf2_64 inverse() const {
    assert(x_ != 0);
    return power(~value_type(0) - 1);
  }

private:
  value_type x_ = 0;
};

namespace convolution_f2_64_internal {

using field = gf2_64;

inline std::vector<field> subset_sums(const std::vector<field> &bas) {
  std::vector<field> res(std::size_t(1) << bas.size());
  for (int bit = 0; bit < int(bas.size()); bit++) {
    for (int mas = 0; mas < (1 << bit); mas++) {
      res[(1 << bit) + mas] = res[mas] + bas[bit];
    }
  }
  return res;
}

struct additive_fft_data {
  std::vector<field> bas;
  std::vector<field> ofs;
  std::vector<field> nba;
  std::vector<field> sof;
  mutable std::vector<field> buf;

  void initialize() {
    int lg = int(bas.size());
    buf.resize(std::size_t(1) << lg);
    field ivl = bas.back().inverse();
    ofs.resize(lg - 1);
    nba.resize(lg - 1);
    for (int idx = 0; idx + 1 < lg; idx++) {
      ofs[idx] = bas[idx] * ivl;
      nba[idx] = ofs[idx] * ofs[idx] + ofs[idx];
    }
    sof = subset_sums(ofs);
  }
};

class additive_fft_cache {
public:
  void prepare(int lg) {
    if (int(d_.size()) > lg) {
      return;
    }
    std::mt19937_64 rng;
    std::vector<field> cha;
    while (int(cha.size()) < lg) {
      cha.clear();
      for (field val(rng()); val != field(); val = val * val + val) {
        cha.push_back(val);
      }
    }
    cha.erase(cha.begin(), cha.end() - lg);
    d_.assign(lg + 1, additive_fft_data{});
    d_[lg].bas = std::move(cha);
    for (int lev = lg; lev > 0; lev--) {
      d_[lev].initialize();
      d_[lev - 1].bas = d_[lev].nba;
    }
  }

  const additive_fft_data &operator[](int lg) const { return d_[lg]; }

private:
  std::vector<additive_fft_data> d_;
};

inline additive_fft_cache fft;

template <bool Inv> void taylor_transform(std::span<field> vs) {
  if constexpr (Inv) {
    for (std::size_t blk = 1; blk * 4 <= vs.size(); blk *= 2) {
      for (std::size_t st = 0; st < vs.size(); st += blk * 4) {
        for (std::size_t idx = 0; idx < blk; idx++) {
          field b = vs[st + blk + idx];
          field thd = vs[st + blk * 2 + idx];
          field fth = vs[st + blk * 3 + idx];
          vs[st + blk + idx] = b + thd;
          vs[st + blk * 2 + idx] = thd + fth;
        }
      }
    }
  } else {
    for (std::size_t blk = vs.size() / 4; blk >= 1; blk /= 2) {
      for (std::size_t st = 0; st < vs.size(); st += blk * 4) {
        for (std::size_t idx = 0; idx < blk; idx++) {
          field b = vs[st + blk + idx];
          field thd = vs[st + blk * 2 + idx];
          field fth = vs[st + blk * 3 + idx];
          vs[st + blk + idx] = b + thd + fth;
          vs[st + blk * 2 + idx] = thd + fth;
        }
      }
    }
  }
}

template <bool Inv = false> void additive_fft(std::span<field> vs) {
  if (vs.size() == 1) {
    return;
  }
  int lg = 63 - __builtin_clzll(vs.size());
  const additive_fft_data &dat = fft[lg];
  if (vs.size() == 2) {
    vs[1] += vs[0];
    return;
  }

  std::size_t hal = vs.size() / 2;
  std::span<field> eve(dat.buf.data(), hal);
  std::span<field> odd(dat.buf.data() + hal, hal);
  if constexpr (!Inv) {
    taylor_transform<false>(vs);
    for (std::size_t idx = 0; idx < hal; idx++) {
      eve[idx] = vs[idx * 2];
      odd[idx] = vs[idx * 2 + 1];
    }
    additive_fft(eve);
    additive_fft(odd);
    for (std::size_t idx = 0; idx < hal; idx++) {
      field a = eve[idx] + dat.sof[idx] * odd[idx];
      vs[idx] = a;
      vs[idx + hal] = a + odd[idx];
    }
  } else {
    for (std::size_t idx = 0; idx < hal; idx++) {
      odd[idx] = vs[idx] + vs[idx + hal];
      eve[idx] = vs[idx] + dat.sof[idx] * odd[idx];
    }
    additive_fft<true>(eve);
    additive_fft<true>(odd);
    for (std::size_t idx = 0; idx < hal; idx++) {
      vs[idx * 2] = eve[idx];
      vs[idx * 2 + 1] = odd[idx];
    }
    taylor_transform<true>(vs);
  }
}

inline std::vector<field> naive_convolution(const std::vector<field> &a,
                                            const std::vector<field> &b) {
  if (a.empty() || b.empty()) {
    return {};
  }
  std::vector<field> res(a.size() + b.size() - 1);
  for (std::size_t i = 0; i < a.size(); i++) {
    for (std::size_t j = 0; j < b.size(); j++) {
      res[i + j] += a[i] * b[j];
    }
  }
  return res;
}

inline std::vector<field> convolve(std::vector<field> a, std::vector<field> b) {
  if (a.empty() || b.empty()) {
    return {};
  }
  std::size_t na = a.size();
  std::size_t nb = b.size();
  int lg = 0;
  while ((std::size_t(1) << lg) < na + nb - 1) {
    lg++;
  }
  std::size_t ntt = std::size_t(1) << lg;
  if (na * nb <= ntt * std::size_t(lg + 1) * (lg + 1)) {
    return naive_convolution(a, b);
  }

  if (lg > 3 && na + nb - 1 == (std::size_t(1) << (lg - 1)) + 1) {
    std::vector<field> tai(nb);
    for (std::size_t idx = 0; idx < nb; idx++) {
      tai[idx] = a.back() * b[idx];
    }
    a.pop_back();
    std::vector<field> res = convolve(std::move(a), std::move(b));
    res.push_back(field());
    for (std::size_t idx = 0; idx < nb; idx++) {
      res[na - 1 + idx] += tai[idx];
    }
    return res;
  }

  fft.prepare(lg);
  a.resize(ntt);
  b.resize(ntt);
  additive_fft(std::span<field>(a));
  additive_fft(std::span<field>(b));
  for (std::size_t idx = 0; idx < ntt; idx++) {
    a[idx] *= b[idx];
  }
  additive_fft<true>(std::span<field>(a));
  a.resize(na + nb - 1);
  return a;
}

} // namespace convolution_f2_64_internal

/// @brief Convolution over GF(2^64). The transform evaluates in a tower basis
/// built from the Artin-Schreier map x -> x^2+x; its subspace polynomials split
/// each evaluation set into two affine halves, giving radix-two butterflies in
/// characteristic two without requiring roots of unity.
inline std::vector<std::uint64_t>
convolution_f2_64(const std::vector<std::uint64_t> &a,
                  const std::vector<std::uint64_t> &b) {
  using namespace convolution_f2_64_internal;
  std::vector<field> fa;
  std::vector<field> fb;
  fa.reserve(a.size());
  fb.reserve(b.size());
  for (std::uint64_t val : a) {
    fa.emplace_back(val);
  }
  for (std::uint64_t val : b) {
    fb.emplace_back(val);
  }
  std::vector<field> fc = convolve(std::move(fa), std::move(fb));
  std::vector<std::uint64_t> res(fc.size());
  for (std::size_t idx = 0; idx < res.size(); idx++) {
    res[idx] = fc[idx].value();
  }
  return res;
}

} // namespace noya

#ifndef NOYA_CONVOLUTION_F2_64_HPP
#define NOYA_CONVOLUTION_F2_64_HPP 1

/// @complexity Time: O((n + m) log(n + m)).
/// Space: O(n + m).

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <cstring>
#include <random>
#include <span>
#include <vector>

#if defined(__x86_64__)
#include <immintrin.h>
#elif defined(__aarch64__)
#include <arm_neon.h>
#endif

namespace noya {

namespace convolution_f2_64_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

#if defined(__x86_64__)
__attribute__((target("pclmul"))) inline u128 carryless_multiply(u64 a, u64 b) {
  __m128i prd =
      _mm_clmulepi64_si128(_mm_cvtsi64_si128(a), _mm_cvtsi64_si128(b), 0);
  u128 res;
  std::memcpy(&res, &prd, sizeof(res));
  return res;
}
#elif defined(__aarch64__)
__attribute__((target("+crypto"))) inline u128 carryless_multiply(u64 a,
                                                                  u64 b) {
  poly128_t prd = vmull_p64(poly64_t(a), poly64_t(b));
  u128 res;
  std::memcpy(&res, &prd, sizeof(res));
  return res;
}
#else
inline u128 carryless_multiply(u64 a, u64 b) {
  u128 res = 0;
  for (int bit = 0; bit < 64; bit++) {
    if ((a >> bit) & 1) {
      res ^= u128(b) << bit;
    }
  }
  return res;
}
#endif

} // namespace convolution_f2_64_internal

/// @brief Element of GF(2^64) represented modulo
/// x^64+x^4+x^3+x+1. Addition is XOR; multiplication uses a carry-less
/// product and folds its high half twice through x^64=x^4+x^3+x+1.
class gf2_64 {
public:
  using value_type = std::uint64_t;

  gf2_64() = default;
  explicit gf2_64(value_type val) : x_(val) {}

  value_type value() const { return x_; }

  gf2_64 &operator+=(gf2_64 oth) {
    x_ ^= oth.x_;
    return *this;
  }
  gf2_64 &operator-=(gf2_64 oth) { return *this += oth; }
  gf2_64 &operator*=(gf2_64 oth) {
    using namespace convolution_f2_64_internal;
    u128 prd = carryless_multiply(x_, oth.x_);
    u128 onc = u64(prd) ^ carryless_multiply(u64(prd >> 64), 0x1b);
    x_ = u64(onc) ^ u64(carryless_multiply(u64(onc >> 64), 0x1b));
    return *this;
  }

  friend gf2_64 operator+(gf2_64 a, gf2_64 b) { return a += b; }
  friend gf2_64 operator-(gf2_64 a, gf2_64 b) { return a -= b; }
  friend gf2_64 operator*(gf2_64 a, gf2_64 b) { return a *= b; }
  friend bool operator==(gf2_64, gf2_64) = default;

  gf2_64 power(value_type exp) const {
    gf2_64 res(1);
    gf2_64 bs = *this;
    while (exp != 0) {
      if (exp & 1) {
        res *= bs;
      }
      bs *= bs;
      exp >>= 1;
    }
    return res;
  }

  gf2_64 inverse() const {
    assert(x_ != 0);
    return power(~value_type(0) - 1);
  }

private:
  value_type x_ = 0;
};

namespace convolution_f2_64_internal {

using field = gf2_64;

inline std::vector<field> subset_sums(const std::vector<field> &bas) {
  std::vector<field> res(std::size_t(1) << bas.size());
  for (int bit = 0; bit < int(bas.size()); bit++) {
    for (int mas = 0; mas < (1 << bit); mas++) {
      res[(1 << bit) + mas] = res[mas] + bas[bit];
    }
  }
  return res;
}

struct additive_fft_data {
  std::vector<field> bas;
  std::vector<field> ofs;
  std::vector<field> nba;
  std::vector<field> sof;
  mutable std::vector<field> buf;

  void initialize() {
    int lg = int(bas.size());
    buf.resize(std::size_t(1) << lg);
    field ivl = bas.back().inverse();
    ofs.resize(lg - 1);
    nba.resize(lg - 1);
    for (int idx = 0; idx + 1 < lg; idx++) {
      ofs[idx] = bas[idx] * ivl;
      nba[idx] = ofs[idx] * ofs[idx] + ofs[idx];
    }
    sof = subset_sums(ofs);
  }
};

class additive_fft_cache {
public:
  void prepare(int lg) {
    if (int(d_.size()) > lg) {
      return;
    }
    std::mt19937_64 rng;
    std::vector<field> cha;
    while (int(cha.size()) < lg) {
      cha.clear();
      for (field val(rng()); val != field(); val = val * val + val) {
        cha.push_back(val);
      }
    }
    cha.erase(cha.begin(), cha.end() - lg);
    d_.assign(lg + 1, additive_fft_data{});
    d_[lg].bas = std::move(cha);
    for (int lev = lg; lev > 0; lev--) {
      d_[lev].initialize();
      d_[lev - 1].bas = d_[lev].nba;
    }
  }

  const additive_fft_data &operator[](int lg) const { return d_[lg]; }

private:
  std::vector<additive_fft_data> d_;
};

inline additive_fft_cache fft;

template <bool Inv> void taylor_transform(std::span<field> vs) {
  if constexpr (Inv) {
    for (std::size_t blk = 1; blk * 4 <= vs.size(); blk *= 2) {
      for (std::size_t st = 0; st < vs.size(); st += blk * 4) {
        for (std::size_t idx = 0; idx < blk; idx++) {
          field b = vs[st + blk + idx];
          field thd = vs[st + blk * 2 + idx];
          field fth = vs[st + blk * 3 + idx];
          vs[st + blk + idx] = b + thd;
          vs[st + blk * 2 + idx] = thd + fth;
        }
      }
    }
  } else {
    for (std::size_t blk = vs.size() / 4; blk >= 1; blk /= 2) {
      for (std::size_t st = 0; st < vs.size(); st += blk * 4) {
        for (std::size_t idx = 0; idx < blk; idx++) {
          field b = vs[st + blk + idx];
          field thd = vs[st + blk * 2 + idx];
          field fth = vs[st + blk * 3 + idx];
          vs[st + blk + idx] = b + thd + fth;
          vs[st + blk * 2 + idx] = thd + fth;
        }
      }
    }
  }
}

template <bool Inv = false> void additive_fft(std::span<field> vs) {
  if (vs.size() == 1) {
    return;
  }
  int lg = 63 - __builtin_clzll(vs.size());
  const additive_fft_data &dat = fft[lg];
  if (vs.size() == 2) {
    vs[1] += vs[0];
    return;
  }

  std::size_t hal = vs.size() / 2;
  std::span<field> eve(dat.buf.data(), hal);
  std::span<field> odd(dat.buf.data() + hal, hal);
  if constexpr (!Inv) {
    taylor_transform<false>(vs);
    for (std::size_t idx = 0; idx < hal; idx++) {
      eve[idx] = vs[idx * 2];
      odd[idx] = vs[idx * 2 + 1];
    }
    additive_fft(eve);
    additive_fft(odd);
    for (std::size_t idx = 0; idx < hal; idx++) {
      field a = eve[idx] + dat.sof[idx] * odd[idx];
      vs[idx] = a;
      vs[idx + hal] = a + odd[idx];
    }
  } else {
    for (std::size_t idx = 0; idx < hal; idx++) {
      odd[idx] = vs[idx] + vs[idx + hal];
      eve[idx] = vs[idx] + dat.sof[idx] * odd[idx];
    }
    additive_fft<true>(eve);
    additive_fft<true>(odd);
    for (std::size_t idx = 0; idx < hal; idx++) {
      vs[idx * 2] = eve[idx];
      vs[idx * 2 + 1] = odd[idx];
    }
    taylor_transform<true>(vs);
  }
}

inline std::vector<field> naive_convolution(const std::vector<field> &a,
                                            const std::vector<field> &b) {
  if (a.empty() || b.empty()) {
    return {};
  }
  std::vector<field> res(a.size() + b.size() - 1);
  for (std::size_t i = 0; i < a.size(); i++) {
    for (std::size_t j = 0; j < b.size(); j++) {
      res[i + j] += a[i] * b[j];
    }
  }
  return res;
}

inline std::vector<field> convolve(std::vector<field> a, std::vector<field> b) {
  if (a.empty() || b.empty()) {
    return {};
  }
  std::size_t na = a.size();
  std::size_t nb = b.size();
  int lg = 0;
  while ((std::size_t(1) << lg) < na + nb - 1) {
    lg++;
  }
  std::size_t ntt = std::size_t(1) << lg;
  if (na * nb <= ntt * std::size_t(lg + 1) * (lg + 1)) {
    return naive_convolution(a, b);
  }

  if (lg > 3 && na + nb - 1 == (std::size_t(1) << (lg - 1)) + 1) {
    std::vector<field> tai(nb);
    for (std::size_t idx = 0; idx < nb; idx++) {
      tai[idx] = a.back() * b[idx];
    }
    a.pop_back();
    std::vector<field> res = convolve(std::move(a), std::move(b));
    res.push_back(field());
    for (std::size_t idx = 0; idx < nb; idx++) {
      res[na - 1 + idx] += tai[idx];
    }
    return res;
  }

  fft.prepare(lg);
  a.resize(ntt);
  b.resize(ntt);
  additive_fft(std::span<field>(a));
  additive_fft(std::span<field>(b));
  for (std::size_t idx = 0; idx < ntt; idx++) {
    a[idx] *= b[idx];
  }
  additive_fft<true>(std::span<field>(a));
  a.resize(na + nb - 1);
  return a;
}

} // namespace convolution_f2_64_internal

/// @brief Convolution over GF(2^64). The transform evaluates in a tower basis
/// built from the Artin-Schreier map x -> x^2+x; its subspace polynomials split
/// each evaluation set into two affine halves, giving radix-two butterflies in
/// characteristic two without requiring roots of unity.
inline std::vector<std::uint64_t>
convolution_f2_64(const std::vector<std::uint64_t> &a,
                  const std::vector<std::uint64_t> &b) {
  using namespace convolution_f2_64_internal;
  std::vector<field> fa;
  std::vector<field> fb;
  fa.reserve(a.size());
  fb.reserve(b.size());
  for (std::uint64_t val : a) {
    fa.emplace_back(val);
  }
  for (std::uint64_t val : b) {
    fb.emplace_back(val);
  }
  std::vector<field> fc = convolve(std::move(fa), std::move(fb));
  std::vector<std::uint64_t> res(fc.size());
  for (std::size_t idx = 0; idx < res.size(); idx++) {
    res[idx] = fc[idx].value();
  }
  return res;
}

} // namespace noya

#endif // NOYA_CONVOLUTION_F2_64_HPP

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <cstring>
#include <random>
#include <span>
#include <vector>

/// @complexity Time: O((n + m) log(n + m)).
/// Space: O(n + m).

#if defined(__x86_64__)
#include <immintrin.h>
#elif defined(__aarch64__)
#include <arm_neon.h>
#endif

namespace noya {

namespace convolution_f2_64_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

#if defined(__x86_64__)
__attribute__((target("pclmul"))) inline u128 carryless_multiply(u64 a, u64 b) {
  __m128i prd =
      _mm_clmulepi64_si128(_mm_cvtsi64_si128(a), _mm_cvtsi64_si128(b), 0);
  u128 res;
  std::memcpy(&res, &prd, sizeof(res));
  return res;
}
#elif defined(__aarch64__)
__attribute__((target("+crypto"))) inline u128 carryless_multiply(u64 a,
                                                                  u64 b) {
  poly128_t prd = vmull_p64(poly64_t(a), poly64_t(b));
  u128 res;
  std::memcpy(&res, &prd, sizeof(res));
  return res;
}
#else
inline u128 carryless_multiply(u64 a, u64 b) {
  u128 res = 0;
  for (int bit = 0; bit < 64; bit++) {
    if ((a >> bit) & 1) {
      res ^= u128(b) << bit;
    }
  }
  return res;
}
#endif

} // namespace convolution_f2_64_internal

/// @brief Element of GF(2^64) represented modulo
/// x^64+x^4+x^3+x+1. Addition is XOR; multiplication uses a carry-less
/// product and folds its high half twice through x^64=x^4+x^3+x+1.
class gf2_64 {
public:
  using value_type = std::uint64_t;

  gf2_64() = default;
  explicit gf2_64(value_type val) : x_(val) {}

  value_type value() const { return x_; }

  gf2_64 &operator+=(gf2_64 oth) {
    x_ ^= oth.x_;
    return *this;
  }
  gf2_64 &operator-=(gf2_64 oth) { return *this += oth; }
  gf2_64 &operator*=(gf2_64 oth) {
    using namespace convolution_f2_64_internal;
    u128 prd = carryless_multiply(x_, oth.x_);
    u128 onc = u64(prd) ^ carryless_multiply(u64(prd >> 64), 0x1b);
    x_ = u64(onc) ^ u64(carryless_multiply(u64(onc >> 64), 0x1b));
    return *this;
  }

  friend gf2_64 operator+(gf2_64 a, gf2_64 b) { return a += b; }
  friend gf2_64 operator-(gf2_64 a, gf2_64 b) { return a -= b; }
  friend gf2_64 operator*(gf2_64 a, gf2_64 b) { return a *= b; }
  friend bool operator==(gf2_64, gf2_64) = default;

  gf2_64 power(value_type exp) const {
    gf2_64 res(1);
    gf2_64 bs = *this;
    while (exp != 0) {
      if (exp & 1) {
        res *= bs;
      }
      bs *= bs;
      exp >>= 1;
    }
    return res;
  }

  gf2_64 inverse() const {
    assert(x_ != 0);
    return power(~value_type(0) - 1);
  }

private:
  value_type x_ = 0;
};

namespace convolution_f2_64_internal {

using field = gf2_64;

inline std::vector<field> subset_sums(const std::vector<field> &bas) {
  std::vector<field> res(std::size_t(1) << bas.size());
  for (int bit = 0; bit < int(bas.size()); bit++) {
    for (int mas = 0; mas < (1 << bit); mas++) {
      res[(1 << bit) + mas] = res[mas] + bas[bit];
    }
  }
  return res;
}

struct additive_fft_data {
  std::vector<field> bas;
  std::vector<field> ofs;
  std::vector<field> nba;
  std::vector<field> sof;
  mutable std::vector<field> buf;

  void initialize() {
    int lg = int(bas.size());
    buf.resize(std::size_t(1) << lg);
    field ivl = bas.back().inverse();
    ofs.resize(lg - 1);
    nba.resize(lg - 1);
    for (int idx = 0; idx + 1 < lg; idx++) {
      ofs[idx] = bas[idx] * ivl;
      nba[idx] = ofs[idx] * ofs[idx] + ofs[idx];
    }
    sof = subset_sums(ofs);
  }
};

class additive_fft_cache {
public:
  void prepare(int lg) {
    if (int(d_.size()) > lg) {
      return;
    }
    std::mt19937_64 rng;
    std::vector<field> cha;
    while (int(cha.size()) < lg) {
      cha.clear();
      for (field val(rng()); val != field(); val = val * val + val) {
        cha.push_back(val);
      }
    }
    cha.erase(cha.begin(), cha.end() - lg);
    d_.assign(lg + 1, additive_fft_data{});
    d_[lg].bas = std::move(cha);
    for (int lev = lg; lev > 0; lev--) {
      d_[lev].initialize();
      d_[lev - 1].bas = d_[lev].nba;
    }
  }

  const additive_fft_data &operator[](int lg) const { return d_[lg]; }

private:
  std::vector<additive_fft_data> d_;
};

inline additive_fft_cache fft;

template <bool Inv> void taylor_transform(std::span<field> vs) {
  if constexpr (Inv) {
    for (std::size_t blk = 1; blk * 4 <= vs.size(); blk *= 2) {
      for (std::size_t st = 0; st < vs.size(); st += blk * 4) {
        for (std::size_t idx = 0; idx < blk; idx++) {
          field b = vs[st + blk + idx];
          field thd = vs[st + blk * 2 + idx];
          field fth = vs[st + blk * 3 + idx];
          vs[st + blk + idx] = b + thd;
          vs[st + blk * 2 + idx] = thd + fth;
        }
      }
    }
  } else {
    for (std::size_t blk = vs.size() / 4; blk >= 1; blk /= 2) {
      for (std::size_t st = 0; st < vs.size(); st += blk * 4) {
        for (std::size_t idx = 0; idx < blk; idx++) {
          field b = vs[st + blk + idx];
          field thd = vs[st + blk * 2 + idx];
          field fth = vs[st + blk * 3 + idx];
          vs[st + blk + idx] = b + thd + fth;
          vs[st + blk * 2 + idx] = thd + fth;
        }
      }
    }
  }
}

template <bool Inv = false> void additive_fft(std::span<field> vs) {
  if (vs.size() == 1) {
    return;
  }
  int lg = 63 - __builtin_clzll(vs.size());
  const additive_fft_data &dat = fft[lg];
  if (vs.size() == 2) {
    vs[1] += vs[0];
    return;
  }

  std::size_t hal = vs.size() / 2;
  std::span<field> eve(dat.buf.data(), hal);
  std::span<field> odd(dat.buf.data() + hal, hal);
  if constexpr (!Inv) {
    taylor_transform<false>(vs);
    for (std::size_t idx = 0; idx < hal; idx++) {
      eve[idx] = vs[idx * 2];
      odd[idx] = vs[idx * 2 + 1];
    }
    additive_fft(eve);
    additive_fft(odd);
    for (std::size_t idx = 0; idx < hal; idx++) {
      field a = eve[idx] + dat.sof[idx] * odd[idx];
      vs[idx] = a;
      vs[idx + hal] = a + odd[idx];
    }
  } else {
    for (std::size_t idx = 0; idx < hal; idx++) {
      odd[idx] = vs[idx] + vs[idx + hal];
      eve[idx] = vs[idx] + dat.sof[idx] * odd[idx];
    }
    additive_fft<true>(eve);
    additive_fft<true>(odd);
    for (std::size_t idx = 0; idx < hal; idx++) {
      vs[idx * 2] = eve[idx];
      vs[idx * 2 + 1] = odd[idx];
    }
    taylor_transform<true>(vs);
  }
}

inline std::vector<field> naive_convolution(const std::vector<field> &a,
                                            const std::vector<field> &b) {
  if (a.empty() || b.empty()) {
    return {};
  }
  std::vector<field> res(a.size() + b.size() - 1);
  for (std::size_t i = 0; i < a.size(); i++) {
    for (std::size_t j = 0; j < b.size(); j++) {
      res[i + j] += a[i] * b[j];
    }
  }
  return res;
}

inline std::vector<field> convolve(std::vector<field> a, std::vector<field> b) {
  if (a.empty() || b.empty()) {
    return {};
  }
  std::size_t na = a.size();
  std::size_t nb = b.size();
  int lg = 0;
  while ((std::size_t(1) << lg) < na + nb - 1) {
    lg++;
  }
  std::size_t ntt = std::size_t(1) << lg;
  if (na * nb <= ntt * std::size_t(lg + 1) * (lg + 1)) {
    return naive_convolution(a, b);
  }

  if (lg > 3 && na + nb - 1 == (std::size_t(1) << (lg - 1)) + 1) {
    std::vector<field> tai(nb);
    for (std::size_t idx = 0; idx < nb; idx++) {
      tai[idx] = a.back() * b[idx];
    }
    a.pop_back();
    std::vector<field> res = convolve(std::move(a), std::move(b));
    res.push_back(field());
    for (std::size_t idx = 0; idx < nb; idx++) {
      res[na - 1 + idx] += tai[idx];
    }
    return res;
  }

  fft.prepare(lg);
  a.resize(ntt);
  b.resize(ntt);
  additive_fft(std::span<field>(a));
  additive_fft(std::span<field>(b));
  for (std::size_t idx = 0; idx < ntt; idx++) {
    a[idx] *= b[idx];
  }
  additive_fft<true>(std::span<field>(a));
  a.resize(na + nb - 1);
  return a;
}

} // namespace convolution_f2_64_internal

/// @brief Convolution over GF(2^64). The transform evaluates in a tower basis
/// built from the Artin-Schreier map x -> x^2+x; its subspace polynomials split
/// each evaluation set into two affine halves, giving radix-two butterflies in
/// characteristic two without requiring roots of unity.
inline std::vector<std::uint64_t>
convolution_f2_64(const std::vector<std::uint64_t> &a,
                  const std::vector<std::uint64_t> &b) {
  using namespace convolution_f2_64_internal;
  std::vector<field> fa;
  std::vector<field> fb;
  fa.reserve(a.size());
  fb.reserve(b.size());
  for (std::uint64_t val : a) {
    fa.emplace_back(val);
  }
  for (std::uint64_t val : b) {
    fb.emplace_back(val);
  }
  std::vector<field> fc = convolve(std::move(fa), std::move(fb));
  std::vector<std::uint64_t> res(fc.size());
  for (std::size_t idx = 0; idx < res.size(); idx++) {
    res[idx] = fc[idx].value();
  }
  return res;
}

} // namespace noya