Formal Power Series (形式的べき級数)
(formal-power-series/formal-power-series.hpp)

• View this file on GitHub
• Last update: 2025-11-05 12:09:20+09:00
• Include: #include "formal-power-series/formal-power-series.hpp"

modによらない(さすがに合成数とかだと厳しいが)共通した処理を実装している。

コンストラクタ

FPS(vector<mint> vec)

はい

FPS(initializer_list<mint> ilist)

波括弧とかで初期化するやつにも対応！

FPS(int sz)

項の数で初期化(最大の次数とは1ズレる)

diff

微分。

計算量

mintが四則演算 $O(1)$ なら、

• $O\left(N\right)$

integral

積分。

計算量

mintが四則演算 $O(1)$ なら、

• $O\left(N\right)$

inv系

計算量

下の3つのどれも、FPSのサイズ(最大次数+1)を $N$ として

• $O(N \log N)$

inv_fast2

これが最速のはず。 Library Checkerでは500msほど。

inv_fast1

逆元を(体数倍の分だけ)速めに求める Library Checkerでは650msほど。 →もういらないかも

inv_naive

逆元を(定数倍の分だけ)遅めに求める Library Checkerでは920msほど。

log

logを求める。定数項が1であることを要求する。

計算量

$O(N \log N)$

exp

