Formal Power Series
(src/polynomial/FormalPowerSeries.hpp)

View this file on GitHub
Last update: 2024-05-22 00:21:24+09:00
Include: #include "src/polynomial/FormalPowerSeries.hpp"
Link: https://arxiv.org/pdf/1301.5804.pdf

概要

形式的冪級数 (Formal Power Series) を効率的に扱うためのライブラリ． AC Library の convolution に準拠している． inv 等の関数について特に $\mathbb{F} _ {998244353} \lbrack x \rbrack$ 等 $\text{mod}$ が NTT-friendly な場合を想定した高速化のもとに実装しているため，$\mathbb{F} _ {10 ^ 9 + 7} \lbrack x \rbrack$ 等を扱う際は注意が必要である．

以下では

\[f(x) = \sum _ {i = 0} ^ {N - 1} a _ i x ^ i\]

に対する操作を想定する．

メンバ関数	効果	時間計算量
`differential()`	$f ^ \prime(x)$ を求める．	$\mathrm{O}(N)$
`integral()`	$\int f(x) \mathrm {d}x$ を求める．	$\mathrm{O}(N)$
`inv()`	$\frac{1}{f(x)}$ を求める．	$\mathrm{O}(N \log N)$
`log()`	$\log(f(x))$ を求める．	$\mathrm{O}(N \log N)$
`sqrt()`	$\sqrt{f(x)}$ を求める．	$\mathrm{O}(N \log N)$
`exp()`	$\exp(f(x))$ を求める．	$\mathrm{O}(N \log N)$
`pow(k)`	$f(x) ^ k$ を求める．	$\mathrm{O}(N \log N)$
`eval(T)`	$f(x)$ について $x = T$ を代入した際の値を計算する．	$\mathrm{O}(N \log N)$

メンバ関数 inv(), log(), sqrt(), exp(), pow() については引数に deg を指定することでその次数で計算を打ち切ることが可能である．これらの操作では高次の項がより次数の低い項に影響を及ぼすことはないため，打ち切りに問題はない．

各種アルゴリズムの詳細

乗算

$f(x) g(x)$ を求める．これは atcoder/convolution により $\mathrm{O}(N \log N)$ で計算可能である．

（多項式としての）除算

多項式 $f(x), g(x)$ について $f(x) = q(x) g(x) + r(x)$ なる多項式 $q(x), r(x)\ (\deg r < \deg g)$ を求める．以下，$\deg f = n - 1, \deg g = m - 1$ とすると $\deg q = n - m$ である．

$f$ について多項式 $f ^ R := f(x ^ {-1}) x ^ {n - 1}$ を定める．商と余りの定義式を $x \gets x ^ {-1}$ として，両辺に $x ^ {n - 1}$ を掛けて，

\[\begin{aligned} &&& f(x ^ {-1}) x ^ {n - 1} = q(x ^ {-1}) x ^ {n - m} \cdot g(x ^ {-1}) x ^ {m - 1} + r(x ^ {-1}) x ^ {m - 2} \cdot x ^ {n - m + 1} \\ && \iff & f ^ R = q ^ R g ^ R + r ^ R x ^ {n - m + 1} \\ && \implies & f ^ R = q ^ R g ^ R \pmod{x ^ {n - m + 1}}. \end{aligned}\]

$\deg q ^ R = n - m < n - m + 1$ より

\[q ^ R = q ^ R \bmod x ^ {n - m + 1} = \frac{f ^ R}{g ^ R}.\]

微分・積分

$f ^ \prime(x), \int f(x) \mathrm{d} x$ を求める．ただし，不定積分の積分定数は $0$ とする．これは naive に計算可能である．ここで，積分については $1$ から $N$ の $\mathbb{F} _ p$ における逆元を求める必要があるが，これは愚直に求めようとすると時間計算量として $\mathrm{O}(N \log \bmod)$ を要するが，$\mathrm{O}(N)$ で求める方法がある．

$x$ の $\mathbb{F} _ p$ における逆元を求めたい． $p$ を $x$ で割ったときの商を $q$，余りを $r (r < x)$ とする． $p = q x + r$ より，

\[\begin{aligned} qx + r \equiv 0 \pmod{p} &\iff q + r x ^ {-1} \equiv 0 \pmod{p} \\ &\iff r ^ {-1} q + x ^ {-1} \equiv 0 \pmod{p} \\ &\iff x ^ {-1} \equiv - r ^ {-1} q \pmod{p} \end{aligned}\]

