PolynomialRing.cpp

• View this file on GitHub
• Last update: 2023-05-10 23:48:48-03:00
• Link: View error logs on GitHub Actions

Depends on

• Epsilon.cpp
• LongMultiplication.cpp
• Math.cpp
• ModularInteger.cpp
• NumberTheory.cpp
• Traits.cpp
• VectorN.cpp
• bits/stdc++.h

Required by

• LinearRecurrence.cpp
• PowerSeries.cpp
• polynomial/MultipointEvaluation.cpp

Verified with

• tests/yosupo/find-lr.test.cpp
• tests/yosupo/find-lr.test.cpp
• tests/yosupo/fps-inv.test.cpp
• tests/yosupo/fps-inv.test.cpp
• tests/yosupo/fps-power.test.cpp
• tests/yosupo/fps-power.test.cpp
• tests/yosupo/kth-term-lr.test.cpp
• tests/yosupo/kth-term-lr.test.cpp
• tests/yosupo/multipoint.test.cpp
• tests/yosupo/multipoint.test.cpp
• tests/yosupo/subset-sum.test.cpp
• tests/yosupo/subset-sum.test.cpp

Code

#ifndef _LIB_POLYNOMIAL_RING
#define _LIB_POLYNOMIAL_RING
#include "Epsilon.cpp"
#include "Math.cpp"
#include "ModularInteger.cpp"
#include "Traits.cpp"
#include "LongMultiplication.cpp"
#include "VectorN.cpp"
#include <bits/stdc++.h>

namespace lib {
using namespace std;
namespace math {
namespace poly {

namespace {
/// keep caide
using traits::IsInputIterator;
/// keep caide
using traits::HasInputIterator;
} // namespace

namespace detail {
  template<class>
  struct sfinae_true : std::true_type{};

  template<class T, class Field, class Func>
  static auto test_transform(int)
      -> sfinae_true<decltype(
    std::declval<T>().template on_transform<Field>(std::declval<int>(), std::declval<Func&>()))>;

  template<class, class Field, class Func>
  static auto test_transform(long) -> std::false_type;
} // detail::

template<class T, class Field, class Func = std::function<VectorN<Field>(int)>>
struct has_transform : decltype(detail::test_transform<T, Field, Func>(0)){};

template <typename P> struct DefaultPowerOp {
  int mod;
  DefaultPowerOp(int mod) : mod(mod) {}
  P operator()() const { return P(1); }
  P operator()(const P &a) const { return a % mod; }
  void operator()(P &x, const P &a, long long cur) const {
    (x *= x) %= mod;
    if (cur & 1)
      (x *= a) %= mod;
  }
};

template <typename P> struct ModPowerOp {
  const P &mod;
  ModPowerOp(const P &p) : mod(p) {}
  P operator()() const { return P(1); }
  P operator()(const P &a) { return a % mod; }
  void operator()(P &x, const P &a, long long cur) const {
    (x *= x) %= mod;
    if (cur & 1)
      (x *= a) %= mod;
  }
};

template <typename P> struct ModShiftPowerOp {
  const P &mod;
  ModShiftPowerOp(const P &p) : mod(p) {}
  P operator()() const { return P(1); }
  P operator()(const P &a) { return a % mod; }
  void operator()(P &x, const P &a, long long cur) const {
    // if(cur < mod.degree())
    // x = P::kth(cur);
    if (cur & 1)
      (x *= (x << 1)) %= mod;
    else
      (x *= x) %= mod;
  }
};

struct DefaultDivmod;
struct NaiveDivmod;

template <typename Field, typename Mult, typename Divmod = DefaultDivmod>
struct Polynomial {
  constexpr static int Magic = 64;
  constexpr static bool NaiveMod = is_same<Divmod, NaiveDivmod>::value;
  constexpr static bool HasTransform = has_transform<Mult, Field>::value;
  using Transform =  typename Mult::template Transform<Field>;

  typedef Polynomial<Field, Mult, Divmod> type;
  typedef Field field;
  vector<Field> p;

  Polynomial() : p(0) {}
  explicit Polynomial(Field x) : p(1, x) {}

  template <
      typename Iterator,
      typename enable_if<IsInputIterator<Iterator>::value>::type * = nullptr>
  Polynomial(Iterator begin, Iterator end) : p(distance(begin, end)) {
    int i = 0;
    for (auto it = begin; it != end; ++it, ++i)
      p[i] = *it;
    normalize();
  }

  template <
      typename Container,
      typename enable_if<HasInputIterator<Container>::value>::type * = nullptr>
  Polynomial(const Container &container)
      : Polynomial(container.begin(), container.end()) {}

  Polynomial(const initializer_list<Field> &v)
      : Polynomial(v.begin(), v.end()) {}

  static type from_root(const Field &root) { return Polynomial({-root, 1}); }

  void normalize() const {
    type *self = const_cast<type *>(this);
    int sz = self->p.size();
    while (sz > 0 && Epsilon<>().null(self->p[sz - 1]))
      sz--;
    if (sz != (int)self->p.size())
      self->p.resize(sz);
  }

  inline int size() const { return p.size(); }
  inline int degree() const { return max((int)p.size() - 1, 0); }
  bool null() const {
    for (Field x : p)
      if (!Epsilon<>().null(x))
        return false;
    return true;
  }

  const vector<Field>& data() const {
    return p;
  }

  Field eval(Field x) const {
    Field pw = 1;
    Field res = 0;
    for(Field c : p) {
      res += pw * c;
      pw *= x;
    }
    return res;
  }

  inline Field operator[](const int i) const {
    if (i >= size())
      return 0;
    return p[i];
  }
  inline Field &operator[](const int i) {
    if (i >= size())
      p.resize(i + 1);
    return p[i];
  }

  Field operator()(const Field &x) const {
    if (null())
      return Field();
    Field acc = p.back();
    for (int i = (int)size() - 2; i >= 0; i--) {
      acc *= x;
      acc += p[i];
    }
    return acc;
  }

  type substr(int i, int sz) const {
    int j = min(sz + i, size());
    i = min(i, size());
    if(i >= j) return type();
    return type(begin(p)+i, begin(p)+j);
  }

