ntt/ntt.hpp

• View this file on GitHub
• Last update: 2026-06-06 19:38:56+09:00
• Include: #include "ntt/ntt.hpp"

Required by

• fps/arbitrary-fps.hpp
• Multipoint Evaluation(高速化版) (fps/fast-multieval.hpp)
• P-recursiveの高速計算 (fps/find-p-recursive.hpp)
• 線形漸化式の高速計算 (fps/kitamasa.hpp)
• Mod-Pow ($f(x)^k \mod g(x)$) (fps/mod-pow.hpp)
• 線形回帰数列の高速計算(Berlekamp-Massey/Bostan-Mori) (fps/nth-term.hpp)
• NTT mod用FPSライブラリ (fps/ntt-friendly-fps.hpp)
• fps/online-fps.hpp
• fps/root-finding.hpp
• fps/sample-point-shift.hpp
• math/bigint-all.hpp
• math/bigint-binary.hpp
• Garner's algorithm for bigint (math/bigint-garner.hpp)
• math/bigint-gcd.hpp
• math/bigint-rational.hpp
• math/bigint-to-hex.hpp
• 多倍長整数 (math/bigint.hpp)
• Black Box Linear Algebra (matrix/black-box-linear-algebra.hpp)
• 多項式行列のprefix product (matrix/polynomial-matrix-prefix-prod.hpp)
• 階乗 $\mod p$ (modulo/factorial.hpp)
• ntt/arbitrary-ntt.hpp
• Cooley-Tukey FFT Algorithm (ntt/cooley-tukey-ntt.hpp)
• 多変数巡回畳み込み (ntt/multivariate-circular-convolution.hpp)
• Rader's FFT Algorithm (ntt/rader-ntt.hpp)
• Wildcard Pattern Matching (string/wildcard-pattern-matching.hpp)

Verified with