$r < x$ だから，上式に基づいて逆元テーブルを昇順に構築できる．

inv

$f(x) g(x) \equiv 1 \pmod{x ^ N}$ なる $g(x) = \sum _ {i = 0} ^ {N - 1} b _ i x ^ i$ を求める．これを $\frac{1}{f(x)}$ とも表すことにする．

補題

$f g _ n \equiv 1 \pmod{x ^ n}$ なる $(n - 1)$ 次多項式 $g _ n$ を定める．このとき，

\[g _ {2 n} = 2 g _ n - f g _ n ^ 2\]

が成立する．

証明

$f g _ n \equiv 1 \pmod{x ^ n}$ より，$f g _ n - 1= h x ^ n$ と表せる．このとき，$(f g _ n - 1) ^ 2 = h ^ 2 x ^ {2 n}$ だから，

\[\begin{aligned} (f g _ n - 1) ^ 2 \equiv 0 \pmod{x ^ {2 n}} &\iff f ^ 2 g _ n ^ 2 - 2 f g _ n + 1 \equiv 0 \pmod{x ^ {2 n}} \\ &\iff f(2 g _ n - f g _ n ^ 2) \equiv 1 \pmod{x ^ {2 n}} \\ &\iff f g _ {2 n} \equiv 1 \pmod{x ^ {2 n}} \quad (g _ {2 n} = 2 g _ n - f g _ n ^ 2) \end{aligned}\]

となる．

$\blacksquare$

$a _ 0 \neq 0$ のとき，$g _ 1 = a _ 0 ^ {-1}$ だから帰納的にダブリングの要領で計算することが可能となる．計算量は

\[T(n) = T(n / 2) + \mathrm{O}(n \log n)\]

より $\mathrm{O}(N \log N)$ を得る．また，$a _ 0 \neq 0$ が逆元の存在の必要条件にもなる．

ここで，以下の性質に注目する．

$g _ {2 n}$ の $\lbrack 0, n )$ 次の係数は $g_n$ のそれに一致する．よって，$g _ {2 n}$ の計算においては $f g _ n ^ 2$ の $\lbrack n, 2 n )$ 次の係数が求まれば良い．
$f g _ n \equiv 1 \pmod{x ^ n}$ より $f g _ n$ の $\lbrack 1, n )$ 次の係数は $0$ であり，$f g _ n$ の $n$ 次未満の部分の多項式を $h$ とすると，$h$ は $0$ 次の多項式となり，$g _ n$ は $n - 1$ 次以下であるから $h g _ n$ は $n - 1$ 次以下となり，$g _ {2 n}$ の計算においては無視して良い．

これらより，以下のより高速なアルゴリズムが提案される．

$f$ の $\lbrack 0, n )$ 次の項を取り出したものを改めて $f$ と定める．
$h = f g _ n \equiv x ^ {2 n}$ を計算する．
$h$ の $\lbrack n, 2 n )$ 次の項を取り出したものを改めて $h$ と定める．
$- h g _ n$ を計算し，これの $\lbrack n, 2 n )$ 次の係数が $g _ {2 n}$ のそれに一致する．

毎回畳み込むのではなく FFT と IFFT をうまく組み合わせることで，高速化が可能である．

log

$\log (f(x))$，すなわち，

\[f(x) \equiv \sum _ {k = 0}^{N - 1} \frac{g(x) ^ k}{k!} \pmod{x ^ N}\]

なる $g(x)$ を求める． $g(x)$ の定数項を $0$ とすると，$k \geq N$ については $g(x) ^ k$ の次数が $N$ 以上となるので打ち切っても問題はない．逆にそうでなければ，$f(x)$ の定数項が不定となり計算不可能である．また，このとき $a _ 0 = 0 ^ 0 = 1$ である．

$\log(f(x))$ と $\int \frac{f ^ \prime(x)}{f(x)} \mathrm{d} x$ の $N$ 次未満の項は一致することが確認できる．これを素直に実装すれば良い．

sqrt

$f(x) \equiv g(x) ^ 2 \pmod{x ^ N}$ なる $g(x)$ を求める．これを $\sqrt{f(x)}$ とも表すことにする．

