Fast Walsh-Hadamard Transform
(convolution/fwht.hpp)

• View this file on GitHub
• Last update: 2023-07-25 00:50:50+07:00
• Include: #include "convolution/fwht.hpp"

Code

template <typename T, typename F>
void abstract_fwht(vector<T> &seq, F f) {
  const int n = seq.size();
  assert(__builtin_popcount(n) == 1);
  for (int w = 1; w < n; w *= 2) {
    for (int i = 0; i < n; i += w * 2) {
      for (int j = 0; j < w; j++) {
        f(seq[i + j], seq[i + j + w]);
      }
    }
  }
}

template <typename T, typename F1, typename F2>
vector<T> bitwise_conv(vector<T> x, vector<T> y, F1 f, F2 finv) {
  const int n = x.size();
  assert(__builtin_popcount(n) == 1);
  assert(x.size() == y.size());
  if (x == y) {
    abstract_fwht(x, f), y = x;
  } else {
    abstract_fwht(x, f), abstract_fwht(y, f);
  }
  for (size_t i = 0; i < x.size(); i++) {
    x[i] *= y[i];
  }
  abstract_fwht(x, finv);
  return x;
}

// bitwise xor convolution (FWHT-based)
// ret[i] = \sum_j x[j] * y[i ^ j]
// if T is integer, ||x||_1 * ||y||_1 * 2 < numeric_limits<T>::max()
template <typename T>
vector<T> xorconv(vector<T> x, vector<T> y) {
  auto f = [](T &lo, T &hi) {
    T c = lo + hi;
    hi = lo - hi, lo = c;
  };
  auto finv = [](T &lo, T &hi) {
    T c = lo + hi;
    hi = (lo - hi) / 2,
    lo =
        c / 2;  // Reconsider HEAVY complexity of division by 2 when T is ModInt
  };
  return bitwise_conv(x, y, f, finv);
}

// bitwise AND conolution
// ret[i] = \sum_{(j & k) == i} x[j] * y[k]
template <typename T>
vector<T> andconv(vector<T> x, vector<T> y) {
  return bitwise_conv(
      x, y, [](T &lo, T &hi) { lo += hi; }, [](T &lo, T &hi) { lo -= hi; });
}

// bitwise OR convolution
// ret[i] = \sum_{(j | k) == i} x[j] * y[k]
template <typename T>
vector<T> orconv(vector<T> x, vector<T> y) {
  return bitwise_conv(
      x, y, [](T &lo, T &hi) { hi += lo; }, [](T &lo, T &hi) { hi -= lo; });
}

#line 1 "convolution/fwht.hpp"
template <typename T, typename F>
void abstract_fwht(vector<T> &seq, F f) {
  const int n = seq.size();
  assert(__builtin_popcount(n) == 1);
  for (int w = 1; w < n; w *= 2) {
    for (int i = 0; i < n; i += w * 2) {
      for (int j = 0; j < w; j++) {
        f(seq[i + j], seq[i + j + w]);
      }
    }
  }
}

template <typename T, typename F1, typename F2>
vector<T> bitwise_conv(vector<T> x, vector<T> y, F1 f, F2 finv) {
  const int n = x.size();
  assert(__builtin_popcount(n) == 1);
  assert(x.size() == y.size());
  if (x == y) {
    abstract_fwht(x, f), y = x;
  } else {
    abstract_fwht(x, f), abstract_fwht(y, f);
  }
  for (size_t i = 0; i < x.size(); i++) {
    x[i] *= y[i];
  }
  abstract_fwht(x, finv);
  return x;
}

// bitwise xor convolution (FWHT-based)
// ret[i] = \sum_j x[j] * y[i ^ j]
// if T is integer, ||x||_1 * ||y||_1 * 2 < numeric_limits<T>::max()
template <typename T>
vector<T> xorconv(vector<T> x, vector<T> y) {
  auto f = [](T &lo, T &hi) {
    T c = lo + hi;
    hi = lo - hi, lo = c;
  };
  auto finv = [](T &lo, T &hi) {
    T c = lo + hi;
    hi = (lo - hi) / 2,
    lo =
        c / 2;  // Reconsider HEAVY complexity of division by 2 when T is ModInt
  };
  return bitwise_conv(x, y, f, finv);
}

// bitwise AND conolution
// ret[i] = \sum_{(j & k) == i} x[j] * y[k]
template <typename T>
vector<T> andconv(vector<T> x, vector<T> y) {
  return bitwise_conv(
      x, y, [](T &lo, T &hi) { lo += hi; }, [](T &lo, T &hi) { lo -= hi; });
}

// bitwise OR convolution
// ret[i] = \sum_{(j | k) == i} x[j] * y[k]
template <typename T>
vector<T> orconv(vector<T> x, vector<T> y) {
  return bitwise_conv(
      x, y, [](T &lo, T &hi) { hi += lo; }, [](T &lo, T &hi) { hi -= lo; });
}

Back to top page