  type &operator+=(const type &rhs) {
    if (rhs.size() > size())
      p.resize(rhs.size());
    int sz = rhs.size();
    for (int i = 0; i < sz; i++)
      p[i] += rhs[i];
    normalize();
    return *this;
  }

  type &operator-=(const type &rhs) {
    if (rhs.size() > size())
      p.resize(rhs.size());
    int sz = rhs.size();
    for (int i = 0; i < sz; i++)
      p[i] -= rhs[i];
    normalize();
    return *this;
  }

  static vector<Field> multiply(const vector<Field>& a, const vector<Field>& b) {
    if(min(a.size(), b.size()) < Magic)
      return NaiveMultiplication()(a, b);
    return Mult()(a, b);
  }

  type &operator*=(const type &rhs) {
    p = multiply(p, rhs.p);
    normalize();
    return *this;
  }

  type &operator*=(const Field &rhs) {
    int sz = size();
    for (int i = 0; i < sz; i++)
      p[i] *= rhs;
    normalize();
    return *this;
  }

  type &operator/=(const Field &rhs) {
    int sz = size();
    for (int i = 0; i < sz; i++)
      p[i] /= rhs;
    normalize();
    return *this;
  }

  type &operator<<=(const int rhs) {
    if (rhs < 0)
      return *this >>= rhs;
    if (rhs == 0)
      return *this;
    int sz = size();
    p.resize(sz + rhs);
    for (int i = sz - 1; i >= 0; i--)
      p[i + rhs] = p[i];
    fill_n(p.begin(), rhs, 0);
    return *this;
  }

  type &operator>>=(const int rhs) {
    if (rhs < 0)
      return *this <<= rhs;
    if (rhs == 0)
      return *this;
    int sz = size();
    if (rhs >= sz) {
      p.clear();
      return *this;
    }
    for (int i = rhs; i < sz; i++)
      p[i - rhs] = p[i];
    p.resize(sz - rhs);
    return *this;
  }

  type &operator%=(const int rhs) {
    if (rhs < size())
      p.resize(rhs);
    normalize();
    return *this;
  }

  type &operator/=(const type &rhs) { return *this = *this / rhs; }

  type operator%=(const type &rhs) { return *this = *this % rhs; }

  type operator+(const type &rhs) const {
    type res = *this;
    return res += rhs;
  }

  type operator-(const type &rhs) const {
    type res = *this;
    return res -= rhs;
  }

  type operator*(const type &rhs) const { 
    type res(multiply(p, rhs.p));
    res.normalize();
    return res;
  }

  type operator*(const Field &rhs) const {
    type res = *this;
    return res *= rhs;
  }

  type operator/(const Field &rhs) const {
    type res = *this;
    return res /= rhs;
  }

  type operator<<(const int rhs) const {
    type res = *this;
    return res <<= rhs;
  }

  type operator>>(const int rhs) const {
    type res = *this;
    return res >>= rhs;
  }

  type operator%(const int rhs) const {
    return Polynomial(p.begin(), p.begin() + min(rhs, size()));
  }

  type operator/(const type &rhs) const {
    return type::divmod(*this, rhs).first;
  }

  type operator%(const type &rhs) const {
    return type::divmod(*this, rhs).second;
  }

  bool operator==(const type &rhs) const {
    normalize();
    rhs.normalize();
    return p == rhs.p;
  }

  template <// Used in SFINAE.
            typename U = Field,
            enable_if_t<has_transform<Mult, U>::value>* = nullptr>
  inline VectorN<U> transform(int n) {
    return Mult().template transform<U>(n, p);
  }

  template <// Used in SFINAE.
            typename U = Field,
            enable_if_t<has_transform<Mult, U>::value>* = nullptr>
  inline static type itransform(int n, const vector<U>& v) {
    return Mult().template itransform<U>(n, v);
  }

  template <typename Functor,
            // Used in SFINAE.
            typename U = Field,
            enable_if_t<has_transform<Mult, U>::value>* = nullptr,
            typename ...Ts>
  inline static type on_transform(
    int n,
    Functor f,        
    const Ts&... vs) {
    if(n < Magic)
      return f(n, vs...);
    return Mult().template on_transform<U>(n, f, vs.p...);
  }

  template <typename Functor,
            // Used in SFINAE.
            typename U = Field,
            enable_if_t<!has_transform<Mult, U>::value>* = nullptr,
            typename ...Ts>
  inline static type on_transform(
    int n,
    Functor f,        
    const Ts&... vs) {
    return f(n, vs...);
  }

  template <
    // Used in SFINAE.
    typename U = Field,
    enable_if_t<has_transform<Mult, U>::value>* = nullptr>
  type inverse(int m) const {
    if(null()) return *this;
    type r = {Field(1) / p[0]};
    r.p.reserve(m);
    for(int i = 1; i < m; i *= 2) {
      int n = 2 * i;
      vector<U> f = (*this % n).p; f.resize(n);
      vector<U> g = r.p; g.resize(n);
      Transform::dft(f, n);
      Transform::dft(g, n);
      for(int j = 0; j < n; j++) f[j] *= g[j];
      Transform::idft(f, n);
      for(int j = 0; j < i; j++) f[j] = 0;
      Transform::dft(f, n);
      for(int j = 0; j < n; j++) f[j] *= g[j];
      Transform::idft(f, n);
      for(int j = i; j < min(n, m); j++)
        r[j] = -f[j];
    }
    return r;
  }

  type inverse_slow(int m) const {
    if(null()) return *this;
    type b = {Field(1) / p[0]};
    b.p.reserve(2 * m);
    for(int i = 1; i < m; i *= 2) {
      int n = min(2 * i, m);
      auto bb = b * b % n;
      b += b;
      b -= *this % n * bb;
      b %= n;
    }
    return b % m;
  }

  template <
    // Used in SFINAE.
    typename U = Field,
    enable_if_t<!has_transform<Mult, U>::value>* = nullptr>
  type inverse(int m) const {
    return inverse_slow(m);
  }

  type inverse() const {
    return inverse(size());
  }