• verify/verify-aoj-ntl/aoj-ntl-2-a.test.cpp
• verify/verify-aoj-ntl/aoj-ntl-2-b.test.cpp
• verify/verify-aoj-ntl/aoj-ntl-2-c.test.cpp
• verify/verify-aoj-ntl/aoj-ntl-2-d.test.cpp
• verify/verify-aoj-ntl/aoj-ntl-2-e.test.cpp
• verify/verify-aoj-ntl/aoj-ntl-2-f.test.cpp
• verify/verify-aoj-other/aoj-2171-bigrational.test.cpp
• verify/verify-unit-test/arbitrary-ntt-mod18446744069414584321.test.cpp
• verify/verify-unit-test/bigint-gcd.test.cpp
• verify/verify-unit-test/bigint.test.cpp
• verify/verify-unit-test/bigint2.test.cpp
• verify/verify-unit-test/bigint3.test.cpp
• verify/verify-unit-test/bigrational.test.cpp
• verify/verify-unit-test/composite-exp.test.cpp
• verify/verify-unit-test/composition.test.cpp
• verify/verify-unit-test/dual-fps.test.cpp
• verify/verify-unit-test/fft2d.test.cpp
• verify/verify-unit-test/fps-sparse.test.cpp
• verify/verify-unit-test/fps.test.cpp
• verify/verify-unit-test/garner-bigint.test.cpp
• verify/verify-unit-test/multieval.test.cpp
• verify/verify-unit-test/multipoint-binomial-sum.test.cpp
• verify/verify-unit-test/p-recursive.test.cpp
• verify/verify-unit-test/partial-fraction-decomposition.test.cpp
• verify/verify-unit-test/polynomial-matrix-prod.test.cpp
• verify/verify-unit-test/relaxed-convolution.test.cpp
• verify/verify-unit-test/stirling-matrix.test.cpp
• verify/verify-yosupo-fps/yosupo-composition-fast.test.cpp
• verify/verify-yosupo-fps/yosupo-composition-large.test.cpp
• verify/verify-yosupo-fps/yosupo-composition.test.cpp
• verify/verify-yosupo-fps/yosupo-compositional-inverse-large.test.cpp
• verify/verify-yosupo-fps/yosupo-compositional-inverse-newton.test.cpp
• verify/verify-yosupo-fps/yosupo-compositional-inverse.test.cpp
• verify/verify-yosupo-fps/yosupo-division-of-polynomials.test.cpp
• verify/verify-yosupo-fps/yosupo-exp-arb.test.cpp
• verify/verify-yosupo-fps/yosupo-exp-newton-method-2.test.cpp
• verify/verify-yosupo-fps/yosupo-exp-newton-method.test.cpp
• verify/verify-yosupo-fps/yosupo-exp-ofps.test.cpp
• verify/verify-yosupo-fps/yosupo-exp-relaxed-convolution.test.cpp
• verify/verify-yosupo-fps/yosupo-exp.test.cpp
• verify/verify-yosupo-fps/yosupo-factorial-p-recursive.test.cpp
• verify/verify-yosupo-fps/yosupo-factorial.test.cpp
• verify/verify-yosupo-fps/yosupo-interpolation.test.cpp
• verify/verify-yosupo-fps/yosupo-inv-arb.test.cpp
• verify/verify-yosupo-fps/yosupo-inv-newton-method.test.cpp
• verify/verify-yosupo-fps/yosupo-inv-of-polynomials.test.cpp
• verify/verify-yosupo-fps/yosupo-inv-ofps.test.cpp
• verify/verify-yosupo-fps/yosupo-inv-relaxed-convolution.test.cpp
• verify/verify-yosupo-fps/yosupo-inv.test.cpp
• verify/verify-yosupo-fps/yosupo-log-arb.test.cpp
• verify/verify-yosupo-fps/yosupo-log.test.cpp
• verify/verify-yosupo-fps/yosupo-multieval-fast.test.cpp
• verify/verify-yosupo-fps/yosupo-multieval.test.cpp
• verify/verify-yosupo-fps/yosupo-polynomial-root-finding.test.cpp
• verify/verify-yosupo-fps/yosupo-pow-arb.test.cpp
• verify/verify-yosupo-fps/yosupo-pow.test.cpp
• verify/verify-yosupo-fps/yosupo-product-of-polynomial-sequence.test.cpp
• verify/verify-yosupo-fps/yosupo-sample-point-shift.test.cpp
• verify/verify-yosupo-fps/yosupo-sparse-exp.test.cpp
• verify/verify-yosupo-fps/yosupo-sparse-inv.test.cpp
• verify/verify-yosupo-fps/yosupo-sparse-log.test.cpp
• verify/verify-yosupo-fps/yosupo-sparse-pow.test.cpp
• verify/verify-yosupo-fps/yosupo-sqrt.test.cpp
• verify/verify-yosupo-fps/yosupo-stirling-1st-row.test.cpp
• verify/verify-yosupo-fps/yosupo-stirling-1st.test.cpp
• verify/verify-yosupo-fps/yosupo-stirling-2nd-row.test.cpp
• verify/verify-yosupo-fps/yosupo-stirling-2nd.test.cpp
• verify/verify-yosupo-fps/yosupo-taylor-shift.test.cpp
• verify/verify-yosupo-math/yosupo-addition-of-big-integers.test.cpp
• verify/verify-yosupo-math/yosupo-addition-of-hex.test.cpp
• verify/verify-yosupo-math/yosupo-determinant-of-matrix-bbla.test.cpp
• verify/verify-yosupo-math/yosupo-determinant-of-sparse-matrix-bbla.test.cpp
• verify/verify-yosupo-math/yosupo-division-of-big-integers.test.cpp
• verify/verify-yosupo-math/yosupo-division-of-hex.test.cpp
• verify/verify-yosupo-math/yosupo-multiplication-of-big-integers.test.cpp
• verify/verify-yosupo-math/yosupo-multiplication-of-hex.test.cpp
• verify/verify-yosupo-math/yosupo-polynomial-composite-set-power-series.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-arbitrarylengthntt.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-arbitraryntt-arbitrarymodint.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-arbitraryntt-arbitraryprimemodint.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-arbitraryntt.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-chirp-z.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-large.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-ntt-normalmodint.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-ntt.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-on-z-pz.test.cpp
• verify/verify-yosupo-ntt/yosupo-convolution-relaxed-convolution.test.cpp
• verify/verify-yosupo-ntt/yosupo-multiplicative-convolution.test.cpp
• verify/verify-yosupo-ntt/yosupo-multipoint-evaluation-chirp-z.test.cpp
• verify/verify-yosupo-ntt/yosupo-multivariate-circular-convolution.test.cpp
• verify/verify-yosupo-other/yosupo-a-plus-b-128bit-bigint.test.cpp
• verify/verify-yosupo-string/yosupo-wildcard-pattern-matching.test.cpp
• verify/verify-yuki/yuki-0214.test.cpp
• verify/verify-yuki/yuki-0215-nth-term.test.cpp
• verify/verify-yuki/yuki-0215.test.cpp
• verify/verify-yuki/yuki-0303.test.cpp
• verify/verify-yuki/yuki-0502.test.cpp
• verify/verify-yuki/yuki-0720.test.cpp
• verify/verify-yuki/yuki-0963-circular.test.cpp
• verify/verify-yuki/yuki-0963.test.cpp
• verify/verify-yuki/yuki-1080.test.cpp
• verify/verify-yuki/yuki-1112-sparse.test.cpp
• verify/verify-yuki/yuki-1112.test.cpp
• verify/verify-yuki/yuki-1145-frac.test.cpp
• verify/verify-yuki/yuki-1145-sum-of-rational-functions.test.cpp
• verify/verify-yuki/yuki-1145.test.cpp
• verify/verify-yuki/yuki-1303.test.cpp
• verify/verify-yuki/yuki-1504.test.cpp
• verify/verify-yuki/yuki-1510.test.cpp
• verify/verify-yuki/yuki-1533.test.cpp
• verify/verify-yuki/yuki-1783.test.cpp
• verify/verify-yuki/yuki-1875.test.cpp
• verify/verify-yuki/yuki-1939-2.test.cpp
• verify/verify-yuki/yuki-1939-sparse-pow.test.cpp
• verify/verify-yuki/yuki-1939.test.cpp
• verify/verify-yuki/yuki-2231.test.cpp
• verify/verify-yuki/yuki-2580.test.cpp
• verify/verify-yuki/yuki-2661.test.cpp

