From bd5b080d42e4fc0ab7f8f9ee80609f6fe8da5a83 Mon Sep 17 00:00:00 2001 From: hitonanode <32937551+hitonanode@users.noreply.github.com> Date: Sun, 1 Mar 2026 22:00:22 +0900 Subject: [PATCH] sparse fps --- formal_power_series/pow_of_sparse_fps.hpp | 47 ----- formal_power_series/pow_of_sparse_fps.md | 43 ----- formal_power_series/sparse_fps.hpp | 174 ++++++++++++++++++ formal_power_series/sparse_fps.md | 65 +++++++ .../test/sparse_fps_exp.test.cpp | 24 +++ .../test/sparse_fps_inv.test.cpp | 24 +++ .../test/sparse_fps_log.test.cpp | 24 +++ ...est.cpp => sparse_fps_pow.stress.test.cpp} | 14 +- .../test/sparse_fps_pow.test.cpp | 25 +++ ...t.cpp => sparse_fps_pow.yuki1939.test.cpp} | 4 +- .../test/sparse_fps_sqrt.test.cpp | 27 +++ 11 files changed, 374 insertions(+), 97 deletions(-) delete mode 100644 formal_power_series/pow_of_sparse_fps.hpp delete mode 100644 formal_power_series/pow_of_sparse_fps.md create mode 100644 formal_power_series/sparse_fps.hpp create mode 100644 formal_power_series/sparse_fps.md create mode 100644 formal_power_series/test/sparse_fps_exp.test.cpp create mode 100644 formal_power_series/test/sparse_fps_inv.test.cpp create mode 100644 formal_power_series/test/sparse_fps_log.test.cpp rename formal_power_series/test/{pow_of_sparse_fps.stress.test.cpp => sparse_fps_pow.stress.test.cpp} (72%) create mode 100644 formal_power_series/test/sparse_fps_pow.test.cpp rename formal_power_series/test/{pow_of_sparse_fps.yuki1939.test.cpp => sparse_fps_pow.yuki1939.test.cpp} (79%) create mode 100644 formal_power_series/test/sparse_fps_sqrt.test.cpp diff --git a/formal_power_series/pow_of_sparse_fps.hpp b/formal_power_series/pow_of_sparse_fps.hpp deleted file mode 100644 index 5f256cca..00000000 --- a/formal_power_series/pow_of_sparse_fps.hpp +++ /dev/null @@ -1,47 +0,0 @@ -#pragma once -#include -#include -#include -#include -#include - -// Calculate f(x)^k up to x^N -// Requirement: k > 0 -// Complexity: O(NM) (M : num. of nonzero coefficients of f(x)) -template std::vector pow_of_sparse_fps(const std::vector &f, int N, long long k) { - assert(k > 0); - std::vector> nonzero_terms; - int d0 = 0; - while (d0 < int(f.size()) and d0 <= N and f[d0] == T()) ++d0; - if (d0 == int(f.size()) or d0 > N) return std::vector(N + 1); - - if (d0 and N / d0 < k) return std::vector(N + 1); - - for (int d = d0; d < std::min(N + 1, f.size()); ++d) { - if (f[d] != T()) nonzero_terms.emplace_back(d - d0, f[d]); - } - - int Ncalc = N - d0 * k; - - std::vector ret(Ncalc + 1); - ret[0] = f[d0].pow(k); - for (int d = 0; d + 1 < int(ret.size()); ++d) { - // Compare [x^d](k f'g - fg') - T tmp = T(); - int i, j; - T fi; - for (auto i_fi : nonzero_terms) { - std::tie(i, fi) = i_fi; - if (!i) continue; - j = d - i; - if (0 <= j) tmp -= fi * ret[j + 1] * (j + 1); - j = d - (i - 1); - if (0 <= j) tmp += fi * i * ret[j] * T(k); - } - ret[d + 1] = tmp / (f[d0] * (d + 1)); - } - - ret.resize(N + 1, T()); - std::rotate(ret.begin(), ret.end() - d0 * k, ret.end()); - return ret; -} diff --git a/formal_power_series/pow_of_sparse_fps.md b/formal_power_series/pow_of_sparse_fps.md deleted file mode 100644 index 2e02eebc..00000000 --- a/formal_power_series/pow_of_sparse_fps.md +++ /dev/null @@ -1,43 +0,0 @@ ---- -title: Power of sparse formal power series （非零な項が疎な形式的冪級数の累乗） -documentation_of: ./pow_of_sparse_fps.hpp ---- - -形式的冪級数 $f(x)$ について，$(f(x))^k$ の $x^N$ までの項を $O(NM)$ で計算する．ここで $M$ は $f(x)$ の $N$ 次以下で非零な項の個数． - -## ナイーブな説明 - -$\displaystyle -g(x) = f(x)^k -$ - -について，両辺 $\log$ をとってから微分することで - -$\displaystyle -\frac{g'(x)}{g(x)} = k \frac{f'(x)}{f(x)} -$ - -を得る．これより - -$\displaystyle -f(x) g'(x) = k f'(x) g(x) -$ - -が成立する．両辺の各項の次数を低い順から順に合わせていけばよい． - -## 使用方法 - -```cpp -vector f; // mint は除算が可能でなければならない -long long k; - -// f(x)^k を x^N の係数まで計算 -vector g = pow_of_sparse_fps(f, N, k) // g.size() == N + 1 -``` - -## 問題例 - -- [No.1939 Numbered Colorful Balls - yukicoder](https://yukicoder.me/problems/no/1939) - -## リンク -- [A problem collection of ODE and differential technique - Codeforces](https://codeforces.com/blog/entry/76447) diff --git a/formal_power_series/sparse_fps.hpp b/formal_power_series/sparse_fps.hpp new file mode 100644 index 00000000..d934d665 --- /dev/null +++ b/formal_power_series/sparse_fps.hpp @@ -0,0 +1,174 @@ +#pragma once +#include +#include +#include +#include +#include +#include + +namespace sparse_fps { +// https://github.com/yosupo06/library-checker-problems/issues/767#issuecomment-1166414020 + +// Calculate f(x)^k up to x^max_deg +template + requires std::derived_from> +Vec pow(const Vec &f, int max_deg, long long k) { + using T = typename Vec::value_type; + assert(k >= 0); + + Vec ret(max_deg + 1); + + if (k == 0) { + ret[0] = T{1}; + return ret; + } + + std::vector> terms; + int d0 = 0; + while (d0 < int(f.size()) and d0 <= max_deg and f[d0] == T()) ++d0; + if (d0 == int(f.size()) or d0 > max_deg) return ret; + + if (d0 and max_deg / d0 < k) return ret; + + for (int d = d0 + 1; d < std::min(max_deg + 1, f.size()); ++d) { + if (f[d] != T{}) terms.emplace_back(d - d0, f[d]); + } + + const int bias = d0 * k; + ret[bias] = f[d0].pow(k); + const T fd0inv = 1 / f[d0]; + for (int d = 0; bias + d + 1 < int(ret.size()); ++d) { + // Compare [x^d](k f'g - fg') + T tmp{0}; + for (auto [i, fi] : terms) { + int j = d - i; + if (0 <= j) tmp -= fi * ret[bias + j + 1] * (j + 1); + j = d - (i - 1); + if (0 <= j) tmp += fi * i * ret[bias + j] * T(k); + } + ret[bias + d + 1] = tmp * fd0inv / (d + 1); + } + + return ret; +} + +template + requires std::derived_from> +Vec inv(const Vec &f, int max_deg) { + using T = typename Vec::value_type; + assert(!f.empty() and f[0] != T{0}); + + Vec ret(max_deg + 1); + + std::vector> terms; + for (int d = 1; d < (int)f.size() and d <= max_deg; ++d) { + if (f[d] != T{}) terms.emplace_back(d, f[d]); + } + + const T f0inv = f[0].inv(); + ret[0] = f0inv; + + for (int d = 1; d <= max_deg; ++d) { + T tmp{0}; + for (auto [i, fi] : terms) { + if (i > d) break; + tmp -= fi * ret[d - i]; + } + ret[d] = tmp * f0inv; + } + + return ret; +} + +template + requires std::derived_from> +Vec log(const Vec &f, int max_deg) { + using T = typename Vec::value_type; + assert(!f.empty() and f[0] != T{0}); + + const auto inv = sparse_fps::inv(f, max_deg); + + std::vector> df_terms; + for (int d = 1; d < (int)f.size() and d <= max_deg; ++d) { + if (f[d] != T{}) df_terms.emplace_back(d - 1, f[d] * T{d}); + } + + Vec ret(max_deg + 1); + + for (int d = 0; d < max_deg; ++d) { + for (auto [i, fi] : df_terms) { + const int j = d + i + 1; + if (j > max_deg) break; + ret[j] += inv[d] * fi * T{j}.inv(); + } + } + + return ret; +} + +template + requires std::derived_from> +Vec exp(const Vec &f, int max_deg) { + using T = typename Vec::value_type; + assert(f.empty() or f[0] == T{0}); + + std::vector> df_terms; + for (int d = 1; d < (int)f.size() and d <= max_deg; ++d) { + if (f[d] != T{}) df_terms.emplace_back(d - 1, f[d] * T{d}); + } + + Vec ret(max_deg + 1); + ret[0] = T{1}; + // F' = F * f' + for (int d = 1; d <= max_deg; ++d) { + T tmp{0}; + for (auto [i, dfi] : df_terms) { + if (i > d - 1) break; + tmp += dfi * ret[d - 1 - i]; + } + ret[d] = tmp * T{d}.inv(); + } + + return ret; +} + +template + requires std::derived_from> +std::optional sqrt(const Vec &f, int max_deg) { + using T = typename Vec::value_type; + + Vec ret(max_deg + 1); + + int d0 = 0; + while (d0 < int(f.size()) and d0 <= max_deg and f[d0] == T{}) ++d0; + if (d0 == int(f.size()) or d0 > max_deg) return ret; + if (d0 & 1) return std::nullopt; + + const T sqrtf0 = f[d0].sqrt(); + if (sqrtf0 == T{}) return std::nullopt; + + std::vector> terms; + const T fd0inv = 1 / f[d0]; + for (int d = d0 + 1; d < std::min(max_deg + 1, f.size()); ++d) { + if (f[d] != T{}) terms.emplace_back(d - d0, f[d] * fd0inv); + } + + const int bias = d0 / 2; + const T inv2 = T{2}.inv(); + ret[bias] = sqrtf0; + for (int d = 0; bias + d + 1 < int(ret.size()); ++d) { + T tmp{0}; + for (auto [i, fi] : terms) { + if (i > d + 1) break; + int j = d - i; + if (0 <= j) tmp -= fi * ret[bias + j + 1] * (j + 1); + j = d - (i - 1); + if (0 <= j) tmp += fi * i * ret[bias + j] * inv2; + } + ret[bias + d + 1] = tmp / (d + 1); + } + + return ret; +} + +} // namespace sparse_fps diff --git a/formal_power_series/sparse_fps.md b/formal_power_series/sparse_fps.md new file mode 100644 index 00000000..4f6c032e --- /dev/null +++ b/formal_power_series/sparse_fps.md @@ -0,0 +1,65 @@ +--- +title: Sparse formal power series （疎な形式的冪級数の演算） +documentation_of: ./sparse_fps.hpp +--- + +疎な形式的冪級数 $f(x)$ に対して，$x^N$ までの累乗・逆元・対数・指数関数・平方根を計算する．$f(x)$ の非零項数を $K$ として，いずれも $O(NK)$ で動作する． + +## 使用方法 + +名前空間 `sparse_fps` に以下の関数が定義されている．テンプレート引数 `Vec` は `std::vector` を継承した型で，`T` は体の元として `+`, `-`, `*`, `/`, `pow`, `inv` 等をサポートする必要がある（典型的には ModInt）． + +### `sparse_fps::pow` + +```cpp +Vec sparse_fps::pow(const Vec &f, int max_deg, long long k); +``` + +$f(x)^k \bmod x^{\mathrm{max\_deg}+1}$ を計算する．$k \ge 0$． + +### `sparse_fps::inv` + +```cpp +Vec sparse_fps::inv(const Vec &f, int max_deg); +``` + +$1 / f(x) \bmod x^{\mathrm{max\_deg}+1}$ を計算する．$f_0 \neq 0$ が必要． + +### `sparse_fps::log` + +```cpp +Vec sparse_fps::log(const Vec &f, int max_deg); +``` + +$\log f(x) \bmod x^{\mathrm{max\_deg}+1}$ を計算する．$f_0 \neq 0$ が必要． + +### `sparse_fps::exp` + +```cpp +Vec sparse_fps::exp(const Vec &f, int max_deg); +``` + +$\exp f(x) \bmod x^{\mathrm{max\_deg}+1}$ を計算する．$f_0 = 0$ が必要． + +### `sparse_fps::sqrt` + +```cpp +std::optional sparse_fps::sqrt(const Vec &f, int max_deg); +``` + +$\sqrt{f(x)} \bmod x^{\mathrm{max\_deg}+1}$ を計算する．最小次数の非零項 $f_{d_0}$ について $d_0$ が偶数かつ $f_{d_0}$ が平方根を持つ必要がある．解が存在しない場合は `std::nullopt` を返す． + +## 問題例 + +- [Pow of Formal Power Series (Sparse) - Library Checker](https://judge.yosupo.jp/problem/pow_of_formal_power_series_sparse) +- [No.1939 Numbered Colorful Balls - yukicoder](https://yukicoder.me/problems/no/1939) +- [Codeforces Round 1058 (Div. 1) E. Super-Short-Polynomial-San](https://codeforces.com/contest/2159/problem/E) +- [Inv of Formal Power Series (Sparse) - Library Checker](https://judge.yosupo.jp/problem/inv_of_formal_power_series_sparse) +- [Log of Formal Power Series (Sparse) - Library Checker](https://judge.yosupo.jp/problem/log_of_formal_power_series_sparse) +- [Exp of Formal Power Series (Sparse) - Library Checker](https://judge.yosupo.jp/problem/exp_of_formal_power_series_sparse) +- [Sqrt of Formal Power Series (Sparse) - Library Checker](https://judge.yosupo.jp/problem/sqrt_of_formal_power_series_sparse) + +## リンク + +- [A problem collection of ODE and differential technique - Codeforces](https://codeforces.com/blog/entry/76447) +- [library-checker-problems #767 (comment)](https://github.com/yosupo06/library-checker-problems/issues/767#issuecomment-1166414020) diff --git a/formal_power_series/test/sparse_fps_exp.test.cpp b/formal_power_series/test/sparse_fps_exp.test.cpp new file mode 100644 index 00000000..cd3f5f95 --- /dev/null +++ b/formal_power_series/test/sparse_fps_exp.test.cpp @@ -0,0 +1,24 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/exp_of_formal_power_series_sparse" +#include "../sparse_fps.hpp" +#include "../../modint.hpp" +#include +#include +using namespace std; +using mint = ModInt<998244353>; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + int N, K; + cin >> N >> K; + vector f(N); + while (K--) { + int i, a; + cin >> i >> a; + f.at(i) = a; + } + + const auto ret = sparse_fps::exp(f, N - 1); + + for (auto e : ret) cout << e << ' '; + cout << '\n'; +} diff --git a/formal_power_series/test/sparse_fps_inv.test.cpp b/formal_power_series/test/sparse_fps_inv.test.cpp new file mode 100644 index 00000000..010f8102 --- /dev/null +++ b/formal_power_series/test/sparse_fps_inv.test.cpp @@ -0,0 +1,24 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/inv_of_formal_power_series_sparse" +#include "../sparse_fps.hpp" +#include "../../modint.hpp" +#include +#include +using namespace std; +using mint = ModInt<998244353>; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + int N, K; + cin >> N >> K; + vector f(N); + while (K--) { + int i, a; + cin >> i >> a; + f.at(i) = a; + } + + const auto ret = sparse_fps::inv(f, N - 1); + + for (auto e : ret) cout << e << ' '; + cout << '\n'; +} diff --git a/formal_power_series/test/sparse_fps_log.test.cpp b/formal_power_series/test/sparse_fps_log.test.cpp new file mode 100644 index 00000000..8d354fc5 --- /dev/null +++ b/formal_power_series/test/sparse_fps_log.test.cpp @@ -0,0 +1,24 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/log_of_formal_power_series_sparse" +#include "../sparse_fps.hpp" +#include "../../modint.hpp" +#include +#include +using namespace std; +using mint = ModInt<998244353>; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + int N, K; + cin >> N >> K; + vector f(N); + while (K--) { + int i, a; + cin >> i >> a; + f.at(i) = a; + } + + const auto ret = sparse_fps::log(f, N - 1); + + for (auto e : ret) cout << e << ' '; + cout << '\n'; +} diff --git a/formal_power_series/test/pow_of_sparse_fps.stress.test.cpp b/formal_power_series/test/sparse_fps_pow.stress.test.cpp similarity index 72% rename from formal_power_series/test/pow_of_sparse_fps.stress.test.cpp rename to formal_power_series/test/sparse_fps_pow.stress.test.cpp index ab651c6b..2f1340ad 100644 --- a/formal_power_series/test/pow_of_sparse_fps.stress.test.cpp +++ b/formal_power_series/test/sparse_fps_pow.stress.test.cpp @@ -1,5 +1,5 @@ #define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A" // DUMMY -#include "../pow_of_sparse_fps.hpp" +#include "../sparse_fps.hpp" #include "../../modint.hpp" #include "../../random/xorshift.hpp" #include "../formal_power_series.hpp" @@ -20,12 +20,16 @@ int main() { vector f(nin); for (int i = d0; i < nin; ++i) f[i] = rand_int(); - auto ret = pow_of_sparse_fps(f, degout, k); - auto ret_fps = fps(f.begin(), f.end()).pow(k, degout + 1); + const fps init = fps(f.begin(), f.end()); + + const vector ret1 = sparse_fps::pow(f, degout, k); + const fps ret2 = sparse_fps::pow(init, degout, k); + const fps ret_fps = init.pow(k, degout + 1); + for (int i = 0; i <= degout; ++i) { - if (ret[i] != ret_fps.coeff(i)) { + if (ret1[i] != ret_fps.coeff(i) or ret2[i] != ret_fps.coeff(i)) { cerr << nin << ' ' << d0 << ' ' << degout << ' ' << k << ' ' - << iter << ' ' << i << ' ' << ret[i] << ' ' + << iter << ' ' << i << ' ' << ret1[i] << ' ' << ret2[i] << ' ' << ret_fps.coeff(i) << endl; exit(1); } diff --git a/formal_power_series/test/sparse_fps_pow.test.cpp b/formal_power_series/test/sparse_fps_pow.test.cpp new file mode 100644 index 00000000..9b7b85a6 --- /dev/null +++ b/formal_power_series/test/sparse_fps_pow.test.cpp @@ -0,0 +1,25 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/pow_of_formal_power_series_sparse" +#include "../sparse_fps.hpp" +#include "../../modint.hpp" +#include +#include +using namespace std; +using mint = ModInt<998244353>; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + int N, K; + long long M; + cin >> N >> K >> M; + vector f(N); + while (K--) { + int i, a; + cin >> i >> a; + f.at(i) = a; + } + + const auto ret = sparse_fps::pow(f, N - 1, M); + + for (auto e : ret) cout << e << ' '; + cout << '\n'; +} diff --git a/formal_power_series/test/pow_of_sparse_fps.yuki1939.test.cpp b/formal_power_series/test/sparse_fps_pow.yuki1939.test.cpp similarity index 79% rename from formal_power_series/test/pow_of_sparse_fps.yuki1939.test.cpp rename to formal_power_series/test/sparse_fps_pow.yuki1939.test.cpp index 09191fd9..bf166b15 100644 --- a/formal_power_series/test/pow_of_sparse_fps.yuki1939.test.cpp +++ b/formal_power_series/test/sparse_fps_pow.yuki1939.test.cpp @@ -1,5 +1,5 @@ #define PROBLEM "https://yukicoder.me/problems/no/1939" -#include "../pow_of_sparse_fps.hpp" +#include "../sparse_fps.hpp" #include "../../modint.hpp" #include #include @@ -17,5 +17,5 @@ int main() { cin >> l; f[l] += 1; } - cout << pow_of_sparse_fps(f, N, N + 1)[N] / (N + 1) << '\n'; + cout << sparse_fps::pow(f, N, N + 1)[N] / (N + 1) << '\n'; } diff --git a/formal_power_series/test/sparse_fps_sqrt.test.cpp b/formal_power_series/test/sparse_fps_sqrt.test.cpp new file mode 100644 index 00000000..d4ca8c80 --- /dev/null +++ b/formal_power_series/test/sparse_fps_sqrt.test.cpp @@ -0,0 +1,27 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/sqrt_of_formal_power_series_sparse" +#include "../sparse_fps.hpp" +#include "../../modint.hpp" +#include +#include +using namespace std; +using mint = ModInt<998244353>; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + int N, K; + cin >> N >> K; + vector f(N); + while (K--) { + int i, a; + cin >> i >> a; + f.at(i) = a; + } + + const auto ret = sparse_fps::sqrt(f, N - 1); + if (!ret) { + cout << -1 << '\n'; + } else { + for (auto e : *ret) cout << e << ' '; + cout << '\n'; + } +}