これも inv と同様の考え方で計算可能である．

$g _ n ^ 2 \equiv f \pmod{x ^ n}$ なる $(n - 1)$ 次多項式 $g _ n$ を定める． $a _ 0 \neq 0$ であれば $g _ n$ の 0 次の係数は $0$ ではないので $g _ n$ には逆元が存在する．このとき，

\[\begin{aligned} (g _ n ^ 2 - f) ^ 2 \equiv 0 \pmod{x ^ {2 n}} &\iff g _ n ^ 4 - 2 f g _ n ^ 2 + f ^ 2 \equiv 0 \pmod{x ^ {2 n}} \\ &\iff g _ n ^ 2 - 2 f + f ^ 2 g _ n ^ {-2} \equiv 0 \pmod{x ^ {2 n}} \\ &\iff 4 f \equiv g _ n ^ 2 + 2 f + f ^ 2 g _ n ^ {-2} \pmod{x ^ {2 n}} \\ &\iff f \equiv \left \lbrack \frac{1}{2}(g _ n + f g _ n ^ {-1}) \right \rbrack ^ 2 \pmod{x ^ {2 n}} \\ \end{aligned}\]

より，

\[g _ {2 n} = \frac{1}{2}(g _ n + f g _ n ^ {-1})\]

とすれば $g _ {2 n} ^ 2 \equiv f \pmod{x ^ {2 n}}$ が成立する． $g _ 1$ については別途離散対数問題を解くことで計算可能である．あとはダブリングの要領で計算でき，計算量は $g _ n ^ {-1}$ の計算がネックとなり，

\[T(n) = T(n / 2) + \mathrm{O}(n \log n)\]

より $\mathrm{O}(N \log N)$ を得る．

次に $a _ 0 = 0$ の場合について考える． $f$ の非零係数をもつ項の次数の最小値を $i$ とすると，$i$ が奇数であれば求めたい平方根が存在しないことは容易にわかる． $i$ が偶数の場合は $\lbrack i, n )$ 次の項を抽出し $x ^ {-i}$ をかけて詰めた多項式を改めて $f$ と定めれば，上記のアルゴリズムを適用することができ，その結果得られた多項式に $x ^ {i / 2}$ をかけたものが求めたい平方根となる．

exp

$\exp (f(x))$ ，すなわち，

\[g(x) \equiv \sum _ {k = 0} ^ {N - 1} \frac{f(x) ^ k}{k!} \pmod{x ^ N}\]

なる $g(x)$ を求める． $a _ 0 = 0$ とすると，$k \geq N$ については $f(x) ^ k$ の次数が $N$ 以上となるので打ち切っても問題はない．逆にそうでなければ，$g(x)$ の定数項が不定となり計算不可能である．

$\exp (f) \equiv g _ n \pmod{x ^ n}$ なる $(n - 1)$ 次多項式 $g _ n$ を定める．このとき，inv や sqrt と同様にして，

\[g _ {2 n} = g _ n(1 + f - \log g _ n)\]

が成立する．計算量は同じく $\mathrm{O}(N \log N)$ である．

これは一般の場合に計算可能な方法だが，$\bmod$ が NTT-friendly な場合には別のより高速な手法を用いることができる [1] ．

pow

$f(x) ^ k \pmod{x ^ N}$ を求める．繰り返し二乗法を用いると $\mathrm{O}(N \log N \log k)$ で計算可能であるが，$\mathrm{O}(N \log N)$ で計算することができる．

$a _ 0 = 1$ のとき，

\[f(x) ^ k \equiv \exp (k \log f(x)) \pmod{x ^ N}\]

が成立することを利用する．定数項が $1$ の多項式 $g(x)$ を用いて $f(x) = c x ^ l g(x)\ (l \geq 0)$ と表せる．このとき，$f(x) ^ k = c ^ k x ^ {lk} g(x) ^ k$ であることを利用すれば良い．

inv（$f$ が sparse な場合）

$f$ の非零の個数を $K$ とすると，$1 / f(x)$ の先頭 $N$ 項を $\mathrm{O}(N K)$ で列挙することができる．

$f g = 1$ より $f _ 0 g _ 0 = 1$ かつ $n \geq 1$ について $\sum _ {i = 0} ^ n f _ i g _ {n - i} = 0$ より $f _ 0 g _ n = - \sum _ {i = 1} ^ n f _ i g _ {n - i}$ で $n$ の昇順に計算可能である．