Code

#pragma once

#include <algorithm>
#include <cstdint>
#include <iterator>
#include <vector>
using namespace std;

template <typename mint>
struct NTT {
  static constexpr uint32_t get_pr() {
    uint32_t _mod = mint::get_mod();
    using u64 = uint64_t;
    u64 ds[32] = {};
    int idx = 0;
    u64 m = _mod - 1;
    for (u64 i = 2; i * i <= m; ++i) {
      if (m % i == 0) {
        ds[idx++] = i;
        while (m % i == 0) m /= i;
      }
    }
    if (m != 1) ds[idx++] = m;

    uint32_t _pr = 2;
    while (1) {
      int flg = 1;
      for (int i = 0; i < idx; ++i) {
        u64 a = _pr, b = (_mod - 1) / ds[i], r = 1;
        while (b) {
          if (b & 1) r = r * a % _mod;
          a = a * a % _mod;
          b >>= 1;
        }
        if (r == 1) {
          flg = 0;
          break;
        }
      }
      if (flg == 1) break;
      ++_pr;
    }
    return _pr;
  };

  static constexpr uint32_t mod = mint::get_mod();
  static constexpr uint32_t pr = get_pr();
  static constexpr int level = __builtin_ctzll(mod - 1);
  mint dw[level], dy[level];