  type reciprocal() const {
    normalize();
    return type(p.rbegin(), p.rend());
  }

  type integral() const {
    int sz = size();
    if (sz == 0)
      return {};
    type res = *this;
    for (int i = sz; i; i--) {
      res[i] = res[i - 1] / i;
    }
    res[0] = 0;
    res.normalize();
    return res;
  }

  type derivative() const {
    int sz = size();
    if (sz == 0)
      return {};
    type res = *this;
    for (int i = 0; i + 1 < sz; i++) {
      res[i] = res[i + 1] * (i + 1);
    }
    res.p.back() = 0;
    res.normalize();
    return res;
  }

  type mulx(field x) const { // component-wise multiplication with x^k
    field cur = 1;
    type res(*this);
    for(auto& c : res.p)
      c *= cur, cur *= x;
    return res;
  }
  type mulx_sq(field x) const { // component-wise multiplication with x^{k^2}
    field cur = x;
    field total = 1;
    field xx = x * x;
    type res(*this);
    for(auto& c : res.p)
      c *= total, total *= cur, cur *= xx;
    return res;
  }
  static pair<type, type> divmod(const type &a, const type &b) {
    if (NaiveMod || min(a.size(), b.size()) < Magic)
      return naive_divmod(a, b);
    a.normalize();
    b.normalize();
    int m = a.size();
    int n = b.size();
    if (m < n)
      return {Polynomial(), a};
    int sz = m - n + 1;
    type ar = a.reciprocal() % sz;
    type br = b.reciprocal() % sz;
    type q = (ar * br.inverse(sz) % sz).reciprocal();
    type r = a - b * q;

    return {q, r % (n-1)};
  }

  static pair<type, type> naive_divmod(const type &a, const type &b) {
    type res = a;
    int a_deg = a.degree();
    int b_deg = b.degree();
    Field normalizer = Field(1) / b[b_deg];
    for (int i = 0; i < a_deg - b_deg + 1; i++) {
      Field coef = (res[a_deg - i] *= normalizer);
      if (coef != 0) {
        for (int j = 1; j <= b_deg; j++) {
          res[a_deg - i - j] += -b[b_deg - j] * coef;
        }
      }
    }
    return {res >> b_deg, res % b_deg};
  }
  vector<Field> czt_even(Field z, int n) const { // P(1), P(z^2), P(z^4), ..., P(z^2(n-1))
    int m = degree();
    if(null()) {
      return vector<Field>(n);
    }
    vector<Field> vv(m + n);
    Field zi = Field(1) / z;
    Field zz = zi * zi;
    Field cur = zi;
    Field total = 1;
    for(int i = 0; i <= max(n - 1, m); i++) {
      if(i <= m) {vv[m - i] = total;}
      if(i < n) {vv[m + i] = total;}
      total *= cur;
      cur *= zz;
    }
    type w = (mulx_sq(z) * vv).substr(m, n).mulx_sq(z);
    vector<Field> res(n);
    for(int i = 0; i < n; i++) {
      res[i] = w[i];
    }
    return res;
  }
  vector<Field> czt(Field z, int n) const {
    auto even = czt_even(z, (n+1)/2);
    auto odd = mulx(z).czt_even(z, n/2);
    vector<Field> ans(n);
    for(int i = 0; i < n/2; i++) {
      ans[2*i] = even[i];
      ans[2*i+1] = odd[i];
    }
    if(n&1) {
      ans.back() = even.back();
    }
    return ans;
  }
  friend type kmul(const vector<type>& polys, int l, int r) {
    if(l == r) return polys[l];
    int mid = (l+r)/2;
    return kmul(polys, l, mid) * kmul(polys, mid+1, r);
  }
  friend type kmul(const vector<type>& polys) {
    if(polys.empty()) return type();
    return kmul(polys, 0, (int)polys.size() - 1);
  }
  static type power(const type &a, long long n, const int mod) {
    return math::generic_power<type>(a, n, DefaultPowerOp<type>(mod));
  }

  static type power(const type &a, long long n, const type &mod) {
    return math::generic_power<type>(a, n, ModPowerOp<type>(mod));
  }

  static type kth(int K) { return type(1) << K; }

  static type kth(long long K, const type &mod) {
    return math::generic_power<type>(type(1) << 1, K,
                                     ModShiftPowerOp<type>(mod));
  }

  friend ostream &operator<<(ostream &output, const type &var) {
    output << "[";
    int sz = var.size();
    for (int i = sz - 1; i >= 0; i--) {
      output << var[i];
      if (i)
        output << " ";
    }
    return output << "]";
  }
};
} // namespace poly
/// keep caide
using poly::Polynomial;
} // namespace math
} // namespace lib

#endif

#line 1 "PolynomialRing.cpp"


#line 1 "Epsilon.cpp"


#include <bits/stdc++.h>

namespace lib {
using namespace std;

template <typename T = double> struct Epsilon {
  T eps;
  constexpr Epsilon(T eps = 1e-9) : eps(eps) {}

  template <typename G,
            typename enable_if<is_floating_point<G>::value>::type * = nullptr>
  int operator()(G a, G b = 0) const {
    return a + eps < b ? -1 : (b + eps < a ? 1 : 0);
  }

  template <typename G,
            typename enable_if<!is_floating_point<G>::value>::type * = nullptr>
  int operator()(G a, G b = 0) const {
    return a < b ? -1 : (a > b ? 1 : 0);
  }

  template <typename G,
            typename enable_if<is_floating_point<G>::value>::type * = nullptr>
  bool null(G a) const {
    return (*this)(a) == 0;
  }

  template <typename G,
            typename enable_if<!is_floating_point<G>::value>::type * = nullptr>
  bool null(G a) const {
    return a == 0;
  }
};
} // namespace lib


#line 1 "Math.cpp"


#line 4 "Math.cpp"

