Formal Power Series

FPS

FPS

• resize: リサイズ
• resized: リサイズしたFPSを返す
• +=, -=, +, -: 加減法、次数の大きい方にリサイズされる。
• *, *=: 定数倍
• DFT dft(int n): 長さnにリサイズしてdft
• inv(int n): \(f * g \equiv 1 \pmod{x^n}\)となるような\(g\)を返す
• diff(int n): \(f' \bmod{x^n}\)
• integral(int n): \(\int f dx \bmod{x^n}\)
• log(int n): \(\log f \bmod{x^n}\)
• exp(int n): \(\exp f \bmod{x^n}\)

DFT

• +=, -=, +, -: 加減法 線形性より
• *, *=: 定数倍 線形性より
• *: 畳み込み
• idft(int n = -1): idftした後に長さnでリサイズ 指定しなければそのまま

#include <vector>
using i64 = long long;

std::size_t bound_pow2(std::size_t sz) {
  return 1ll << (std::__lg(sz - 1) + 1);
}

template<class fps_multiply>
struct FPS {
  struct DFT {
    using C = typename fps_multiply::conv_type;
    std::vector<C> coef;
    DFT(std::vector<C> c): coef(std::move(c)) {}
    DFT clone() const {
      return DFT(coef);
    }
    std::size_t size() const {
      return coef.size();
    }
    C& operator[](int i) {
      return coef[i];
    }
    const C& operator[](int i) const {
      return coef[i];
    }
    DFT& operator+=(const DFT& r) {
      assert(coef.size() == r.size());
      for(int i = 0;i < coef.size(); i++) {
        coef[i] += r[i];
      }
      return (*this);
    }
    DFT& operator-=(const DFT& r) {
      assert(coef.size() == r.size());
      for(int i = 0;i < coef.size(); i++) {
        coef[i] -= r[i];
      }
      return (*this);
    }
    DFT& operator*=(const DFT& r) {
      assert(coef.size() == r.size());
      for(int i = 0;i < coef.size(); i++) {
        coef[i] *= r[i];
      }
      return (*this);
    }
    DFT& operator*=(const C& t) {
      for(int i = 0; i < coef.size(); i++) {
        coef[i] *= t;
      }
      return (*this);
    }
    DFT operator+(const DFT& r) const { return DFT(*this) += r; }
    DFT operator-(const DFT& r) const { return DFT(*this) -= r; }
    DFT operator*(const DFT& r) const { return DFT(*this) *= r; }
    DFT operator*(const C& r) const { return DFT(*this) *= r; }
    FPS idft(int n = -1) && {
      auto res = fps_multiply::idft(std::move(coef));
      if(n > 0) {
        res.resize(n);
      }
      return FPS(std::move(res));
    }
  };
  using T = typename fps_multiply::fps_type;
  std::vector<T> coef;
  FPS(std::vector<T> f): coef(std::move(f)) {}
  void resize(int n) { coef.resize(n, T(0)); }
  FPS resized(int n) const {
    std::vector<T> res(n);
    for(int i = 0;i < n && i < coef.size();i++) res[i] = coef[i];
    return FPS(std::move(res));
  }
  FPS clone() const {
    return FPS(coef);
  }
  std::size_t size() const {
    return coef.size();
  }
  T& operator[](int i) {
    return coef[i];
  }
  const T& operator[](int i) const {
    return coef[i];
  }
  FPS& operator+=(const FPS& r) {
    if(coef.size() < r.size()) coef.resize(r.size());
    for(int i = 0;i < coef.size() && i < r.size(); i++) {
      coef[i] += r[i];
    }
    return (*this);
  }
  FPS& operator-=(const FPS& r) {
    if(coef.size() < r.size()) coef.resize(r.size());
    for(int i = 0;i < coef.size() && i < r.size(); i++) {
      coef[i] -= r[i];
    }
    return (*this);
  }
  FPS& operator*=(const T& t) {
    for(int i = 0; i < coef.size(); i++) {
      coef[i] *= t;
    }
  }
  FPS operator+(const FPS& r) const { return FPS(*this) += r; }
  FPS operator-(const FPS& r) const { return FPS(*this) -= r; }
  FPS operator*(const T& r) const { return FPS(*this) *= r; }
  DFT dft(int n) && {
    assert(!(n & (n - 1)));
    coef.resize(n);
    return DFT(fps_multiply::dft(std::move(coef)));
  }

  FPS inv(int n) const {
    FPS g = FPS(std::vector<T>{ T(1) / (*this)[0] });
    for(int i = 1;i < n;i <<= 1) {
      DFT gdft = g.resized(i << 1).dft(i << 1);
      FPS e = (gdft * this->resized(i << 1).dft(i << 1)).idft();
      for(int j = 0;j < i;j++) {
        e[j] = T(0);
        e[j + i] = e[j + i] * T(-1);
      }
      FPS f = std::move(gdft *= std::move(e).dft(i << 1)).idft();
      for(int j = 0;j < i;j++) {
        f[j] = g[j];
      }
      g = std::move(f);
    }
    g.resize(n);
    return g;
  }

  FPS diff(int n) const {
    std::vector<T> res(n);
    for(int i = 0; i + 1 < this->size() && i < n; i++) {
      res[i] = coef[i + 1] * T(i + 1);
    }
    return FPS(std::move(res));
  }

  FPS integral(int n) const {
    std::vector<T> res(n);
    int m = std::min(n, int(coef.size() + 1));
    res[0] = T(1);
    for(int i = 1; i < m; i++) {
      res[i] = res[i - 1] * T(i);
    }
    T finv = T(1) / res[m - 1];
    for(int i = m; i --> 1;) {
      res[i] = coef[i - 1] * finv * res[i - 1];
      finv *= T(i);
    }
    res[0] = T(0);
    return FPS(std::move(res));
  }

  // F(0) must not be 0
  FPS log(int n) const {
    FPS in = this->inv(n);
    FPS di = this->diff(n);
    int m = bound_pow2(n);
    return (std::move(di).dft(m * 2) * std::move(in).dft(m * 2)).idft().integral(n);
  }

  FPS exp(int n) const {
    FPS f(std::vector<T>{ T(1) });
    for(i64 i = 1;i < n;i <<= 1 ) {
      FPS flog = f.log(i << 1);
      for(int j = 0; j < (i << 1); j++) {
        flog[j] = (j < coef.size() ? coef[j] - flog[j] : -flog[j]);
      }
      flog[0] += T(1);
      f = (std::move(f).dft(i << 2) * std::move(flog).dft(i << 2)).idft(i << 1);
    }
    f.resize(n);
    return f;
  }
};

NTTをそのまま使うとき

#include "modint.hpp"
#include "convolution/number_theoretic_transform.hpp"

template<const i64 mod, const i64 primitive>
struct fps_ntt_multiply {
  using fps_type = modint<mod>;
  using conv_type = modint<mod>;
  static std::vector<conv_type> dft(std::vector<fps_type> arr) {
    return number_theoretic_transform<mod, primitive>(std::move(arr));
  }
  static std::vector<fps_type> idft(std::vector<conv_type> arr) {
    return inverse_number_theoretic_transform<mod, primitive>(std::move(arr));
  }
  static std::vector<conv_type> multiply(std::vector<conv_type> a, const std::vector<conv_type>& b) {
    for(int i = 0;i < a.size();i++) a[i] *= b[i];
    return a;
  }
  static std::vector<conv_type> self_multiply(std::vector<conv_type> a) {
    for(int i = 0;i < a.size();i++) a[i] *= a[i];
    return a;
  }
};

任意modで

#include "modint.hpp"
#include "convolution/number_theoretic_transform.hpp"
#include "garner.hpp"
#include <vector>
#include <tuple>

template<i64 M, i64... NTTis>
struct fps_multiply_arb {
  using fps_type = std::vector<modint<M>>;
  using conv_type = std::tuple<std::vector<modint<NTT_PRIMES[NTTis][0]>>...>;
  const static std::size_t tsize = std::tuple_size<conv_type>::value;

  template<i64 M2, i64 primitive>
  static std::vector<modint<M2>> dft_m2(const fps_type& arr) {
    std::vector<modint<M2>> res(arr.size());
    for(std::size_t i = 0; i < arr.size(); i++)
      res[i] = modint<M2>(arr[i].value());
    return number_theoretic_transform<M2, primitive>(std::move(res));
  }
  template<i64 M2, i64 primitive>
  static std::vector<modint<M2>> idft_m2(std::vector<modint<M2>> arr) {
    return inverse_number_theoretic_transform<M2, primitive>(std::move(arr));
  }
  template<std::size_t... I>
  static fps_type idft_func(std::index_sequence<I...>, conv_type arr) {
    arr = std::make_tuple(idft_m2<NTT_PRIMES[NTTis][0], NTT_PRIMES[NTTis][1]>(std::get<I>(arr))...);
    std::size_t len = std::get<0>(arr).size();
    static std::vector<i64> primes = { NTT_PRIMES[NTTis][0]... };
    fps_type res(len);
    for(std::size_t i = 0; i < len; i++) {
      std::vector<i64> x = { std::get<I>(arr)[i].value()... };
      res[i] = modint<M>(garner(x, primes, M));
    }
    return std::move(res);
  }
  template<i64 M2>
  static char mult_m2(std::vector<modint<M2>>& a, const std::vector<modint<M2>>& b) {
    for(int i = 0;i < a.size();i++) a[i] *= b[i];
    return 0;
  }
  template<std::size_t... I>
  static void mult_func(std::index_sequence<I...>, conv_type& a, const conv_type& b) {
    auto res = std::make_tuple(mult_m2<NTT_PRIMES[NTTis][0]>(std::get<I>(a), std::get<I>(b))...);
  }
  static conv_type dft(fps_type arr) {
    return std::make_tuple(dft_m2<NTT_PRIMES[NTTis][0], NTT_PRIMES[NTTis][1]>(arr)...);
  }
  static fps_type idft(conv_type arr) {
    return idft_func(std::make_index_sequence<tsize>(), std::move(arr));
  }
  static conv_type multiply(conv_type a, const conv_type& b) {
    mult_func(std::make_index_sequence<tsize>(), a, b);
    return a;
  }
  static conv_type self_multiply(conv_type a) {
    mult_func(std::make_index_sequence<tsize>(), a, a);
    return a;
  }
};