より一般に $g$ を疎な $f$ で割った $g / f$ を求めることも同様の計算量で可能である．

log（$f$ が sparse な場合）

$f$ の非零の個数を $K$ とすると，$\log f$ の先頭 $N$ 項を $\mathrm{O}(N K)$ で列挙することができる．

$\log f = f ^ \prime / f$ であるから inv の場合と同様にして計算できる．

exp（$f$ が sparse な場合）

$f$ の非零の個数を $K$ とすると，$\exp f$ の先頭 $N$ 項を $\mathrm{O}(N K)$ で列挙することができる．

疎でない場合の exp と同様に $f _ 0 = 0$ である． $F = \exp f$ とすると，$F _ 0 = 1$ で，両辺微分して $F ^ \prime = F f ^ \prime$． $f ^ \prime$ も疎であるから，$F _ n$ までわかっているとき，右辺の $n$ 次の係数，すなわち $F ^ \prime$ の $n$ 次の係数が計算できるから $F$ の $n + 1$ 次の係数がわかる．

pow（$f$ が sparse な場合）

$f$ の非零の個数を $K$ とすると，$f(x) ^ M$ の先頭 $N$ 項を $\mathrm{O}(N K + \log M)$ で列挙することができる．

$f$ が定数項をもつ場合に帰着する． $F = f ^ M$ とすると，両辺微分して $F ^ \prime = M f ^ {M - 1} f ^ \prime$ より $f F ^ \prime = M f ^ \prime F$． $g = F ^ \prime$ として $n$ 次の係数について，$\sum _ {i = 0} ^ n f _ i g _ {n - i} = M \sum _ {i = 1} ^ n i F _ {n - i + 1} f _ i$ であり，$f _ 0 g _ n = - \sum _ {i = 1} ^ n f _ i g _ {n - i} + M \sum _ {i = 1} ^ n i F _ {n - i + 1} f _ i$． $F _ n$ までわかっているとき，右辺が計算可能であるから $g _ n$ がわかり，$g = F ^ \prime$ より $F _ {n + 1} = \frac{g _ n}{(n + 1) f _ 0}$．

sqrt（$f$ が sparse な場合）

$f$ の非零の個数を $K$ とすると，$\sqrt{f(x)}$ の先頭 $N$ 項を $\mathrm{O}(N K)$ で列挙することができる．

$f$ の定数項が $1$ の場合に帰着すると，$\sqrt{f} = f ^ {1 / 2}$ が成立し，pow と同様にして計算可能である．

以上の $f$ が sparse な場合の各種演算は Number Theoretic Transform による積の高速化を必要としないため，$\mathbb{F} _ {10 ^ 9 + 7}[x]$ 等，法が NTT-friendly でない場合にも適用可能である．

Polynomial Taylor Shift

与えられた多項式 $f(x)$ 及び定数 $c$ について，$f(x + c)$ を求める．

\[\begin{aligned} f(x + c) & = \sum _ {i = 0} ^ {N - 1} a _ i (x + c) ^ i \\ & = \sum _ {i = 0} ^ {N - 1} x ^ i \sum _ {j = i} ^ {N - 1} a _ j \binom{j}{i} c ^ {j - i} \\ & = \sum _ {i = 0} ^ {N - 1} \frac{x ^ i}{i!} \sum _ {j = i} ^ {N - 1} \frac{c ^ {j - i}}{(j - i)!} a _ j j! \\ & = \sum _ {i = 0} ^ {N - 1} \frac{x ^ i}{i!} \sum _ {j = 0} ^ {N - 1 - i} \frac{c ^ j}{j!} a _ {i + j} (i + j)! \\ & = \sum _ {k = 0} ^ {N - 1} a _ k k! \sum _ {j = 0} ^ {k} \frac{c ^ j}{j!} \frac{x ^ {k - j}}{(k - j)!} \end{aligned}\]

となり，畳込みに帰着される．計算量は $\mathrm{O}(N \log N)$．

出題例