  void setwy(int k) {
    mint w[level], y[level];
    w[k - 1] = mint(pr).pow((mod - 1) / (1 << k));
    y[k - 1] = w[k - 1].inverse();
    for (int i = k - 2; i > 0; --i)
      w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1];
    dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2];
    for (int i = 3; i < k; ++i) {
      dw[i] = dw[i - 1] * y[i - 2] * w[i];
      dy[i] = dy[i - 1] * w[i - 2] * y[i];
    }
  }

  NTT() { setwy(level); }

  void fft4(vector<mint> &a, int k) {
    if ((int)a.size() <= 1) return;
    if (k == 1) {
      mint a1 = a[1];
      a[1] = a[0] - a[1];
      a[0] = a[0] + a1;
      return;
    }
    if (k & 1) {
      int v = 1 << (k - 1);
      for (int j = 0; j < v; ++j) {
        mint ajv = a[j + v];
        a[j + v] = a[j] - ajv;
        a[j] += ajv;
      }
    }
    int u = 1 << (2 + (k & 1));
    int v = 1 << (k - 2 - (k & 1));
    mint one = mint(1);
    mint imag = dw[1];
    while (v) {
      // jh = 0
      {
        int j0 = 0;
        int j1 = v;
        int j2 = j1 + v;
        int j3 = j2 + v;
        for (; j0 < v; ++j0, ++j1, ++j2, ++j3) {
          mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3];
          mint t0p2 = t0 + t2, t1p3 = t1 + t3;
          mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
          a[j0] = t0p2 + t1p3, a[j1] = t0p2 - t1p3;
          a[j2] = t0m2 + t1m3, a[j3] = t0m2 - t1m3;
        }
      }
      // jh >= 1
      mint ww = one, xx = one * dw[2], wx = one;
      for (int jh = 4; jh < u;) {
        ww = xx * xx, wx = ww * xx;
        int j0 = jh * v;
        int je = j0 + v;
        int j2 = je + v;
        for (; j0 < je; ++j0, ++j2) {
          mint t0 = a[j0], t1 = a[j0 + v] * xx, t2 = a[j2] * ww,
               t3 = a[j2 + v] * wx;
          mint t0p2 = t0 + t2, t1p3 = t1 + t3;
          mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
          a[j0] = t0p2 + t1p3, a[j0 + v] = t0p2 - t1p3;
          a[j2] = t0m2 + t1m3, a[j2 + v] = t0m2 - t1m3;
        }
        xx *= dw[__builtin_ctzll((jh += 4))];
      }
      u <<= 2;
      v >>= 2;
    }
  }

  void ifft4(vector<mint> &a, int k) {
    if ((int)a.size() <= 1) return;
    if (k == 1) {
      mint a1 = a[1];
      a[1] = a[0] - a[1];
      a[0] = a[0] + a1;
      return;
    }
    int u = 1 << (k - 2);
    int v = 1;
    mint one = mint(1);
    mint imag = dy[1];
    while (u) {
      // jh = 0
      {
        int j0 = 0;
        int j1 = v;
        int j2 = v + v;
        int j3 = j2 + v;
        for (; j0 < v; ++j0, ++j1, ++j2, ++j3) {
          mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3];
          mint t0p1 = t0 + t1, t2p3 = t2 + t3;
          mint t0m1 = t0 - t1, t2m3 = (t2 - t3) * imag;
          a[j0] = t0p1 + t2p3, a[j2] = t0p1 - t2p3;
          a[j1] = t0m1 + t2m3, a[j3] = t0m1 - t2m3;
        }
      }
      // jh >= 1
      mint ww = one, xx = one * dy[2], yy = one;
      u <<= 2;
      for (int jh = 4; jh < u;) {
        ww = xx * xx, yy = xx * imag;
        int j0 = jh * v;
        int je = j0 + v;
        int j2 = je + v;
        for (; j0 < je; ++j0, ++j2) {
          mint t0 = a[j0], t1 = a[j0 + v], t2 = a[j2], t3 = a[j2 + v];
          mint t0p1 = t0 + t1, t2p3 = t2 + t3;
          mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy;
          a[j0] = t0p1 + t2p3, a[j2] = (t0p1 - t2p3) * ww;
          a[j0 + v] = t0m1 + t2m3, a[j2 + v] = (t0m1 - t2m3) * ww;
        }
        xx *= dy[__builtin_ctzll(jh += 4)];
      }
      u >>= 4;
      v <<= 2;
    }
    if (k & 1) {
      u = 1 << (k - 1);
      for (int j = 0; j < u; ++j) {
        mint ajv = a[j] - a[j + u];
        a[j] += a[j + u];
        a[j + u] = ajv;
      }
    }
  }

  void ntt(vector<mint> &a) {
    if ((int)a.size() <= 1) return;
    fft4(a, __builtin_ctz(a.size()));
  }

  void intt(vector<mint> &a) {
    if ((int)a.size() <= 1) return;
    ifft4(a, __builtin_ctz(a.size()));
    mint iv = mint(a.size()).inverse();
    for (auto &x : a) x *= iv;
  }

  vector<mint> multiply(const vector<mint> &a, const vector<mint> &b) {
    int l = a.size() + b.size() - 1;
    if (min<int>(a.size(), b.size()) <= 40) {
      vector<mint> s(l);
      for (int i = 0; i < (int)a.size(); ++i)
        for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j];
      return s;
    }
    int k = 2, M = 4;
    while (M < l) M <<= 1, ++k;
    setwy(k);
    vector<mint> s(M);
    for (int i = 0; i < (int)a.size(); ++i) s[i] = a[i];
    fft4(s, k);
    if (a.size() == b.size() && a == b) {
      for (int i = 0; i < M; ++i) s[i] *= s[i];
    } else {
      vector<mint> t(M);
      for (int i = 0; i < (int)b.size(); ++i) t[i] = b[i];
      fft4(t, k);
      for (int i = 0; i < M; ++i) s[i] *= t[i];
    }
    ifft4(s, k);
    s.resize(l);
    mint invm = mint(M).inverse();
    for (int i = 0; i < l; ++i) s[i] *= invm;
    return s;
  }

  void ntt_doubling(vector<mint> &a) {
    int M = (int)a.size();
    auto b = a;
    intt(b);
    mint r = 1, zeta = mint(pr).pow((mint::get_mod() - 1) / (M << 1));
    for (int i = 0; i < M; i++) b[i] *= r, r *= zeta;
    ntt(b);
    copy(begin(b), end(b), back_inserter(a));
  }
};