namespace lib {
using namespace std;
namespace math {

/// caide keep
template <typename Type> struct DefaultPowerOp {
  Type operator()() const { return Type(1); }
  Type operator()(const Type &a) const { return a; }
  void operator()(Type &x, const Type &a, long long cur) const {
    x *= x;
    if (cur & 1)
      x *= a;
  }
};

template <typename Type, typename Op>
Type generic_power(const Type &a, long long n, Op op) {
  if (n == 0)
    return op();
  Type res = op(a);
  int hi = 63 - __builtin_clzll(n);
  for (int i = hi - 1; ~i; i--) {
    op(res, a, n >> i);
  }
  return res;
}

template <typename Type> Type generic_power(const Type &a, long long n) {
  return generic_power(a, n, DefaultPowerOp<Type>());
}
} // namespace math
} // namespace lib


#line 1 "ModularInteger.cpp"


#line 1 "NumberTheory.cpp"


#line 4 "NumberTheory.cpp"

namespace lib {
using namespace std;
namespace nt {
int64_t inverse(int64_t a, int64_t b) {
  long long b0 = b, t, q;
  long long x0 = 0, x1 = 1;
  if (b == 1)
    return 1;
  while (a > 1) {
    q = a / b;
    t = b, b = a % b, a = t;
    t = x0, x0 = x1 - q * x0, x1 = t;
  }
  if (x1 < 0)
    x1 += b0;
  return x1;
}
template<typename T, typename U>
T powmod (T a, U b, U p) {
    int res = 1;
    while (b)
        if (b & 1)
            res = (int) (res * 1ll * a % p),  --b;
        else
            a = (int) (a * 1ll * a % p),  b >>= 1;
    return res;
}
template<typename T>
vector<T> factors(T n) {
  vector<T> f;
  for(T i = 2; i*i <= n; i++) {
    if(n % i == 0) f.push_back(i);
    while(n % i == 0) n /= i;
  }
  if(n > 1) f.push_back(n);
  return f;
}
} // namespace nt
} // namespace lib


#line 5 "ModularInteger.cpp"

#if __cplusplus < 201300
#error required(c++14)
#endif

namespace lib {
using namespace std;
namespace {
template <typename T, T... Mods> struct ModularIntegerBase {
  typedef ModularIntegerBase<T, Mods...> type;

  T x[sizeof...(Mods)];
  friend ostream &operator<<(ostream &output, const type &var) {
    output << "(";
    for (int i = 0; i < sizeof...(Mods); i++) {
      if (i)
        output << ", ";
      output << var.x[i];
    }
    return output << ")";
  }
};

template <typename T, T Mod> struct ModularIntegerBase<T, Mod> {
  typedef ModularIntegerBase<T, Mod> type;
  constexpr static T mod = Mod;

  T x[1];

  T& data() { return this->x[0]; }
  T data() const { return this->x[0]; }
  explicit operator int() const { return this->x[0]; }
  explicit operator int64_t() const { return this->x[0]; }
  explicit operator double() const { return this->x[0]; }
  explicit operator long double() const { return this->x[0]; }
  friend ostream &operator<<(ostream &output, const type &var) {
    return output << var.x[0];
  }
};

template<typename T, typename U, T... Mods>
struct InversesTable {
  constexpr static size_t n_mods = sizeof...(Mods);
  constexpr static T mods[sizeof...(Mods)] = {Mods...};
  constexpr static int n_inverses = 1e6 + 10;

  T v[n_inverses][n_mods];
  T max_x;

  InversesTable() : v(), max_x(n_inverses) {
    for(int j = 0; j < sizeof...(Mods); j++)
      v[1][j] = 1, max_x = min(max_x, mods[j]);
    for(int i = 2; i < max_x; i++) {
      for(int j = 0; j < sizeof...(Mods); j++) {
        v[i][j] = mods[j] - (T)((U)(mods[j] / i) * v[mods[j] % i][j] % mods[j]);
      }
    }
  }
};

// Make available for linkage.
template <typename T, class U, T... Mods>
constexpr T InversesTable<T, U, Mods...>::mods[];

template <typename T, class Enable, T... Mods>
struct ModularIntegerImpl : ModularIntegerBase<T, Mods...> {
  typedef ModularIntegerImpl<T, Enable, Mods...> type;
  typedef T type_int;
  typedef uint64_t large_int;
  constexpr static size_t n_mods = sizeof...(Mods);
  constexpr static T mods[sizeof...(Mods)] = {Mods...};
  using ModularIntegerBase<T, Mods...>::x;
  using Inverses = InversesTable<T, large_int, Mods...>;

  struct Less {
    bool operator()(const type &lhs, const type &rhs) const {
      for (size_t i = 0; i < sizeof...(Mods); i++)
        if (lhs.x[i] != rhs.x[i])
          return lhs.x[i] < rhs.x[i];
      return false;
    };
  };
  typedef Less less;


  constexpr ModularIntegerImpl() {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      x[i] = T();
  }
  constexpr ModularIntegerImpl(large_int y) {
    for (size_t i = 0; i < sizeof...(Mods); i++) {
      x[i] = y % mods[i];
      if (x[i] < 0)
        x[i] += mods[i];
    }
  }
  static type with_remainders(T y[sizeof...(Mods)]) {
    type res;
    for (size_t i = 0; i < sizeof...(Mods); i++)
      res.x[i] = y[i];
    res.normalize();
    return res;
  }

  inline void normalize() {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((x[i] %= mods[i]) < 0)
        x[i] += mods[i];
  }

  inline T operator[](int i) const { return x[i]; }

  inline T multiply(T a, T b, T mod) const { return (large_int)a * b % mod; }

  inline T inv(T a, T mod) const { return static_cast<T>(nt::inverse(a, mod)); }

  inline T invi(T a, int i) const {
    const static Inverses inverses = Inverses();
    if(a < inverses.max_x)
      return inverses.v[a][i];
    return inv(a, mods[i]);
  }

  type inverse() const {
    T res[sizeof...(Mods)];
    for (size_t i = 0; i < sizeof...(Mods); i++)
      res[i] = invi(x[i], i);
    return type::with_remainders(res);
  }