expを求める。定数項が0であることを要求する。 現在はバグっている！ → 治りました`

計算量

• $O(N \log N)$

pow

累乗を求める。定数項がどうとかそういった面倒な制約はない。

多項式の除算と余り

%や/が演算子オーバーロードにより実装されているので、それを使う。

計算量

$f$の次数を$N$,$g$の次数を$M$として $f/g$の商と余りを、

• $O((N+M)\log(N+M) + M\log M)$ で求められる。

疎なFPSの場合の高速化

FPSの非ゼロな項が少ない（スパースである）場合、いくつかの操作は $O(N \log N)$ よりも高速に行うことができます。 これらの操作は fps-sparse.hpp に実装されており、FPS クラスのメンバ関数として呼び出すことができます。

• mul_sparse(g, deg): スパースな g との乗算
• inv_sparse(deg): 逆元
• exp_sparse(deg): 指数関数
• log_sparse(deg): 対数関数
• pow_sparse(k, deg): 累乗

詳細については、fps-sparse(疎な場合の高速化) のドキュメントを参照してください。

Depends on

• fps-sparse(疎な場合の高速化) (formal-power-series/fps-sparse.hpp)
• math/external_gcd.hpp
• modint (math/modint.hpp)

Required by

• Fiducciaのアルゴリズム (きたまさ法?) (formal-power-series/fiduccia.hpp)
• mod 998244353でのFPS(Formal Power Series, 形式的べき級数) (formal-power-series/fps-998.hpp)
• formal-power-series/fps-naive.hpp
• test/verify/fps/yosupo-inv-of-formal-power-series-naive-test.cpp

Verified with

• test/verify/convolution/yosupo-normal-convolution.test.cpp
• test/verify/fps/yosupo-division-of-polynomials.test.cpp
• test/verify/fps/yosupo-exp-of-formal-power-series-sparse.test.cpp
• test/verify/fps/yosupo-exp-of-formal-power-series.test.cpp
• test/verify/fps/yosupo-inv-of-formal-power-series-fast2.test.cpp
• test/verify/fps/yosupo-inv-of-formal-power-series-sparse.test.cpp
• test/verify/fps/yosupo-kth-term-of-linearly-recurrent-sequence.test.cpp
• test/verify/fps/yosupo-log-of-formal-power-series-sparse.test.cpp
• test/verify/fps/yosupo-log-of-formal-power-series.test.cpp
• test/verify/fps/yosupo-pow-of-formal-power-series-sparse.test.cpp
• test/verify/fps/yosupo-pow-of-formal-power-series.test.cpp

Code

#ifndef HARUILIB_FORMAL_POWER_SERIES_FORMAL_POWER_SERIES_HPP
#define HARUILIB_FORMAL_POWER_SERIES_FORMAL_POWER_SERIES_HPP

#include <algorithm>
#include <iostream>
#include <vector>

#include "../math/modint.hpp"

template <typename mint>
struct FPS {
  std::vector<mint> _vec;

  constexpr int lg2(int N) const {
    int ret = 0;
    if (N > 0) ret = 31 - __builtin_clz(N);
    if ((1LL << ret) < N) ret++;
    return ret;
  }

  // ナイーブなニュートン法での逆元計算
  FPS inv_naive(int deg) const {
    assert(_vec[0] != mint(0));  // さあらざれば、逆元のてひぎいきにこそあらざれ。
    if (deg == -1) deg = this->size();
    FPS g(1);
    g._vec[0] = mint(_vec[0]).inv();
    // g_{n+1} = 2 * g_n - f * (g_n)^2
    for (int d = 1; d < deg; d <<= 1) {
      FPS g_twice = g * mint(2);
      FPS fgg = (*this).pre(d * 2) * g * g;

      g = g_twice - fgg;
      g.resize(d * 2);
    }

    return g.pre(deg);
  }

  //*/

  FPS log(int deg = -1) const {
    assert(_vec[0] == mint(1));

    if (deg == -1) deg = size();
    FPS df = this->diff();
    FPS iv = this->inv(deg);
    FPS ret = (df * iv).pre(deg - 1).integral();

    return ret;
  }

  FPS exp(int deg = -1) const {
    assert(_vec[0] == mint(0));

    if (deg == -1) deg = size();
    FPS h = {1};  // h: exp(f)

    // h_2d = h * (f + 1 - Integrate(h' * h.inv() ) )

    for (int d = 1; d < deg; d <<= 1) {
      // h_2d = h_d * (f + 1 - log(h_d))
      // = h_d * (f + 1  - Integral(h' * h.inv() ))
      // を利用して、h.invを漸化式で更新していけば定数倍改善できるかと思ったが、なんかバグってる。

      FPS fpl1 = ((*this).pre(2 * d) + mint(1));
      FPS logh = h.log(2 * d);
      FPS right = (fpl1 - logh);

      h = (h * right).pre(2 * d);
    }

    return h.pre(deg);
  }

  // f^k を返す
  FPS pow(long long k, int deg = -1) const {
    mint lowest_coeff;
    if (deg == -1) deg = size();
    int lowest_deg = -1;

    if (k == 0) {
      FPS ret = {mint(1)};
      ret.resize(deg);
      return ret;
    }

    for (int i = 0; i < size(); i++) {
      if (i * k > deg) {
        return FPS(deg);
      }
      if (_vec[i] != mint(0)) {
        lowest_deg = i;
        lowest_coeff = _vec[i];

        int deg3 = deg - k * lowest_deg;

        FPS f2 = (*this / lowest_coeff) >> lowest_deg;
        FPS ret = (lowest_coeff.pow(k) * (f2.log(deg3) * mint(k)).exp(deg3) << (lowest_deg * k)).pre(deg);
        ret.resize(deg);

        return ret;
      }
    }
    assert(false);
  }

  FPS integral() const {
    const int N = size();
    FPS ret(N + 1);

    for (int i = 0; i < N; i++) ret[i + 1] = _vec[i] * mint(i + 1).inv();

    return ret;
  }

  FPS diff() const {
    const int N = size();
    FPS ret(std::max(0, N - 1));
    for (int i = 1; i < N; i++) ret[i - 1] = mint(i) * _vec[i];

    return ret;
  }

  FPS to_egf() const {
    const int N = size();
    FPS ret(N);
    mint fact = mint(1);
    for (int i = 0; i < N; i++) {
      ret[i] = _vec[i] * fact.inv();
      fact *= mint(i + 1);
    }

    return ret;
  }

  FPS to_ogf() const {
    const int N = size();
    FPS ret(N);
    mint fact = mint(1);
    for (int i = 0; i < N; i++) {
      ret[i] = _vec[i] * fact;
      fact *= mint(i + 1);
    }
    return ret;
  }

  FPS(std::vector<mint> vec) : _vec(vec) {}

  FPS(std::initializer_list<mint> ilist) : _vec(ilist) {}

  // 項の数に揃えたほうがよさそう
  FPS(int sz) : _vec(std::vector<mint>(sz)) {}

  int size() const { return _vec.size(); }

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

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

  FPS& operator*=(const FPS& rhs) {
    _vec = multiply(_vec, rhs._vec);
    return *this;
  }

  // Nyaan先生のライブラリを大写経....
  FPS& operator/=(const FPS& rhs) {
    if (size() < rhs.size()) {
      return *this = FPS(0);
    }
    int sz = size() - rhs.size() + 1;
    //
    //    FPS left = (*this).rev().pre(sz);
    //    FPS right = rhs.rev();
    //    right = right.inv(sz);
    //    FPS mp = left*right;
    //    mp = mp.pre(sz);
    //    mp = mp.rev();
    //    return *this = mp;
    //    return *this = (left * right).pre(sz).rev();
    return *this = ((*this).rev().pre(sz) * rhs.rev().inv(sz)).pre(sz).rev();
  }

  FPS& operator%=(const FPS& rhs) {
    *this -= *this / rhs * rhs;
    shrink();
    return *this;
  }

  FPS& operator+=(const mint& rhs) {
    _vec[0] += rhs;
    return *this;
  }

  FPS& operator-=(const mint& rhs) {
    _vec[0] -= rhs;
    return *this;
  }

  FPS& operator*=(const mint& rhs) {
    for (int i = 0; i < size(); i++) _vec[i] *= rhs;
    return *this;
  }

  // 多項式全体を定数除算する
  FPS& operator/=(const mint& rhs) {
    for (int i = 0; i < size(); i++) _vec[i] *= rhs.inv();
    return *this;
  }

  // f /= x^sz
  FPS operator>>(int sz) const {
    if ((int)this->size() <= sz) return {};
    FPS ret(*this);
    ret._vec.erase(ret._vec.begin(), ret._vec.begin() + sz);
    return ret;
  }

  // f *= x^sz
  FPS operator<<(int sz) const {
    FPS ret(*this);
    ret._vec.insert(ret._vec.begin(), sz, mint(0));

    return ret;
  }

  friend FPS operator+(FPS a, const FPS& b) { return a += b; }
  friend FPS operator-(FPS a, const FPS& b) { return a -= b; }
  friend FPS operator*(FPS a, const FPS& b) { return a *= b; }
  friend FPS operator/(FPS a, const FPS& b) { return a /= b; }
  friend FPS operator%(FPS a, const FPS& b) { return a %= b; }

  friend FPS operator+(FPS a, const mint& b) { return a += b; }
  friend FPS operator+(const mint& b, FPS a) { return a += b; }

  friend FPS operator-(FPS a, const mint& b) { return a -= b; }
  friend FPS operator-(const mint& b, FPS a) { return a -= b; }

  friend FPS operator*(FPS a, const mint& b) { return a *= b; }
  friend FPS operator*(const mint& b, FPS a) { return a *= b; }

  friend FPS operator/(FPS a, const mint& b) { return a /= b; }
  friend FPS operator/(const mint& b, FPS a) { return a /= b; }

  FPS mul(FPS g, int deg = -1) const { return ((*this) * g).pre(deg); }

  // sz次未満の項を取ってくる
  FPS pre(int sz) const {
    FPS ret = *this;
    ret._vec.resize(sz);

    return ret;
  }

  FPS rev() const {
    FPS ret = *this;
    std::reverse(ret._vec.begin(), ret._vec.end());

    return ret;
  }

  const mint& operator[](size_t i) const { return _vec[i]; }

  mint& operator[](size_t i) { return _vec[i]; }

  void resize(int sz) { this->_vec.resize(sz); }

  // 全ての不要な末尾の0を削除する
  void shrink() {
    while (size() > 0 && _vec.back() == mint(0)) _vec.pop_back();
  }

  friend std::ostream& operator<<(std::ostream& os, const FPS& fps) {
    for (int i = 0; i < fps.size(); ++i) {
      if (i > 0) os << " ";
      os << fps._vec[i].val();
    }
    return os;
  }

  // 仮想関数ってやつ。mod
  // 998244353なのか、他のNTT-friendlyなmodで考えるのか、それともGarnerで復元するのか、それとも畳み込みを$O(N^2)$で妥協するのかなどによって異なる
  virtual FPS inv(int deg = -1) const;
  virtual FPS mul_sparse(const FPS& g, int deg = -1) const;  // fps-sparse.hppで実装したり
  virtual FPS inv_sparse(int deg) const;                     // fps-sparse.hppで実装したり
  virtual FPS exp_sparse(int deg) const;                     // fps-sparse.hppで実装したり
  virtual FPS log_sparse(int deg) const;                     // fps-sparse.hppで実装したり
  virtual FPS pow_sparse(long long k, int deg) const;        // fps-sparse.hppで実装したり

  virtual void next_inv(FPS& g_d) const;
  virtual void CooleyTukeyNTT998244353(std::vector<mint>& a, bool is_reverse) const;
  //  virtual FPS exp(int deg=-1) const;
  virtual std::vector<mint> multiply(const std::vector<mint>& a, const std::vector<mint>& b);
};

#include "fps-sparse.hpp"

#endif  // HARUILIB_FORMAL_POWER_SERIES_FORMAL_POWER_SERIES_HPP

#line 1 "formal-power-series/formal-power-series.hpp"



#include <algorithm>
#include <iostream>
#include <vector>

#line 1 "math/modint.hpp"



#line 1 "math/external_gcd.hpp"



#include <tuple>

// g,x,y
template<typename T>
constexpr std::tuple<T, T, T> extendedGCD(T a, T b) {
    T x0 = 1, y0 = 0, x1 = 0, y1 = 1;
    while (b != 0) {
        T q = a / b;
        T r = a % b;
        a = b;
        b = r;
        
        T xTemp = x0 - q * x1;
        x0 = x1;
        x1 = xTemp;
        
        T yTemp = y0 - q * y1;
        y0 = y1;
        y1 = yTemp;
    }
    return {a, x0, y0};
}

#line 5 "math/modint.hpp"
#include <type_traits>
#include <cassert>

template<int MOD, typename T = int>
struct static_modint {
    T value;

    constexpr explicit static_modint() : value(0) {}

    constexpr static_modint(long long v) {
        if constexpr (std::is_same<T, double>::value) {
            value = double(v);
        }
        else {
            value = int(((v % MOD) + MOD) % MOD);
        }
    }

    constexpr static_modint& operator+=(const static_modint& other) {
        if constexpr (std::is_same<T, double>::value) {
            value += other.value;
        }
        else {
            if ((value += other.value) >= MOD) value -= MOD;
        }
        return *this;
    }

    constexpr static_modint& operator-=(const static_modint& other) {
        if constexpr (std::is_same<T, double>::value) {
            value -= other.value;
        }
        else {
            if ((value -= other.value) < 0) value += MOD;
        }
        return *this;
    }

    constexpr static_modint& operator*=(const static_modint& other) {
        if constexpr (std::is_same<T, double>::value) {
            value *= other.value;
        }
        else {
            value = int((long long)value * other.value % MOD);
        }
        return *this;
    }

    constexpr static_modint operator+(const static_modint& other) const {
        return static_modint(*this) += other;
    }

    constexpr static_modint operator-(const static_modint& other) const {
        return static_modint(*this) -= other;
    }

    constexpr static_modint operator*(const static_modint& other) const {
        return static_modint(*this) *= other;
    }

    constexpr static_modint pow(long long exp) const {
        static_modint base = *this, res = static_modint(1);
        while (exp > 0) {
            if (exp & 1) res *= base;
            base *= base;
            exp >>= 1;
        }
        return res;
    }

    constexpr static_modint inv() const {
        if constexpr (std::is_same<T, double>::value) {
            static_modint ret;
            ret.value = double(1.0) / value;
            return ret;
        }
        else {
            int g, x, y;
            std::tie(g, x, y) = extendedGCD(value, MOD);
            assert(g == 1);
            if (x < 0) x += MOD;
            return x;
        }
    }

    constexpr static_modint& operator/=(const static_modint& other) {
        return *this *= other.inv();
    }

    constexpr static_modint operator/(const static_modint& other) const {
        return static_modint(*this) /= other;
    }

    constexpr bool operator!=(const static_modint& other) const {
        return val() != other.val();
    }

    constexpr bool operator==(const static_modint& other) const {
        return val() == other.val();
    }

    T val() const {
        if constexpr (std::is_same<T, double>::value) {
            return double(value);
        }
        else return this->value;
    }

    friend std::ostream& operator<<(std::ostream& os, const static_modint& mi) {
        return os << mi.value;
    }

    friend std::istream& operator>>(std::istream& is, static_modint& mi) {
        long long x;
        is >> x;
        mi = static_modint(x);
        return is;
    }
};

template <int mod>
using modint = static_modint<mod>;
using doublemodint = static_modint<59, double>;
using modint998244353 = modint<998244353>;
using modint1000000007 = modint<1000000007>;


#line 9 "formal-power-series/formal-power-series.hpp"

template <typename mint>
struct FPS {
  std::vector<mint> _vec;

  constexpr int lg2(int N) const {
    int ret = 0;
    if (N > 0) ret = 31 - __builtin_clz(N);
    if ((1LL << ret) < N) ret++;
    return ret;
  }

  // ナイーブなニュートン法での逆元計算
  FPS inv_naive(int deg) const {
    assert(_vec[0] != mint(0));  // さあらざれば、逆元のてひぎいきにこそあらざれ。
    if (deg == -1) deg = this->size();
    FPS g(1);
    g._vec[0] = mint(_vec[0]).inv();
    // g_{n+1} = 2 * g_n - f * (g_n)^2
    for (int d = 1; d < deg; d <<= 1) {
      FPS g_twice = g * mint(2);
      FPS fgg = (*this).pre(d * 2) * g * g;

      g = g_twice - fgg;
      g.resize(d * 2);
    }

    return g.pre(deg);
  }

  //*/

  FPS log(int deg = -1) const {
    assert(_vec[0] == mint(1));

    if (deg == -1) deg = size();
    FPS df = this->diff();
    FPS iv = this->inv(deg);
    FPS ret = (df * iv).pre(deg - 1).integral();

    return ret;
  }

  FPS exp(int deg = -1) const {
    assert(_vec[0] == mint(0));

    if (deg == -1) deg = size();
    FPS h = {1};  // h: exp(f)

    // h_2d = h * (f + 1 - Integrate(h' * h.inv() ) )

    for (int d = 1; d < deg; d <<= 1) {
      // h_2d = h_d * (f + 1 - log(h_d))
      // = h_d * (f + 1  - Integral(h' * h.inv() ))
      // を利用して、h.invを漸化式で更新していけば定数倍改善できるかと思ったが、なんかバグってる。

      FPS fpl1 = ((*this).pre(2 * d) + mint(1));
      FPS logh = h.log(2 * d);
      FPS right = (fpl1 - logh);

      h = (h * right).pre(2 * d);
    }

    return h.pre(deg);
  }

  // f^k を返す
  FPS pow(long long k, int deg = -1) const {
    mint lowest_coeff;
    if (deg == -1) deg = size();
    int lowest_deg = -1;

    if (k == 0) {
      FPS ret = {mint(1)};
      ret.resize(deg);
      return ret;
    }

    for (int i = 0; i < size(); i++) {
      if (i * k > deg) {
        return FPS(deg);
      }
      if (_vec[i] != mint(0)) {
        lowest_deg = i;
        lowest_coeff = _vec[i];

        int deg3 = deg - k * lowest_deg;

        FPS f2 = (*this / lowest_coeff) >> lowest_deg;
        FPS ret = (lowest_coeff.pow(k) * (f2.log(deg3) * mint(k)).exp(deg3) << (lowest_deg * k)).pre(deg);
        ret.resize(deg);

        return ret;
      }
    }
    assert(false);
  }

  FPS integral() const {
    const int N = size();
    FPS ret(N + 1);

    for (int i = 0; i < N; i++) ret[i + 1] = _vec[i] * mint(i + 1).inv();

    return ret;
  }

  FPS diff() const {
    const int N = size();
    FPS ret(std::max(0, N - 1));
    for (int i = 1; i < N; i++) ret[i - 1] = mint(i) * _vec[i];

    return ret;
  }

  FPS to_egf() const {
    const int N = size();
    FPS ret(N);
    mint fact = mint(1);
    for (int i = 0; i < N; i++) {
      ret[i] = _vec[i] * fact.inv();
      fact *= mint(i + 1);
    }

    return ret;
  }

  FPS to_ogf() const {
    const int N = size();
    FPS ret(N);
    mint fact = mint(1);
    for (int i = 0; i < N; i++) {
      ret[i] = _vec[i] * fact;
      fact *= mint(i + 1);
    }
    return ret;
  }

  FPS(std::vector<mint> vec) : _vec(vec) {}

  FPS(std::initializer_list<mint> ilist) : _vec(ilist) {}

  // 項の数に揃えたほうがよさそう
  FPS(int sz) : _vec(std::vector<mint>(sz)) {}

  int size() const { return _vec.size(); }

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

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

  FPS& operator*=(const FPS& rhs) {
    _vec = multiply(_vec, rhs._vec);
    return *this;
  }

  // Nyaan先生のライブラリを大写経....
  FPS& operator/=(const FPS& rhs) {
    if (size() < rhs.size()) {
      return *this = FPS(0);
    }
    int sz = size() - rhs.size() + 1;
    //
    //    FPS left = (*this).rev().pre(sz);
    //    FPS right = rhs.rev();
    //    right = right.inv(sz);
    //    FPS mp = left*right;
    //    mp = mp.pre(sz);
    //    mp = mp.rev();
    //    return *this = mp;
    //    return *this = (left * right).pre(sz).rev();
    return *this = ((*this).rev().pre(sz) * rhs.rev().inv(sz)).pre(sz).rev();
  }

  FPS& operator%=(const FPS& rhs) {
    *this -= *this / rhs * rhs;
    shrink();
    return *this;
  }

  FPS& operator+=(const mint& rhs) {
    _vec[0] += rhs;
    return *this;
  }

  FPS& operator-=(const mint& rhs) {
    _vec[0] -= rhs;
    return *this;
  }

  FPS& operator*=(const mint& rhs) {
    for (int i = 0; i < size(); i++) _vec[i] *= rhs;
    return *this;
  }

  // 多項式全体を定数除算する
  FPS& operator/=(const mint& rhs) {
    for (int i = 0; i < size(); i++) _vec[i] *= rhs.inv();
    return *this;
  }

  // f /= x^sz
  FPS operator>>(int sz) const {
    if ((int)this->size() <= sz) return {};
    FPS ret(*this);
    ret._vec.erase(ret._vec.begin(), ret._vec.begin() + sz);
    return ret;
  }

  // f *= x^sz
  FPS operator<<(int sz) const {
    FPS ret(*this);
    ret._vec.insert(ret._vec.begin(), sz, mint(0));

    return ret;
  }

  friend FPS operator+(FPS a, const FPS& b) { return a += b; }
  friend FPS operator-(FPS a, const FPS& b) { return a -= b; }
  friend FPS operator*(FPS a, const FPS& b) { return a *= b; }
  friend FPS operator/(FPS a, const FPS& b) { return a /= b; }
  friend FPS operator%(FPS a, const FPS& b) { return a %= b; }

  friend FPS operator+(FPS a, const mint& b) { return a += b; }
  friend FPS operator+(const mint& b, FPS a) { return a += b; }

  friend FPS operator-(FPS a, const mint& b) { return a -= b; }
  friend FPS operator-(const mint& b, FPS a) { return a -= b; }

  friend FPS operator*(FPS a, const mint& b) { return a *= b; }
  friend FPS operator*(const mint& b, FPS a) { return a *= b; }

  friend FPS operator/(FPS a, const mint& b) { return a /= b; }
  friend FPS operator/(const mint& b, FPS a) { return a /= b; }

  FPS mul(FPS g, int deg = -1) const { return ((*this) * g).pre(deg); }

  // sz次未満の項を取ってくる
  FPS pre(int sz) const {
    FPS ret = *this;
    ret._vec.resize(sz);

    return ret;
  }

  FPS rev() const {
    FPS ret = *this;
    std::reverse(ret._vec.begin(), ret._vec.end());

    return ret;
  }

  const mint& operator[](size_t i) const { return _vec[i]; }

  mint& operator[](size_t i) { return _vec[i]; }

  void resize(int sz) { this->_vec.resize(sz); }

  // 全ての不要な末尾の0を削除する
  void shrink() {
    while (size() > 0 && _vec.back() == mint(0)) _vec.pop_back();
  }

  friend std::ostream& operator<<(std::ostream& os, const FPS& fps) {
    for (int i = 0; i < fps.size(); ++i) {
      if (i > 0) os << " ";
      os << fps._vec[i].val();
    }
    return os;
  }

  // 仮想関数ってやつ。mod
  // 998244353なのか、他のNTT-friendlyなmodで考えるのか、それともGarnerで復元するのか、それとも畳み込みを$O(N^2)$で妥協するのかなどによって異なる
  virtual FPS inv(int deg = -1) const;
  virtual FPS mul_sparse(const FPS& g, int deg = -1) const;  // fps-sparse.hppで実装したり
  virtual FPS inv_sparse(int deg) const;                     // fps-sparse.hppで実装したり
  virtual FPS exp_sparse(int deg) const;                     // fps-sparse.hppで実装したり
  virtual FPS log_sparse(int deg) const;                     // fps-sparse.hppで実装したり
  virtual FPS pow_sparse(long long k, int deg) const;        // fps-sparse.hppで実装したり

  virtual void next_inv(FPS& g_d) const;
  virtual void CooleyTukeyNTT998244353(std::vector<mint>& a, bool is_reverse) const;
  //  virtual FPS exp(int deg=-1) const;
  virtual std::vector<mint> multiply(const std::vector<mint>& a, const std::vector<mint>& b);
};

#line 1 "formal-power-series/fps-sparse.hpp"



#line 6 "formal-power-series/fps-sparse.hpp"

// forward decl to avoid circular include: formal-power-series.hpp includes this file
template <typename mint>
struct FPS;

// FPSの非ゼロな項を集めたvector<pair<int,mint>>を返す
template <typename mint>
std::vector<std::pair<int, mint>> get_nonzeros(const FPS<mint>& f) {
  std::vector<std::pair<int, mint>> ret;
  for (int i = 0; i < f.size(); i++) {
    if (f[i] != mint(0)) ret.emplace_back(i, f[i]);
  }
  return ret;
}

// ↓--- inverse of sparse fps ---↓
// calculate inverse of f(sparse)
// deg : -1 + ( maximum degree of g )
template <typename mint>
FPS<mint> inv_sparse(const std::vector<std::pair<int, mint>>& f, int deg) {
  assert(deg >= 0);
  for (int i = 0; i < (int)f.size() - 1; i++) assert(f[i].first < f[i + 1].first);
  assert(f[0].first == 0 && f[0].second != mint(0));

  mint f0inv = f[0].second.inv();
  std::vector<mint> g(deg);
  g[0] = f0inv;
  for (int i = 0; i < deg - 1; i++) {
    for (std::pair<int, mint> pim : f) {
      if (i + 1 - pim.first >= 0)
        g[i + 1] -= pim.second * g[i + 1 - pim.first];
      else
        continue;
    }
    g[i + 1] *= f0inv;
  }

  return g;
}

template <typename mint>
FPS<mint> inv_sparse(const FPS<mint>& f, int deg) {
  return inv_sparse(get_nonzeros(f), deg);
}

template <typename mint>
FPS<mint> FPS<mint>::inv_sparse(int deg) const {
  return ::inv_sparse(*this, deg);
}

// ↑--- inverse of sparse fps ----↑

// exp(f)のdeg次未満の部分を求める。
// F := exp(f) = F_0 + F_1 x + F_2 x^2 + ... とする。
// F' = F * f' なので
// F_1 + 2F_2 x + 3F_3 x^3 + ... = f' F.
// 0以上の整数iについて、i次の項に注目すると、
// (i+1) * F_{i+1} = [x^i] (f' * F)
// とわかる。Fは0,1,...,i次までわかってればF_{i+1}もわかるということになる。f'はスパースだからF_{i+1}はたかだかK回の計算で求められる.

template <typename mint>
FPS<mint> exp_sparse(const FPS<mint>& f, int deg) {
  FPS<mint> F(deg);
  F[0] = mint(1);

  std::vector<std::pair<int, mint>> nonzero_fdiff = get_nonzeros(f.diff());

  for (int i = 0; i + 1 < deg; i++) {
    // F[i+1]を求める
    // (i+1) * F_{i+1} = [x^i] (f' * F)

    for (std::pair<int, mint> pim : nonzero_fdiff) {
      int a = pim.first;
      // Fのi-a次の項を足していく
      if (i - a < 0) continue;
      assert(i - a >= 0);
      assert(i + 1 > i - a);
      F[i + 1] += pim.second * F[i - a];
    }
    F[i + 1] /= mint(i + 1);
  }
  return F;
}

template <typename mint>
FPS<mint> FPS<mint>::exp_sparse(int deg) const {
  return ::exp_sparse(*this, deg);
}

template <typename mint>
FPS<mint> log_sparse(const FPS<mint>& f, int deg) {
  FPS<mint> f_inv = inv_sparse(f, deg);
  return multiply_sparse(f_inv, f.diff(), deg).integral().pre(deg);
}

template <typename mint>
FPS<mint> FPS<mint>::log_sparse(int deg) const {
  return ::log_sparse(*this, deg);
}

// g := f ^ k
// g' = k * f^{k-1} * f'
// fg' = k * f^k * f'
// fg' = k * g * f'

template <typename mint>
FPS<mint> pow_sparse(const FPS<mint>& f, long long k, int deg) {
  if (k == 0) {
    FPS ret = {mint(1)};
    ret.resize(deg);
    return ret;
  }

  if (f[0] == mint(0)) {
    int mindeg = 0;
    while (mindeg < deg && f[mindeg] == mint(0)) mindeg++;

    // (x^{mindeg})^k = x^{mindeg * k}
    // mindeg * k >= deg ⇔ k >= floor(deg / mindeg) である。
    // →: 自明 (k >= deg / mindeg >= floor(deg / mindeg) なので)

    // ←について:  h1: k >= floor(deg / mindeg) を仮定して Goal: k >= degを示す。
    // deg = mindeg * q + r (0 <= r < mindeg)と表す。これを使うと
    // h1: k >= q.
    // Goal: k >= mindeg * q + r.

    // mindeg * k > LLINF
    // mindeg > LLINF / k
    constexpr long long INF = 4450000000011100000;
    if (mindeg > INF / k || mindeg * k >= deg) {
      FPS<mint> ret(deg);
      assert(ret[0] == mint(0));
      return ret;
    }
    return pow_sparse(f >> mindeg, k, deg - mindeg * k) << k * mindeg;
  }

  FPS<mint> g(deg);
  assert(f[0] != mint(0));
  g[0] = f[0].pow(k);

  std::vector<std::pair<int, mint>> nonzero_f = get_nonzeros(f);

  for (int i = 0; i + 1 < deg; i++) {
    // g[0], g[1], ..., g[i]が判っている状態で,x^iに注目してg[i+1]を求めにいく。

    // fg' = (f[0] + f[1]x + f[2]x^2 + ... + f[i]x^i)(g[1] + 2g[2]x + 3g[3]x^2 + ... + ig[i]x^{i-1} + (i+1)g[i+1]x^i)
    // (左) kgf' = k(g[0] + g[1]x + g[2]x^2 + ... + g[i]x^i) * (f[1] + 2*f[2]x 3*f[3]x^2 + ... + i*f[i]x^{i-1} +
    // (i+1)*f[i+1]x^i) (右)

    // 左のx^iの係数は f[0](i+1)g[i+1] + f[1]ig[i] + f[2](i-1)g[i-1] + ... + f[i]g[1]
    // 右のx^iの係数は k * ( g[0]*(i+1)f[i+1] +  g[1]if[i] + ... + g[i]*1f[1])

    // f[0](i+1)g[i+1] = k * (g[0]*(i+1)f[i+1] + g[1]if[i] + ... + g[i]*1f[1]) - f[1]ig[i] - f[2](i-1)g[i-1] - ... -
    // f[i]g[1]

    mint sum(0);
    for (std::pair<int, mint> pim : nonzero_f) {
      // f[pim.first]: pim.second
      // 左では 0 <= pim.first <= i
      // 右では 1 <= pim.first <= i+1の部分を見る

      if (0 <= pim.first && pim.first <= i) sum -= pim.second * mint(i + 1 - pim.first) * g[i + 1 - pim.first];
      if (1 <= pim.first && pim.first <= i + 1) sum += mint(k) * g[i + 1 - pim.first] * mint(pim.first) * pim.second;
    }

    // for (int j=0; j<=i; j++) sum += mint(k) * g[j] * mint(i+1-j) * f[i+1-j];
    // for (int j=1; j<=i; j++) sum -= f[j] * mint(i+1-j) * g[i+1-j];

    g[i + 1] = sum / f[0] / mint(i + 1);
  }

  return g;
}

template <typename mint>
FPS<mint> FPS<mint>::pow_sparse(long long k, int deg) const {
  return ::pow_sparse(*this, k, deg);
}

// tabun baggute masu.
template <typename mint>
FPS<mint> multiply_sparse(const FPS<mint>& f, const std::vector<std::pair<int, mint>>& g, int deg = -1) {
  if (deg == -1) deg = f.size() - 1 + g.back().first + 1;

  FPS<mint> ret(deg);
  for (std::pair<int, mint> pim : g) {
    assert(pim.second != 0);
    if (pim.second == 0) continue;

    for (int i = 0; i < f.size(); i++) {
      if (i + pim.first >= ret.size()) continue;
      if (f[i] != mint(0) && pim.second != mint(0)) ret[i + pim.first] += pim.second * f[i];
    }
  }

  return ret;
}

template <typename mint>
FPS<mint> multiply_sparse(const FPS<mint>& f, const FPS<mint>& g, int deg = -1) {
  std::vector<std::pair<int, mint>> vpmi;

  for (int i = 0; i < g.size(); i++)
    if (g[i] != mint(0)) vpmi.emplace_back(i, g[i]);

  return multiply_sparse(f, vpmi, deg);
}

template <typename mint>
FPS<mint> FPS<mint>::mul_sparse(const FPS<mint>& g, int deg) const {
  return multiply_sparse(*this, g, deg);
}


#line 304 "formal-power-series/formal-power-series.hpp"

Back to top page