多项式

hj23308 · 2024-05-05 15:09:14 · 算法·理论

FFT A

[c-array] (FFT) Discrete Fourier transform

namespace FFT {
  template<typename T>
  void DFT(std::complex<T> y[], int len) {
    const T pi = std::acos(T(-1));
    std::vector<int> rev(len);
    for (int i = 1; i < len; i++) 
      rev[i] = (rev[i >> 1] + len * (i & 1)) >> 1;
    for (int i = 1; i < len; i++) 
      if (rev[i] < i) 
        std::swap(y[i], y[rev[i]]);
    for (int h = 1; h < len; h <<= 1) {
      std::complex<T> step(std::cos(pi / h), std::sin(pi / h)), w, u, v;
      for (int i = 0; w = {1, 0}, i < len; i += h << 1)
        for (int j = i; j < i + h; j++, w *= step)
          u = y[j], v = w * y[j + h], y[j] = u + v, y[j + h] = u - v;
    }
  }

  template<typename T>
  void IDFT(std::complex<T> y[], int len) {
    DFT(y, len), std::reverse(y + 1, y + len);
    for (int i = 0; i < len; i++) y[i] /= len;
  }
}

FFT V

[std::vector] (FFT) Discrete Fourier transform

namespace FFT {
  template<typename T>
  void DFT(std::vector<std::complex<T>> &y, int len) {
    const T pi = std::acos(T(-1));
    std::vector<int> rev(len);
    for (int i = 1; i < len; i++) 
      rev[i] = (rev[i >> 1] + len * (i & 1)) >> 1;
    for (int i = 1; i < len; i++) 
      if (rev[i] < i) 
        std::swap(y[i], y[rev[i]]);
    for (int h = 1; h < len; h <<= 1) {
      std::complex<T> step(std::cos(pi / h), std::sin(pi / h)), w, u, v;
      for (int i = 0; w = {1, 0}, i < len; i += h << 1)
        for (int j = i; j < i + h; j++, w *= step)
          u = y[j], v = w * y[j + h], y[j] = u + v, y[j + h] = u - v;
    }
  }

  template<typename T>
  void IDFT(std::vector<std::complex<T>> &y, int len) {
    DFT(y, len), std::reverse(y.begin() + 1, y.begin() + len);
    for (int i = 0; i < len; i++) y[i] /= len;
  }
}

FFT I

[iterator] (FFT) Discrete Fourier transform

namespace FFT {
  template<typename Iterator>
  void DFT(Iterator begin, int len) {
    using T = typename std::iterator_traits<Iterator>::value_type;
    using D = decltype(begin->real());
    const D pi = acos(D(-1));
    std::vector<int> rev(len);

    for (int i = 1; i < len; i++)
      rev[i] = (rev[i >> 1] + len * (i & 1)) >> 1;
    for (int i = 1; i < len; i++)
      if (rev[i] < i)
        std::swap(*(begin + i), *(begin + rev[i]));

    for (int h = 1; h < len; h <<= 1) {
      T step(std::cos(pi / h), std::sin(pi / h)), w, u, v;
      for (int i = 0; w = {1, 0}, i < len; i += h << 1)
        for (int j = i; j < i + h; j++, w *= step)
          u = *(begin + j), v = w * *(begin + j + h),
          *(begin + j) = u + v, *(begin + j + h) = u - v;
    }
  }

  template<typename Iterator>
  void IDFT(Iterator begin, int len) {
    DFT(begin, len), std::reverse(begin + 1, begin + len);
    for (int i = 0; i < len; i++) *(begin + i) /= len;
  }

  template<typename Iterator>
  void DFT(Iterator begin, Iterator end) {
    DFT(begin, std::distance(end - begin));
  }

  template<typename Iterator>
  void IDFT(Iterator begin, Iterator end) {
    IDFT(begin, std::distance(end - begin));
  }
}

NTT A

[c-array] (NTT) Number Theoretic transform

namespace NTT {
  void DFT(int y[], int len) {
    std::vector<int> rev(len);
    for (int i = 1; i < len; i++) 
      rev[i] = (rev[i >> 1] + len * (i & 1)) >> 1;
    for (int i = 1; i < len; i++) 
      if (rev[i] < i) 
        std::swap(y[i], y[rev[i]]);
    for (int h = 1; h < len; h <<= 1) {
      int step = Power(3, (pMod - 1) / h / 2), u, v;
      for (int i = 0; i < len; i += h << 1)
        for (int j = i, w = 1; j < i + h; j++, MulC(w, step)) {
          u = y[j], v = Mul(y[j + h], w);
          y[j] = Add(u, v), y[j + h] = Add(u, -v);
        }
    }
  }