  template <typename U> T power_(T a, U p, T mod) {
    if (mod == 1)
      return T();
    if (p < 0) {
      if (a == 0)
        throw domain_error("0^p with negative p is invalid");
      p = -p;
      a = inv(a, mod);
    }
    if (p == 0)
      return T(1);
    if (p == 1)
      return a;
    T res = 1;
    while (p > 0) {
      if (p & 1)
        res = multiply(res, a, mod);
      p >>= 1;
      a = multiply(a, a, mod);
    }
    return res;
  }

  inline type &operator+=(const type &rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((x[i] += rhs.x[i]) >= mods[i])
        x[i] -= mods[i];
    return *this;
  }
  inline type &operator-=(const type &rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((x[i] -= rhs.x[i]) < 0)
        x[i] += mods[i];
    return *this;
  }
  inline type &operator*=(const type &rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      x[i] = multiply(x[i], rhs.x[i], mods[i]);
    return *this;
  }
  inline type &operator/=(const type &rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      x[i] = multiply(x[i], invi(rhs.x[i], i), mods[i]);
    return *this;
  }

  inline type &operator+=(T rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((x[i] += rhs) >= mods[i])
        x[i] -= mods[i];
    return *this;
  }

  type &operator-=(T rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((x[i] -= rhs) < 0)
        x[i] += mods[i];
    return *this;
  }

  type &operator*=(T rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      x[i] = multiply(x[i], rhs, mods[i]);
    return *this;
  }

  type &operator/=(T rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      x[i] = multiply(invi(rhs, i), x[i], mods[i]);
    return *this;
  }

  type &operator^=(large_int p) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      x[i] = power_(x[i], p, mods[i]);
    return *this;
  }

  type &operator++() {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((++x[i]) >= mods[i])
        x[i] -= mods[i];
    return *this;
  }
  type &operator--() {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if ((--x[i]) < 0)
        x[i] += mods[i];
    return *this;
  }
  type operator++(int unused) {
    type res = *this;
    ++(*this);
    return res;
  }
  type operator--(int unused) {
    type res = *this;
    --(*this);
    return res;
  }

  friend type operator+(const type &lhs, const type &rhs) {
    type res = lhs;
    return res += rhs;
  }
  friend type operator-(const type &lhs, const type &rhs) {
    type res = lhs;
    return res -= rhs;
  }
  friend type operator*(const type &lhs, const type &rhs) {
    type res = lhs;
    return res *= rhs;
  }
  friend type operator/(const type &lhs, const type &rhs) {
    type res = lhs;
    return res /= rhs;
  }

  friend type operator+(const type &lhs, T rhs) {
    type res = lhs;
    return res += rhs;
  }

  friend type operator-(const type &lhs, T rhs) {
    type res = lhs;
    return res -= rhs;
  }

  friend type operator*(const type &lhs, T rhs) {
    type res = lhs;
    return res *= rhs;
  }

  friend type operator/(const type &lhs, T rhs) {
    type res = lhs;
    return res /= rhs;
  }

  friend type operator^(const type &lhs, large_int rhs) {
    type res = lhs;
    return res ^= rhs;
  }

  friend type power(const type &lhs, large_int rhs) { return lhs ^ rhs; }

  type operator-() const {
    type res = *this;
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if (res.x[i])
        res.x[i] = mods[i] - res.x[i];
    return res;
  }

  friend bool operator==(const type &lhs, const type &rhs) {
    for (size_t i = 0; i < sizeof...(Mods); i++)
      if (lhs.x[i] != rhs.x[i])
        return false;
    return true;
  }
  friend bool operator!=(const type &lhs, const type &rhs) {
    return !(lhs == rhs);
  }

  friend istream &operator>>(istream &input, type &var) {
    T y;
    cin >> y;
    var = y;
    return input;
  }
};
} // namespace

// Explicitly make constexpr available for linkage.
template <typename T, class Enable, T... Mods>
constexpr T ModularIntegerImpl<T, Enable, Mods...>::mods[];

template <typename T, T... Mods>
using ModularInteger =
    ModularIntegerImpl<T, typename enable_if<is_integral<T>::value>::type,
                       Mods...>;

template <int32_t... Mods> using Mint32 = ModularInteger<int32_t, Mods...>;

template <int64_t... Mods> using Mint64 = ModularInteger<int64_t, Mods...>;

using MintP = Mint32<(int32_t)1e9+7>;
using MintNTT = Mint32<998244353>;
} // namespace lib


#line 1 "Traits.cpp"


#line 4 "Traits.cpp"

namespace lib {
using namespace std;
namespace traits {

template <typename...> struct make_void { using type = void; };

template <typename... T> using void_t = typename make_void<T...>::type;

/// keep caide
template <typename Iterator>
using IteratorCategory = typename iterator_traits<Iterator>::iterator_category;

/// keep caide
template <typename Container>
using IteratorCategoryOf = IteratorCategory<typename Container::iterator>;

/// keep caide
template <typename Iterator>
using IteratorValue = typename iterator_traits<Iterator>::value_type;

/// keep caide
template <typename Container>
using IteratorValueOf = IteratorValue<typename Container::iterator>;

/// keep caide
template <typename Iterator>
using IsRandomIterator =
    is_base_of<random_access_iterator_tag, IteratorCategory<Iterator>>;

/// keep caide
template <typename Iterator>
using IsInputIterator =
    is_base_of<input_iterator_tag, IteratorCategory<Iterator>>;

/// keep caide
template <typename Iterator>
using IsBidirectionalIterator =
    is_base_of<bidirectional_iterator_tag, IteratorCategory<Iterator>>;

/// keep caide
template <typename Container>
using HasRandomIterator =
    is_base_of<random_access_iterator_tag, IteratorCategoryOf<Container>>;

/// keep caide
template <typename Container>
using HasInputIterator =
    is_base_of<input_iterator_tag, IteratorCategoryOf<Container>>;

/// keep caide
template <typename Container>
using HasBidirectionalIterator =
    is_base_of<bidirectional_iterator_tag, IteratorCategoryOf<Container>>;
} // namespace traits
} // namespace lib


#line 1 "LongMultiplication.cpp"


#line 4 "LongMultiplication.cpp"

namespace lib {
using namespace std;
namespace math {
struct NaiveMultiplication {
  template<typename T>
  using Transform = void;