#line 2 "ntt/ntt.hpp"

#include <algorithm>
#include <cstdint>
#include <iterator>
#include <vector>
using namespace std;

template <typename mint>
struct NTT {
  static constexpr uint32_t get_pr() {
    uint32_t _mod = mint::get_mod();
    using u64 = uint64_t;
    u64 ds[32] = {};
    int idx = 0;
    u64 m = _mod - 1;
    for (u64 i = 2; i * i <= m; ++i) {
      if (m % i == 0) {
        ds[idx++] = i;
        while (m % i == 0) m /= i;
      }
    }
    if (m != 1) ds[idx++] = m;

    uint32_t _pr = 2;
    while (1) {
      int flg = 1;
      for (int i = 0; i < idx; ++i) {
        u64 a = _pr, b = (_mod - 1) / ds[i], r = 1;
        while (b) {
          if (b & 1) r = r * a % _mod;
          a = a * a % _mod;
          b >>= 1;
        }
        if (r == 1) {
          flg = 0;
          break;
        }
      }
      if (flg == 1) break;
      ++_pr;
    }
    return _pr;
  };

  static constexpr uint32_t mod = mint::get_mod();
  static constexpr uint32_t pr = get_pr();
  static constexpr int level = __builtin_ctzll(mod - 1);
  mint dw[level], dy[level];