  void IDFT(int y[], int len) {
    DFT(y, len), std::reverse(y + 1, y + len);
    int inv = Power(len, pMod - 2);
    for (int i = 0; i < len; i++) MulC(y[i], inv);
  }
}

NNT V

[std::vector] (NTT) Number Theoretic transform

namespace NTT {
  void DFT(std::vector<int> &y, int len) {
    std::vector<int> rev(len);
    for (int i = 1; i < len; i++)
      rev[i] = (rev[i >> 1] + len * (i & 1)) >> 1;
    for (int i = 1; i < len; i++)
      if (rev[i] < i)
        std::swap(y[i], y[rev[i]]);
    for (int h = 1; h < len; h <<= 1) {
      int step = Power(3, (pMod - 1) / h / 2), u, v;
      for (int i = 0; i < len; i += h << 1)
        for (int j = i, w = 1; j < i + h; j++, MulC(w, step)) {
          u = y[j], v = Mul(y[j + h], w);
          y[j] = Add(u, v), y[j + h] = Add(u, pMod - v);
        }
    }
  }

  void IDFT(std::vector<int> &y, int len) {
    DFT(y, len), std::reverse(y.begin() + 1, y.begin() + len);
    int inv = Power(len, pMod - 2);
    for (int i = 0; i < len; i++) MulC(y[i], inv);
  }
}

LagrangeInterpolation

[n^2] lagrange interpolation

int LagrangeInterpolation(const std::vector<int> &x,
                          const std::vector<int> &y,
                          int k) {
  if (x.size() != y.size())
    throw std::runtime_error("x.size() != y.size()");
  int n = (int) x.size(), res = 0;
  for (int i = 0; i < n; i++) {
    int mul = y[i], inv = 1;
    for (int j = 0; j < n; j++)
      if (i != j) {
        MulC(mul, -x[j] + k);
        MulC(inv, -x[j] + x[i]);
      }
    AddC(res, Mul(mul, Inv(inv)));
  }
  return res;
}

LagrangeInterpolation-Continuous

[0-index] continuous lagrange interpolation

int LagrangeInterpolation(const std::vector<int> &y, int k) {
  int n = (int) y.size(), res = 0;
  if (k < n) return y[k];
  std::vector<int> pre(n), suf(n);
  for (int i = 0; i < n; i++) pre[i] = suf[i] = k - i;
  for (int i = 1; i < n; i++) MulC(pre[i], pre[i - 1]);
  for (int i = n; i > 1; i--) MulC(suf[i - 2], suf[i - 1]);
  for (int i = 0; i < n; i++) {
    int mul = (n - i) & 1 ? y[i] : -y[i];
    int inv = Mul(invFac[i], invFac[n - 1 - i]);
    if (i > 0) MulC(mul, pre[i - 1]);
    if (i < n - 1) MulC(mul, suf[i + 1]);
    AddC(res, Mul(mul, inv));
  }
  return res;
}

LagrangeInterpolation-Poly

[n^2] (return polynomial) lagrangian interpolation

std::vector<int> LagrangeInterpolation(const std::vector<int> &x,
                                       const std::vector<int> &y) {
  int n = (int) y.size();
  std::vector<int> ret(n), f(n + 1);
  f[0] = 1;
  for (int i = 0; i < n; i++) {
    for (int j = i + 1; j > 0; j--)
      f[j] = Add(f[j - 1], -Mul(f[j], x[i]));
    MulC(f[0], -x[i]);
  }
  for (int i = 0; i < n; i++) {
    int mul = 1;
    for (int j = 0; j < n; j++)
      if (i != j)
        MulC(mul, x[i] - x[j]);
    mul = Mul(y[i], Inv(mul));
    auto g = f;
    for (int j = n; j > 0; j--) {
      if (x[i])
        AddC(g[j - 1], Mul(g[j], x[i]));
      AddC(ret[j - 1], Mul(g[j], mul));
    }
  }
  return ret;
}

std::vector<int> LagrangeInterpolation(const std::vector<int> &y) {
  int n = (int) y.size();
  std::vector<int> x(n);
  for (int i = 0; i < n; i++) x[i] = i;
  return LagrangeInterpolation(x, y);
}