  template <typename Field>
  vector<Field> operator()(const vector<Field> &a,
                           const vector<Field> &b) const {
    vector<Field> res(a.size() + b.size());
    for (size_t i = 0; i < a.size(); i++) {
      for (size_t j = 0; j < b.size(); j++) {
        res[i + j] += a[i] * b[j];
      }
    }
    return res;
  }
};

template <typename Mult, typename Field>
vector<Field> shift_conv(const vector<Field> &a, vector<Field> b) {
  if (b.empty())
    return {};
  reverse(b.begin(), b.end());
  int n = a.size();
  int m = b.size();

  auto res = Mult()(a, b);
  return vector<Field>(res.begin() + m - 1, res.end());
}
} // namespace math
} // namespace lib


#line 1 "VectorN.cpp"


#line 5 "VectorN.cpp"

#define VEC_CONST_OP(op, typ) \
  type operator op(const typ rhs) const { \
    auto res = *this; \
    return res op##= rhs; \
  }

#define VEC_BIN_OP(op) \
  type& operator op##=(const type& rhs) { \
    if(rhs.size() > this->size()) \
      this->resize(rhs.size()); \
    int sz = this->size(); \
    for(int i = 0; i < (int)rhs.size(); i++) \
      (*this)[i] op##= rhs[i]; \
    for(int i = rhs.size(); i < sz; i++) \
      (*this)[i] op##= 0; \
    return *this; \
  } \
  VEC_CONST_OP(op, type)

#define VEC_SINGLE_OP(op, typ) \
  type& operator op##=(const typ rhs) { \
    for(auto& x : *this) \
      x op##= rhs; \
    return *this; \
  } \
  VEC_CONST_OP(op, typ)
  

namespace lib {
using namespace std;
template<typename T>
struct VectorN : vector<T> {
  using type = VectorN<T>;

  template <
      typename Container,
      typename enable_if<traits::HasInputIterator<Container>::value>::type * = nullptr>
  VectorN(const Container &container)
      : vector<T>(container.begin(), container.end()) {}

  VectorN(const initializer_list<T> &v)
      : vector<T>(v.begin(), v.end()) {}

  template<typename... Args>
  VectorN( Args&&... args ) 
      : vector<T>(std::forward<Args>(args)...) {}

  VEC_BIN_OP(+)
  VEC_BIN_OP(-)
  VEC_BIN_OP(*)

  VEC_SINGLE_OP(+, T&)
  VEC_SINGLE_OP(-, T&)
  VEC_SINGLE_OP(*, T&)
  VEC_SINGLE_OP(/, T&)
  VEC_SINGLE_OP(^, int64_t)

  type operator-() const {
    auto res = *this;
    for(auto& x : res) x = -x;
    return res;
  }

  type operator%(int n) const {
    // TODO: get rid of this
    // return *const_cast<type*>(this);
    return *this;
  }
};
} // namespace lib


#line 10 "PolynomialRing.cpp"

namespace lib {
using namespace std;
namespace math {
namespace poly {

namespace {
/// keep caide
using traits::IsInputIterator;
/// keep caide
using traits::HasInputIterator;
} // namespace

namespace detail {
  template<class>
  struct sfinae_true : std::true_type{};

  template<class T, class Field, class Func>
  static auto test_transform(int)
      -> sfinae_true<decltype(
    std::declval<T>().template on_transform<Field>(std::declval<int>(), std::declval<Func&>()))>;

  template<class, class Field, class Func>
  static auto test_transform(long) -> std::false_type;
} // detail::

template<class T, class Field, class Func = std::function<VectorN<Field>(int)>>
struct has_transform : decltype(detail::test_transform<T, Field, Func>(0)){};

template <typename P> struct DefaultPowerOp {
  int mod;
  DefaultPowerOp(int mod) : mod(mod) {}
  P operator()() const { return P(1); }
  P operator()(const P &a) const { return a % mod; }
  void operator()(P &x, const P &a, long long cur) const {
    (x *= x) %= mod;
    if (cur & 1)
      (x *= a) %= mod;
  }
};

template <typename P> struct ModPowerOp {
  const P &mod;
  ModPowerOp(const P &p) : mod(p) {}
  P operator()() const { return P(1); }
  P operator()(const P &a) { return a % mod; }
  void operator()(P &x, const P &a, long long cur) const {
    (x *= x) %= mod;
    if (cur & 1)
      (x *= a) %= mod;
  }
};

template <typename P> struct ModShiftPowerOp {
  const P &mod;
  ModShiftPowerOp(const P &p) : mod(p) {}
  P operator()() const { return P(1); }
  P operator()(const P &a) { return a % mod; }
  void operator()(P &x, const P &a, long long cur) const {
    // if(cur < mod.degree())
    // x = P::kth(cur);
    if (cur & 1)
      (x *= (x << 1)) %= mod;
    else
      (x *= x) %= mod;
  }
};

struct DefaultDivmod;
struct NaiveDivmod;

template <typename Field, typename Mult, typename Divmod = DefaultDivmod>
struct Polynomial {
  constexpr static int Magic = 64;
  constexpr static bool NaiveMod = is_same<Divmod, NaiveDivmod>::value;
  constexpr static bool HasTransform = has_transform<Mult, Field>::value;
  using Transform =  typename Mult::template Transform<Field>;

  typedef Polynomial<Field, Mult, Divmod> type;
  typedef Field field;
  vector<Field> p;

  Polynomial() : p(0) {}
  explicit Polynomial(Field x) : p(1, x) {}

  template <
      typename Iterator,
      typename enable_if<IsInputIterator<Iterator>::value>::type * = nullptr>
  Polynomial(Iterator begin, Iterator end) : p(distance(begin, end)) {
    int i = 0;
    for (auto it = begin; it != end; ++it, ++i)
      p[i] = *it;
    normalize();
  }

  template <
      typename Container,
      typename enable_if<HasInputIterator<Container>::value>::type * = nullptr>
  Polynomial(const Container &container)
      : Polynomial(container.begin(), container.end()) {}

  Polynomial(const initializer_list<Field> &v)
      : Polynomial(v.begin(), v.end()) {}