Codeforces Round #250 (Div. 1) E. The Child and Binary Tree
yukicoder No.1145 Sums of Powers
- 答えの母関数を考えると，
  \[\begin{align*} \sum _ {n = 0} ^ \infty x ^ n \sum _ {i = 1} ^ N A _ i ^ n &= \sum _ {i = 1} ^ N \sum _ {n = 0} ^ \infty (A _ i x) ^ n \\ &= \sum _ {i = 1} ^ N \frac{1}{1 - A _ i x} \end{align*}\]
  となり，これは分子と分母を二分木のように計算していくことで求めることができる．
  
  もしくは，
  \[- \log (1 - x) = \sum _ {n = 1} ^ \infty \frac{x ^ n}{n}\]
  を用いると
  \[\begin{align*} - \log \prod _ {i = 1} ^ N (1 - A _ i x) &= \sum _ {i = 1} ^ N - \log (1 - A _ i x) \\ &= \sum _ {i = 1} ^ N \sum _ {n = 0} ^ \infty \frac{(A _ i x) ^ n}{n} \\ &= \sum _ {n = 0} ^ \infty \left(\sum _ {i = 1} ^ N A _ i ^ n \right) \frac{x ^ n}{n} \end{align*}\]
  となることからも計算可能である．
AtCoder Beginner Contest 222 H - Beautiful Binary Tree
AtCoder Beginner Contest 230 H - Bullion
AtCoder Beginner Contest 260 Ex - Colorfulness
- 包除原理により $C(P) = n$ である $P$ の個数を求めた後，Sums of Powers と同様の考察により解くことができる．
灘校文化祭コンテスト 2022 Day1 K - Li
AtCoder Regular Contest 154 F - Dice Game
- $\sum _ k a _ k e ^ {b _ k x}$ が計算できる．具体的には求める多項式を $f$ として $f$ の $i$ 次の係数を $i!$ 倍した多項式を $g$ とすると，これは $e ^ {b x} = \sum _ {i \geq 0} \frac{(b x) ^ i}{i!}$ より $\sum _ k a _ k \sum _ {i \geq 0} (b _ k x) ^ i = \sum _ k \frac{a _ k}{1 - b _ k x}$ に等しい．
- このような計算テクニックはここに詳しい．
yukicoder No.2514 Twelvefold Way Returns
- $k$ の倍数次の係数を取り出したい際には $1$ の $k$ 乗根に注目すると良い（参照 1 (Roots of Unity Filter)，参照 2）．
AtCoder Beginner Contest 387 G - Prime Circuit
- 畳み込みにおいて二項係数が関わってくる際には指数型母関数が有効．

Reference

[1] A simple and fast algorithm for computing exponentials of power series

Required by

Verified with

Code

#pragma once
#include <algorithm>
#include <cassert>
#include <functional>
#include <queue>
#include <utility>
#include <vector>

#include "../atcoder/convolution"