  void setwy(int k) {
    mint w[level], y[level];
    w[k - 1] = mint(pr).pow((mod - 1) / (1 << k));
    y[k - 1] = w[k - 1].inverse();
    for (int i = k - 2; i > 0; --i)
      w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1];
    dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2];
    for (int i = 3; i < k; ++i) {
      dw[i] = dw[i - 1] * y[i - 2] * w[i];
      dy[i] = dy[i - 1] * w[i - 2] * y[i];
    }
  }

  NTT() { setwy(level); }

  void fft4(vector<mint> &a, int k) {
    if ((int)a.size() <= 1) return;
    if (k == 1) {
      mint a1 = a[1];
      a[1] = a[0] - a[1];
      a[0] = a[0] + a1;
      return;
    }
    if (k & 1) {
      int v = 1 << (k - 1);
      for (int j = 0; j < v; ++j) {
        mint ajv = a[j + v];
        a[j + v] = a[j] - ajv;
        a[j] += ajv;
      }
    }
    int u = 1 << (2 + (k & 1));
    int v = 1 << (k - 2 - (k & 1));
    mint one = mint(1);
    mint imag = dw[1];
    while (v) {
      // jh = 0
      {
        int j0 = 0;
        int j1 = v;
        int j2 = j1 + v;
        int j3 = j2 + v;
        for (; j0 < v; ++j0, ++j1, ++j2, ++j3) {
          mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3];
          mint t0p2 = t0 + t2, t1p3 = t1 + t3;
          mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
          a[j0] = t0p2 + t1p3, a[j1] = t0p2 - t1p3;
          a[j2] = t0m2 + t1m3, a[j3] = t0m2 - t1m3;
        }
      }
      // jh >= 1
      mint ww = one, xx = one * dw[2], wx = one;
      for (int jh = 4; jh < u;) {
        ww = xx * xx, wx = ww * xx;
        int j0 = jh * v;
        int je = j0 + v;
        int j2 = je + v;
        for (; j0 < je; ++j0, ++j2) {
          mint t0 = a[j0], t1 = a[j0 + v] * xx, t2 = a[j2] * ww,
               t3 = a[j2 + v] * wx;
          mint t0p2 = t0 + t2, t1p3 = t1 + t3;
          mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
          a[j0] = t0p2 + t1p3, a[j0 + v] = t0p2 - t1p3;
          a[j2] = t0m2 + t1m3, a[j2 + v] = t0m2 - t1m3;
        }
        xx *= dw[__builtin_ctzll((jh += 4))];
      }
      u <<= 2;
      v >>= 2;
    }
  }

  void ifft4(vector<mint> &a, int k) {
    if ((int)a.size() <= 1) return;
    if (k == 1) {
      mint a1 = a[1];
      a[1] = a[0] - a[1];
      a[0] = a[0] + a1;
      return;
    }
    int u = 1 << (k - 2);
    int v = 1;
    mint one = mint(1);
    mint imag = dy[1];
    while (u) {
      // jh = 0
      {
        int j0 = 0;
        int j1 = v;
        int j2 = v + v;
        int j3 = j2 + v;
        for (; j0 < v; ++j0, ++j1, ++j2, ++j3) {
          mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3];
          mint t0p1 = t0 + t1, t2p3 = t2 + t3;
          mint t0m1 = t0 - t1, t2m3 = (t2 - t3) * imag;
          a[j0] = t0p1 + t2p3, a[j2] = t0p1 - t2p3;
          a[j1] = t0m1 + t2m3, a[j3] = t0m1 - t2m3;
        }
      }
      // jh >= 1
      mint ww = one, xx = one * dy[2], yy = one;
      u <<= 2;
      for (int jh = 4; jh < u;) {
        ww = xx * xx, yy = xx * imag;
        int j0 = jh * v;
        int je = j0 + v;
        int j2 = je + v;
        for (; j0 < je; ++j0, ++j2) {
          mint t0 = a[j0], t1 = a[j0 + v], t2 = a[j2], t3 = a[j2 + v];
          mint t0p1 = t0 + t1, t2p3 = t2 + t3;
          mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy;
          a[j0] = t0p1 + t2p3, a[j2] = (t0p1 - t2p3) * ww;
          a[j0 + v] = t0m1 + t2m3, a[j2 + v] = (t0m1 - t2m3) * ww;
        }
        xx *= dy[__builtin_ctzll(jh += 4)];
      }
      u >>= 4;
      v <<= 2;
    }
    if (k & 1) {
      u = 1 << (k - 1);
      for (int j = 0; j < u; ++j) {
        mint ajv = a[j] - a[j + u];
        a[j] += a[j + u];
        a[j + u] = ajv;
      }
    }
  }

  void ntt(vector<mint> &a) {
    if ((int)a.size() <= 1) return;
    fft4(a, __builtin_ctz(a.size()));
  }

  void intt(vector<mint> &a) {
    if ((int)a.size() <= 1) return;
    ifft4(a, __builtin_ctz(a.size()));
    mint iv = mint(a.size()).inverse();
    for (auto &x : a) x *= iv;
  }

  vector<mint> multiply(const vector<mint> &a, const vector<mint> &b) {
    int l = a.size() + b.size() - 1;
    if (min<int>(a.size(), b.size()) <= 40) {
      vector<mint> s(l);
      for (int i = 0; i < (int)a.size(); ++i)
        for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j];
      return s;
    }
    int k = 2, M = 4;
    while (M < l) M <<= 1, ++k;
    setwy(k);
    vector<mint> s(M);
    for (int i = 0; i < (int)a.size(); ++i) s[i] = a[i];
    fft4(s, k);
    if (a.size() == b.size() && a == b) {
      for (int i = 0; i < M; ++i) s[i] *= s[i];
    } else {
      vector<mint> t(M);
      for (int i = 0; i < (int)b.size(); ++i) t[i] = b[i];
      fft4(t, k);
      for (int i = 0; i < M; ++i) s[i] *= t[i];
    }
    ifft4(s, k);
    s.resize(l);
    mint invm = mint(M).inverse();
    for (int i = 0; i < l; ++i) s[i] *= invm;
    return s;
  }

  void ntt_doubling(vector<mint> &a) {
    int M = (int)a.size();
    auto b = a;
    intt(b);
    mint r = 1, zeta = mint(pr).pow((mint::get_mod() - 1) / (M << 1));
    for (int i = 0; i < M; i++) b[i] *= r, r *= zeta;
    ntt(b);
    copy(begin(b), end(b), back_inserter(a));
  }
};

Back to top page