  static type from_root(const Field &root) { return Polynomial({-root, 1}); }

  void normalize() const {
    type *self = const_cast<type *>(this);
    int sz = self->p.size();
    while (sz > 0 && Epsilon<>().null(self->p[sz - 1]))
      sz--;
    if (sz != (int)self->p.size())
      self->p.resize(sz);
  }

  inline int size() const { return p.size(); }
  inline int degree() const { return max((int)p.size() - 1, 0); }
  bool null() const {
    for (Field x : p)
      if (!Epsilon<>().null(x))
        return false;
    return true;
  }

  const vector<Field>& data() const {
    return p;
  }

  Field eval(Field x) const {
    Field pw = 1;
    Field res = 0;
    for(Field c : p) {
      res += pw * c;
      pw *= x;
    }
    return res;
  }

  inline Field operator[](const int i) const {
    if (i >= size())
      return 0;
    return p[i];
  }
  inline Field &operator[](const int i) {
    if (i >= size())
      p.resize(i + 1);
    return p[i];
  }

  Field operator()(const Field &x) const {
    if (null())
      return Field();
    Field acc = p.back();
    for (int i = (int)size() - 2; i >= 0; i--) {
      acc *= x;
      acc += p[i];
    }
    return acc;
  }

  type substr(int i, int sz) const {
    int j = min(sz + i, size());
    i = min(i, size());
    if(i >= j) return type();
    return type(begin(p)+i, begin(p)+j);
  }

  type &operator+=(const type &rhs) {
    if (rhs.size() > size())
      p.resize(rhs.size());
    int sz = rhs.size();
    for (int i = 0; i < sz; i++)
      p[i] += rhs[i];
    normalize();
    return *this;
  }

  type &operator-=(const type &rhs) {
    if (rhs.size() > size())
      p.resize(rhs.size());
    int sz = rhs.size();
    for (int i = 0; i < sz; i++)
      p[i] -= rhs[i];
    normalize();
    return *this;
  }

  static vector<Field> multiply(const vector<Field>& a, const vector<Field>& b) {
    if(min(a.size(), b.size()) < Magic)
      return NaiveMultiplication()(a, b);
    return Mult()(a, b);
  }

  type &operator*=(const type &rhs) {
    p = multiply(p, rhs.p);
    normalize();
    return *this;
  }

  type &operator*=(const Field &rhs) {
    int sz = size();
    for (int i = 0; i < sz; i++)
      p[i] *= rhs;
    normalize();
    return *this;
  }

  type &operator/=(const Field &rhs) {
    int sz = size();
    for (int i = 0; i < sz; i++)
      p[i] /= rhs;
    normalize();
    return *this;
  }

  type &operator<<=(const int rhs) {
    if (rhs < 0)
      return *this >>= rhs;
    if (rhs == 0)
      return *this;
    int sz = size();
    p.resize(sz + rhs);
    for (int i = sz - 1; i >= 0; i--)
      p[i + rhs] = p[i];
    fill_n(p.begin(), rhs, 0);
    return *this;
  }

  type &operator>>=(const int rhs) {
    if (rhs < 0)
      return *this <<= rhs;
    if (rhs == 0)
      return *this;
    int sz = size();
    if (rhs >= sz) {
      p.clear();
      return *this;
    }
    for (int i = rhs; i < sz; i++)
      p[i - rhs] = p[i];
    p.resize(sz - rhs);
    return *this;
  }

  type &operator%=(const int rhs) {
    if (rhs < size())
      p.resize(rhs);
    normalize();
    return *this;
  }

  type &operator/=(const type &rhs) { return *this = *this / rhs; }

  type operator%=(const type &rhs) { return *this = *this % rhs; }

  type operator+(const type &rhs) const {
    type res = *this;
    return res += rhs;
  }

  type operator-(const type &rhs) const {
    type res = *this;
    return res -= rhs;
  }

  type operator*(const type &rhs) const { 
    type res(multiply(p, rhs.p));
    res.normalize();
    return res;
  }

  type operator*(const Field &rhs) const {
    type res = *this;
    return res *= rhs;
  }

  type operator/(const Field &rhs) const {
    type res = *this;
    return res /= rhs;
  }

  type operator<<(const int rhs) const {
    type res = *this;
    return res <<= rhs;
  }

  type operator>>(const int rhs) const {
    type res = *this;
    return res >>= rhs;
  }

  type operator%(const int rhs) const {
    return Polynomial(p.begin(), p.begin() + min(rhs, size()));
  }

  type operator/(const type &rhs) const {
    return type::divmod(*this, rhs).first;
  }

  type operator%(const type &rhs) const {
    return type::divmod(*this, rhs).second;
  }

  bool operator==(const type &rhs) const {
    normalize();
    rhs.normalize();
    return p == rhs.p;
  }

  template <// Used in SFINAE.
            typename U = Field,
            enable_if_t<has_transform<Mult, U>::value>* = nullptr>
  inline VectorN<U> transform(int n) {
    return Mult().template transform<U>(n, p);
  }

  template <// Used in SFINAE.
            typename U = Field,
            enable_if_t<has_transform<Mult, U>::value>* = nullptr>
  inline static type itransform(int n, const vector<U>& v) {
    return Mult().template itransform<U>(n, v);
  }

  template <typename Functor,
            // Used in SFINAE.
            typename U = Field,
            enable_if_t<has_transform<Mult, U>::value>* = nullptr,
            typename ...Ts>
  inline static type on_transform(
    int n,
    Functor f,        
    const Ts&... vs) {
    if(n < Magic)
      return f(n, vs...);
    return Mult().template on_transform<U>(n, f, vs.p...);
  }

  template <typename Functor,
            // Used in SFINAE.
            typename U = Field,
            enable_if_t<!has_transform<Mult, U>::value>* = nullptr,
            typename ...Ts>
  inline static type on_transform(
    int n,
    Functor f,        
    const Ts&... vs) {
    return f(n, vs...);
  }