template <typename T> struct FormalPowerSeries : std::vector<T> {
  private:
    using std::vector<T>::vector;
    using FPS = FormalPowerSeries;
    void shrink() {
        while (this->size() and this->back() == T(0)) this->pop_back();
    }

    FPS pre(size_t sz) const { return FPS(this->begin(), this->begin() + std::min(this->size(), sz)); }

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

    FPS operator>>(size_t sz) const {
        if (this->size() <= sz) return {};
        return FPS(this->begin() + sz, this->end());
    }

    FPS operator<<(size_t sz) const {
        if (this->empty()) return {};
        FPS ret(*this);
        ret.insert(ret.begin(), sz, T(0));
        return ret;
    }

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

    FPS& operator+=(const T& v) {
        if (this->empty()) this->resize(1);
        (*this)[0] += v;
        shrink();
        return *this;
    }

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

    FPS& operator-=(const T& v) {
        if (this->empty()) this->resize(1);
        (*this)[0] -= v;
        shrink();
        return *this;
    }

    FPS& operator*=(const FPS& r) {
        auto res = atcoder::convolution(*this, r);
        return *this = {res.begin(), res.end()};
    }

    FPS& operator*=(const T& v) {
        for (auto& x : (*this)) x *= v;
        shrink();
        return *this;
    }

    FPS& operator/=(const FPS& r) {
        if (this->size() < r.size()) {
            this->clear();
            return *this;
        }
        int n = this->size() - r.size() + 1;
        return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev();
    }

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

    FPS operator+(const FPS& r) const { return FPS(*this) += r; }

    FPS operator+(const T& v) const { return FPS(*this) += v; }

    FPS operator-(const FPS& r) const { return FPS(*this) -= r; }

    FPS operator-(const T& v) const { return FPS(*this) -= v; }

    FPS operator*(const FPS& r) const { return FPS(*this) *= r; }

    FPS operator*(const T& v) const { return FPS(*this) *= v; }

    FPS operator/(const FPS& r) const { return FPS(*this) /= r; }

    FPS operator%(const FPS& r) const { return FPS(*this) %= r; }

    FPS operator-() const {
        FPS ret = *this;
        for (auto& v : ret) v = -v;
        return ret;
    }

    FPS differential() const {
        const int n = (int)this->size();
        FPS ret(std::max(0, n - 1));
        for (int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i);
        return ret;
    }

    FPS integral() const {
        const int n = (int)this->size();
        FPS ret(n + 1);
        ret[0] = T(0);
        if (n > 0) ret[1] = T(1);
        auto mod = T::mod();
        for (int i = 2; i <= n; i++) ret[i] = -ret[mod % i] * (mod / i);
        for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i];
        return ret;
    }

    FPS inv(int deg = -1) const {
        assert((*this)[0] != T(0));
        const int n = (int)this->size();
        if (deg == -1) deg = n;
        FPS ret{(*this)[0].inv()};
        ret.reserve(deg);
        for (int d = 1; d < deg; d <<= 1) {
            FPS f(d << 1), g(d << 1);
            std::copy(this->begin(), this->begin() + std::min(n, d << 1), f.begin());
            std::copy(ret.begin(), ret.end(), g.begin());
            atcoder::internal::butterfly(f);
            atcoder::internal::butterfly(g);
            for (int i = 0; i < (d << 1); i++) f[i] *= g[i];
            atcoder::internal::butterfly_inv(f);
            std::fill(f.begin(), f.begin() + d, T(0));
            atcoder::internal::butterfly(f);
            for (int i = 0; i < (d << 1); i++) f[i] *= g[i];
            atcoder::internal::butterfly_inv(f);
            T iz = T(d << 1).inv();
            iz *= -iz;
            for (int i = d; i < std::min(d << 1, deg); i++) ret.push_back(f[i] * iz);
        }
        return ret.pre(deg);
    }

    FPS log(int deg = -1) const {
        assert((*this)[0] == T(1));
        if (deg == -1) deg = (int)this->size();
        return (differential() * inv(deg)).pre(deg - 1).integral();
    }

    FPS sqrt(const std::function<T(T)>& get_sqrt, int deg = -1) const {
        const int n = this->size();
        if (deg == -1) deg = n;
        if (this->empty()) return FPS(deg, 0);
        if ((*this)[0] == T(0)) {
            for (int i = 1; i < n; i++) {
                if ((*this)[i] != T(0)) {
                    if (i & 1) return {};
                    if (deg - i / 2 <= 0) break;
                    auto ret = (*this >> i).sqrt(get_sqrt, deg - i / 2);
                    if (ret.empty()) return {};
                    ret = ret << (i / 2);
                    if ((int)ret.size() < deg) ret.resize(deg, T(0));
                    return ret;
                }
            }
            return FPS(deg, T(0));
        }
        auto sqrtf0 = T(get_sqrt((*this)[0]));
        if (sqrtf0 * sqrtf0 != (*this)[0]) return {};
        FPS ret{sqrtf0};
        T inv2 = T(2).inv();
        for (int i = 1; i < deg; i <<= 1) ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2;
        return ret.pre(deg);
    }

    /**
     * @brief Exp of Formal Power Series
     *
     * @see https://arxiv.org/pdf/1301.5804.pdf
     */
    FPS exp(int deg = -1) const {
        assert(this->empty() or (*this)[0] == T(0));
        if (this->size() <= 1) return {T(1)};
        if (deg == -1) deg = (int)this->size();
        FPS inv;
        inv.reserve(deg + 1);
        inv.push_back(T(0));
        inv.push_back(T(1));
        auto inplace_integral = [&](FPS& F) -> void {
            const int n = (int)F.size();
            auto mod = T::mod();
            while ((int)inv.size() <= n) {
                int i = inv.size();
                inv.push_back(-inv[mod % i] * (mod / i));
            }
            F.insert(F.begin(), T(0));
            for (int i = 1; i <= n; i++) F[i] *= inv[i];
        };
        auto inplace_differential = [](FPS& F) -> void {
            if (F.empty()) return;
            F.erase(F.begin());
            for (size_t i = 0; i < F.size(); i++) F[i] *= T(i + 1);
        };
        FPS f{1, (*this)[1]}, g{T(1)}, g_fft{T(1), T(1)};
        for (int m = 2; m < deg; m <<= 1) {
            const T iz1 = T(m).inv(), iz2 = T(m << 1).inv();
            auto f_fft = f;
            f_fft.resize(m << 1);
            atcoder::internal::butterfly(f_fft);
            {
                // Step 2.a'
                FPS _g(m);
                for (int i = 0; i < m; i++) _g[i] = f_fft[i] * g_fft[i];
                atcoder::internal::butterfly_inv(_g);
                std::fill(_g.begin(), _g.begin() + (m >> 1), T(0));
                atcoder::internal::butterfly(_g);
                for (int i = 0; i < m; i++) _g[i] *= -g_fft[i] * iz1 * iz1;
                atcoder::internal::butterfly_inv(_g);
                g.insert(g.end(), _g.begin() + (m >> 1), _g.end());

                g_fft = g;
                g_fft.resize(m << 1);
                atcoder::internal::butterfly(g_fft);
            }
            FPS x(this->begin(), this->begin() + std::min((int)this->size(), m));
            {
                // Step 2.b'
                x.resize(m);
                inplace_differential(x);
                x.push_back(T(0));
                atcoder::internal::butterfly(x);
            }
            {
                // Step 2.c'
                for (int i = 0; i < m; i++) x[i] *= f_fft[i] * iz1;
                atcoder::internal::butterfly_inv(x);
            }
            {
                // Step 2.d' and 2.e'
                x -= f.differential();
                x.resize(m << 1);
                for (int i = 0; i < m - 1; i++) x[m + i] = x[i], x[i] = T(0);
                atcoder::internal::butterfly(x);
                for (int i = 0; i < (m << 1); i++) x[i] *= g_fft[i] * iz2;
                atcoder::internal::butterfly_inv(x);
            }
            {
                // Step 2.f'
                x.pop_back();
                inplace_integral(x);
                for (int i = m; i < std::min((int)this->size(), m << 1); i++) x[i] += (*this)[i];
                std::fill(x.begin(), x.begin() + m, T(0));
            }
            {
                // Step 2.g' and 2.h'
                atcoder::internal::butterfly(x);
                for (int i = 0; i < (m << 1); i++) x[i] *= f_fft[i] * iz2;
                atcoder::internal::butterfly_inv(x);
                f.insert(f.end(), x.begin() + m, x.end());
            }
        }
        return FPS{f.begin(), f.begin() + deg};
    }

    FPS pow(int64_t k, int deg = -1) const {
        const int n = (int)this->size();
        if (deg == -1) deg = n;
        if (k == 0) {
            auto res = FPS(deg, T(0));
            res[0] = T(1);
            return res;
        }
        for (int i = 0; i < n; i++) {
            if ((*this)[i] != T(0)) {
                if (i >= (deg + k - 1) / k) return FPS(deg, T(0));
                T rev = (*this)[i].inv();
                FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg) * ((*this)[i].pow(k));
                ret = (ret << (i * k)).pre(deg);
                if ((int)ret.size() < deg) ret.resize(deg, T(0));
                return ret;
            }
        }
        return FPS(deg, T(0));
    }

    T eval(T x) const {
        T ret = 0, w = 1;
        for (const auto& v : *this) ret += w * v, w *= x;
        return ret;
    }

    static FPS product_of_polynomial_sequence(const std::vector<FPS>& fs) {
        if (fs.empty()) return {T(1)};
        auto comp = [](const FPS& f, const FPS& g) { return f.size() > g.size(); };
        std::priority_queue<FPS, std::vector<FPS>, decltype(comp)> pq{comp};
        for (const auto& f : fs) pq.emplace(f);
        while (pq.size() > 1) {
            auto f = pq.top();
            pq.pop();
            auto g = pq.top();
            pq.pop();
            pq.emplace(f * g);
        }
        return pq.top();
    }

    static FPS pow_sparse(const std::vector<std::pair<int, T>>& f, int64_t k, int n) {
        assert(k >= 0);
        int d = f.size(), offset = 0;
        while (offset < d and f[offset].second == 0) offset++;
        FPS res(n, 0);
        if (offset == d) {
            if (k == 0) res[0]++;
            return res;
        }
        if (f[offset].first > 0) {
            int deg = f[offset].first;
            if (k > (n - 1) / deg) return res;
            std::vector<std::pair<int, T>> g(f.begin() + offset, f.end());
            for (auto& p : g) p.first -= deg;
            auto tmp = pow_sparse(g, k, n - k * deg);
            for (int i = 0; i < n - k * deg; i++) res[k * deg + i] = tmp[i];
            return res;
        }
        std::vector<T> invs(n + 1);
        invs[0] = T(0);
        invs[1] = T(1);
        auto mod = T::mod();
        for (int i = 2; i <= n; i++) invs[i] = -invs[mod % i] * (mod / i);
        res[0] = f[0].second.pow(k);
        T coef = f[0].second.inv();
        for (int i = 1; i < n; i++) {
            for (int j = 1; j < d; j++) {
                if (i - f[j].first < 0) break;
                res[i] += f[j].second * res[i - f[j].first] * (T(k) * f[j].first - (i - f[j].first));
            }
            res[i] *= invs[i] * coef;
        }
        return res;
    }

    FPS taylor_shift(T c) const {
        FPS f(*this);
        const int n = f.size();
        std::vector<T> fac(n), finv(n);
        fac[0] = 1;
        for (int i = 1; i < n; i++) {
            fac[i] = fac[i - 1] * i;
            f[i] *= fac[i];
        }
        finv[n - 1] = fac[n - 1].inv();
        for (int i = n - 1; i > 0; i--) finv[i - 1] = finv[i] * i;
        std::reverse(f.begin(), f.end());
        FPS g(n);
        g[0] = T(1);
        for (int i = 1; i < n; i++) g[i] = g[i - 1] * c * finv[i] * fac[i - 1];
        f = (f * g).pre(n);
        std::reverse(f.begin(), f.end());
        for (int i = 0; i < n; i++) f[i] *= finv[i];
        return f;
    }
};

Traceback (most recent call last):
  File "/opt/hostedtoolcache/Python/3.13.3/x64/lib/python3.13/site-packages/onlinejudge_verify/documentation/build.py", line 71, in _render_source_code_stat
    bundled_code = language.bundle(stat.path, basedir=basedir, options={'include_paths': [basedir]}).decode()
                   ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.13.3/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/cplusplus.py", line 187, in bundle
    bundler.update(path)
    ~~~~~~~~~~~~~~^^^^^^
  File "/opt/hostedtoolcache/Python/3.13.3/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/cplusplus_bundle.py", line 401, in update
    self.update(self._resolve(pathlib.Path(included), included_from=path))
    ~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.13.3/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/cplusplus_bundle.py", line 401, in update
    self.update(self._resolve(pathlib.Path(included), included_from=path))
                ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.13.3/x64/lib/python3.13/site-packages/onlinejudge_verify/languages/cplusplus_bundle.py", line 260, in _resolve
    raise BundleErrorAt(path, -1, "no such header")
onlinejudge_verify.languages.cplusplus_bundle.BundleErrorAt: atcoder/convolution.hpp: line -1: no such header

cp-library

Formal Power Series
(src/polynomial/FormalPowerSeries.hpp)

概要

各種アルゴリズムの詳細

乗算

（多項式としての）除算

微分・積分

inv

補題

証明

log

sqrt

exp

pow

inv（$f$ が sparse な場合）

log（$f$ が sparse な場合）

exp（$f$ が sparse な場合）

pow（$f$ が sparse な場合）

sqrt（$f$ が sparse な場合）

Polynomial Taylor Shift

出題例

Reference

Links

Required by

Verified with

Code

Formal Power Series (src/polynomial/FormalPowerSeries.hpp)

概要

各種アルゴリズムの詳細

乗算

（多項式としての）除算

微分・積分

inv

補題

証明

log

sqrt

exp

pow

inv（$f$ が sparse な場合）

log（$f$ が sparse な場合）

exp（$f$ が sparse な場合）

pow（$f$ が sparse な場合）

sqrt（$f$ が sparse な場合）

Polynomial Taylor Shift

出題例

Reference

Links

Required by

Verified with

Code

Formal Power Series
(src/polynomial/FormalPowerSeries.hpp)