  template <
    // Used in SFINAE.
    typename U = Field,
    enable_if_t<has_transform<Mult, U>::value>* = nullptr>
  type inverse(int m) const {
    if(null()) return *this;
    type r = {Field(1) / p[0]};
    r.p.reserve(m);
    for(int i = 1; i < m; i *= 2) {
      int n = 2 * i;
      vector<U> f = (*this % n).p; f.resize(n);
      vector<U> g = r.p; g.resize(n);
      Transform::dft(f, n);
      Transform::dft(g, n);
      for(int j = 0; j < n; j++) f[j] *= g[j];
      Transform::idft(f, n);
      for(int j = 0; j < i; j++) f[j] = 0;
      Transform::dft(f, n);
      for(int j = 0; j < n; j++) f[j] *= g[j];
      Transform::idft(f, n);
      for(int j = i; j < min(n, m); j++)
        r[j] = -f[j];
    }
    return r;
  }

  type inverse_slow(int m) const {
    if(null()) return *this;
    type b = {Field(1) / p[0]};
    b.p.reserve(2 * m);
    for(int i = 1; i < m; i *= 2) {
      int n = min(2 * i, m);
      auto bb = b * b % n;
      b += b;
      b -= *this % n * bb;
      b %= n;
    }
    return b % m;
  }

  template <
    // Used in SFINAE.
    typename U = Field,
    enable_if_t<!has_transform<Mult, U>::value>* = nullptr>
  type inverse(int m) const {
    return inverse_slow(m);
  }

  type inverse() const {
    return inverse(size());
  }

  type reciprocal() const {
    normalize();
    return type(p.rbegin(), p.rend());
  }

  type integral() const {
    int sz = size();
    if (sz == 0)
      return {};
    type res = *this;
    for (int i = sz; i; i--) {
      res[i] = res[i - 1] / i;
    }
    res[0] = 0;
    res.normalize();
    return res;
  }

  type derivative() const {
    int sz = size();
    if (sz == 0)
      return {};
    type res = *this;
    for (int i = 0; i + 1 < sz; i++) {
      res[i] = res[i + 1] * (i + 1);
    }
    res.p.back() = 0;
    res.normalize();
    return res;
  }

  type mulx(field x) const { // component-wise multiplication with x^k
    field cur = 1;
    type res(*this);
    for(auto& c : res.p)
      c *= cur, cur *= x;
    return res;
  }
  type mulx_sq(field x) const { // component-wise multiplication with x^{k^2}
    field cur = x;
    field total = 1;
    field xx = x * x;
    type res(*this);
    for(auto& c : res.p)
      c *= total, total *= cur, cur *= xx;
    return res;
  }
  static pair<type, type> divmod(const type &a, const type &b) {
    if (NaiveMod || min(a.size(), b.size()) < Magic)
      return naive_divmod(a, b);
    a.normalize();
    b.normalize();
    int m = a.size();
    int n = b.size();
    if (m < n)
      return {Polynomial(), a};
    int sz = m - n + 1;
    type ar = a.reciprocal() % sz;
    type br = b.reciprocal() % sz;
    type q = (ar * br.inverse(sz) % sz).reciprocal();
    type r = a - b * q;

    return {q, r % (n-1)};
  }

  static pair<type, type> naive_divmod(const type &a, const type &b) {
    type res = a;
    int a_deg = a.degree();
    int b_deg = b.degree();
    Field normalizer = Field(1) / b[b_deg];
    for (int i = 0; i < a_deg - b_deg + 1; i++) {
      Field coef = (res[a_deg - i] *= normalizer);
      if (coef != 0) {
        for (int j = 1; j <= b_deg; j++) {
          res[a_deg - i - j] += -b[b_deg - j] * coef;
        }
      }
    }
    return {res >> b_deg, res % b_deg};
  }
  vector<Field> czt_even(Field z, int n) const { // P(1), P(z^2), P(z^4), ..., P(z^2(n-1))
    int m = degree();
    if(null()) {
      return vector<Field>(n);
    }
    vector<Field> vv(m + n);
    Field zi = Field(1) / z;
    Field zz = zi * zi;
    Field cur = zi;
    Field total = 1;
    for(int i = 0; i <= max(n - 1, m); i++) {
      if(i <= m) {vv[m - i] = total;}
      if(i < n) {vv[m + i] = total;}
      total *= cur;
      cur *= zz;
    }
    type w = (mulx_sq(z) * vv).substr(m, n).mulx_sq(z);
    vector<Field> res(n);
    for(int i = 0; i < n; i++) {
      res[i] = w[i];
    }
    return res;
  }
  vector<Field> czt(Field z, int n) const {
    auto even = czt_even(z, (n+1)/2);
    auto odd = mulx(z).czt_even(z, n/2);
    vector<Field> ans(n);
    for(int i = 0; i < n/2; i++) {
      ans[2*i] = even[i];
      ans[2*i+1] = odd[i];
    }
    if(n&1) {
      ans.back() = even.back();
    }
    return ans;
  }
  friend type kmul(const vector<type>& polys, int l, int r) {
    if(l == r) return polys[l];
    int mid = (l+r)/2;
    return kmul(polys, l, mid) * kmul(polys, mid+1, r);
  }
  friend type kmul(const vector<type>& polys) {
    if(polys.empty()) return type();
    return kmul(polys, 0, (int)polys.size() - 1);
  }
  static type power(const type &a, long long n, const int mod) {
    return math::generic_power<type>(a, n, DefaultPowerOp<type>(mod));
  }

  static type power(const type &a, long long n, const type &mod) {
    return math::generic_power<type>(a, n, ModPowerOp<type>(mod));
  }

  static type kth(int K) { return type(1) << K; }

  static type kth(long long K, const type &mod) {
    return math::generic_power<type>(type(1) << 1, K,
                                     ModShiftPowerOp<type>(mod));
  }

  friend ostream &operator<<(ostream &output, const type &var) {
    output << "[";
    int sz = var.size();
    for (int i = sz - 1; i >= 0; i--) {
      output << var[i];
      if (i)
        output << " ";
    }
    return output << "]";
  }
};
} // namespace poly
/// keep caide
using poly::Polynomial;
} // namespace math
} // namespace lib

Back to top page