sparse_formal_power_series.hpp¶
当非零项很少时求形式幂级数的逆，避免使用完整 NTT。
\[
\displaystyle f(x)g(x)=1
\]
Complexity: Time: O(nk), where k is the number of nonzero input terms. Space: O(n + k).
AC 记录：exp_of_formal_power_series_sparse, inv_of_formal_power_series_sparse, log_of_formal_power_series_sparse, pow_of_formal_power_series_sparse, sqrt_of_formal_power_series_sparse。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(nk), where k is the number of nonzero input terms.
/// Space: O(n + k).
#include "noya/mod_sqrt.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <optional>
#include <utility>
#include <vector>
namespace noya {
template <class Mint>
using sparse_power_series = std::vector<std::pair<int, Mint>>;
namespace sparse_fps_detail {
template <class Mint> std::vector<Mint> coefficient_inverses(int sz) {
std::vector<Mint> inv(sz);
if (sz > 1) {
inv[1] = Mint(1);
}
for (int vl = 2; vl < sz; vl++) {
inv[vl] = -Mint(Mint::mod() / vl) * inv[Mint::mod() % vl];
}
return inv;
}
template <class Mint>
void check_terms(int sz, const sparse_power_series<Mint> &ter) {
int pre = -1;
for (auto [deg, cf] : ter) {
assert(pre < deg && deg < sz && cf != Mint{});
pre = deg;
}
}
template <class Mint>
std::vector<Mint> unit_power(int sz, const sparse_power_series<Mint> &ter,
Mint exp) {
std::vector<Mint> res(sz);
if (sz == 0)
return res;
res[0] = Mint(1);
std::vector<Mint> inv = coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
vl += cf * res[deg - td] * (exp * Mint(td) - Mint(deg - td));
}
res[deg] = vl * inv[deg];
}
return res;
}
template <class Mint> Mint scalar_power(Mint vl, std::uint64_t exp) {
Mint res = Mint(1);
while (exp > 0) {
if (exp & 1)
res *= vl;
vl *= vl;
exp >>= 1;
}
return res;
}
} // namespace sparse_fps_detail
/// @brief Invert a sparse formal power series. From `f g = 1`, each new
/// coefficient of g is a dot product against the nonconstant terms of f.
template <class Mint>
std::vector<Mint> sparse_fps_inverse(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(!ter.empty() && ter.front().first == 0);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
Mint ic = Mint(1) / ter.front().second;
res[0] = ic;
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
vl += cf * res[deg - td];
}
res[deg] = -vl * ic;
}
return res;
}
/// @brief Compute log(f) for sparse f with f(0)=1. Comparing coefficients in
/// `f (log f)' = f'` gives one linear recurrence per output coefficient.
template <class Mint>
std::vector<Mint> sparse_fps_logarithm(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(!ter.empty() && ter.front().first == 0 &&
ter.front().second == Mint(1));
std::vector<Mint> res(sz);
std::vector<Mint> inv = sparse_fps_detail::coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
if (td == deg) {
vl += Mint(deg) * cf;
} else {
vl -= cf * Mint(deg - td) * res[deg - td];
}
}
res[deg] = vl * inv[deg];
}
return res;
}
/// @brief Compute exp(f) for sparse f with f(0)=0. The differential identity
/// `g' = f' g` gives the recurrence directly.
template <class Mint>
std::vector<Mint> sparse_fps_exponential(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(ter.empty() || ter.front().first > 0);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
res[0] = Mint(1);
std::vector<Mint> inv = sparse_fps_detail::coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td > deg)
break;
vl += Mint(td) * cf * res[deg - td];
}
res[deg] = vl * inv[deg];
}
return res;
}
/// @brief Raise a sparse series to a nonnegative integer exponent. The first
/// nonzero monomial determines the output shift and scale; after normalization
/// `f g' = exp f' g` yields an O(nk) recurrence.
template <class Mint>
std::vector<Mint> sparse_fps_power(int sz, const sparse_power_series<Mint> &ter,
std::uint64_t exp) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
if (exp == 0) {
res[0] = Mint(1);
return res;
}
if (ter.empty())
return res;
int fd = ter.front().first;
if (fd > 0 && exp > std::uint64_t((sz - 1) / fd)) {
return res;
}
int shf = int(std::uint64_t(fd) * exp);
int tar = sz - shf;
Mint ld = ter.front().second;
sparse_power_series<Mint> nrm;
nrm.reserve(ter.size());
for (auto [deg, cf] : ter) {
int nd = deg - fd;
if (nd >= tar)
break;
nrm.emplace_back(nd, cf / ld);
}
std::vector<Mint> uni = sparse_fps_detail::unit_power(tar, nrm, Mint(exp));
Mint scl = sparse_fps_detail::scalar_power(ld, exp);
for (int idx = 0; idx < tar; idx++) {
res[shf + idx] = scl * uni[idx];
}
return res;
}
/// @brief Return one square root of a sparse series when it exists. Its first
/// degree must be even and its leading coefficient a quadratic residue; after
/// removing both, the unit series is raised to the field exponent 1/2.
template <class Mint>
std::optional<std::vector<Mint>>
sparse_fps_square_root(int sz, const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
if (ter.empty())
return std::vector<Mint>(sz);
int fd = ter.front().first;
if (fd % 2 != 0)
return std::nullopt;
auto rt = mod_sqrt(std::uint64_t(ter.front().second.val()),
std::uint64_t(Mint::mod()));
if (!rt)
return std::nullopt;
int shf = fd / 2;
int tar = sz - fd;
Mint ld = ter.front().second;
sparse_power_series<Mint> nrm;
nrm.reserve(ter.size());
for (auto [deg, cf] : ter) {
int nd = deg - fd;
if (nd >= tar)
break;
nrm.emplace_back(nd, cf / ld);
}
std::vector<Mint> uni =
sparse_fps_detail::unit_power(tar, nrm, Mint(1) / Mint(2));
std::vector<Mint> res(sz);
for (int idx = 0; idx < tar; idx++) {
res[shf + idx] = Mint(*rt) * uni[idx];
}
return res;
}
} // namespace noya
#ifndef NOYA_SPARSE_FORMAL_POWER_SERIES_HPP
#define NOYA_SPARSE_FORMAL_POWER_SERIES_HPP 1
/// @complexity Time: O(nk), where k is the number of nonzero input terms.
/// Space: O(n + k).
#include "noya/mod_sqrt.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <optional>
#include <utility>
#include <vector>
namespace noya {
template <class Mint>
using sparse_power_series = std::vector<std::pair<int, Mint>>;
namespace sparse_fps_detail {
template <class Mint> std::vector<Mint> coefficient_inverses(int sz) {
std::vector<Mint> inv(sz);
if (sz > 1) {
inv[1] = Mint(1);
}
for (int vl = 2; vl < sz; vl++) {
inv[vl] = -Mint(Mint::mod() / vl) * inv[Mint::mod() % vl];
}
return inv;
}
template <class Mint>
void check_terms(int sz, const sparse_power_series<Mint> &ter) {
int pre = -1;
for (auto [deg, cf] : ter) {
assert(pre < deg && deg < sz && cf != Mint{});
pre = deg;
}
}
template <class Mint>
std::vector<Mint> unit_power(int sz, const sparse_power_series<Mint> &ter,
Mint exp) {
std::vector<Mint> res(sz);
if (sz == 0)
return res;
res[0] = Mint(1);
std::vector<Mint> inv = coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
vl += cf * res[deg - td] * (exp * Mint(td) - Mint(deg - td));
}
res[deg] = vl * inv[deg];
}
return res;
}
template <class Mint> Mint scalar_power(Mint vl, std::uint64_t exp) {
Mint res = Mint(1);
while (exp > 0) {
if (exp & 1)
res *= vl;
vl *= vl;
exp >>= 1;
}
return res;
}
} // namespace sparse_fps_detail
/// @brief Invert a sparse formal power series. From `f g = 1`, each new
/// coefficient of g is a dot product against the nonconstant terms of f.
template <class Mint>
std::vector<Mint> sparse_fps_inverse(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(!ter.empty() && ter.front().first == 0);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
Mint ic = Mint(1) / ter.front().second;
res[0] = ic;
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
vl += cf * res[deg - td];
}
res[deg] = -vl * ic;
}
return res;
}
/// @brief Compute log(f) for sparse f with f(0)=1. Comparing coefficients in
/// `f (log f)' = f'` gives one linear recurrence per output coefficient.
template <class Mint>
std::vector<Mint> sparse_fps_logarithm(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(!ter.empty() && ter.front().first == 0 &&
ter.front().second == Mint(1));
std::vector<Mint> res(sz);
std::vector<Mint> inv = sparse_fps_detail::coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
if (td == deg) {
vl += Mint(deg) * cf;
} else {
vl -= cf * Mint(deg - td) * res[deg - td];
}
}
res[deg] = vl * inv[deg];
}
return res;
}
/// @brief Compute exp(f) for sparse f with f(0)=0. The differential identity
/// `g' = f' g` gives the recurrence directly.
template <class Mint>
std::vector<Mint> sparse_fps_exponential(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(ter.empty() || ter.front().first > 0);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
res[0] = Mint(1);
std::vector<Mint> inv = sparse_fps_detail::coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td > deg)
break;
vl += Mint(td) * cf * res[deg - td];
}
res[deg] = vl * inv[deg];
}
return res;
}
/// @brief Raise a sparse series to a nonnegative integer exponent. The first
/// nonzero monomial determines the output shift and scale; after normalization
/// `f g' = exp f' g` yields an O(nk) recurrence.
template <class Mint>
std::vector<Mint> sparse_fps_power(int sz, const sparse_power_series<Mint> &ter,
std::uint64_t exp) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
if (exp == 0) {
res[0] = Mint(1);
return res;
}
if (ter.empty())
return res;
int fd = ter.front().first;
if (fd > 0 && exp > std::uint64_t((sz - 1) / fd)) {
return res;
}
int shf = int(std::uint64_t(fd) * exp);
int tar = sz - shf;
Mint ld = ter.front().second;
sparse_power_series<Mint> nrm;
nrm.reserve(ter.size());
for (auto [deg, cf] : ter) {
int nd = deg - fd;
if (nd >= tar)
break;
nrm.emplace_back(nd, cf / ld);
}
std::vector<Mint> uni = sparse_fps_detail::unit_power(tar, nrm, Mint(exp));
Mint scl = sparse_fps_detail::scalar_power(ld, exp);
for (int idx = 0; idx < tar; idx++) {
res[shf + idx] = scl * uni[idx];
}
return res;
}
/// @brief Return one square root of a sparse series when it exists. Its first
/// degree must be even and its leading coefficient a quadratic residue; after
/// removing both, the unit series is raised to the field exponent 1/2.
template <class Mint>
std::optional<std::vector<Mint>>
sparse_fps_square_root(int sz, const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
if (ter.empty())
return std::vector<Mint>(sz);
int fd = ter.front().first;
if (fd % 2 != 0)
return std::nullopt;
auto rt = mod_sqrt(std::uint64_t(ter.front().second.val()),
std::uint64_t(Mint::mod()));
if (!rt)
return std::nullopt;
int shf = fd / 2;
int tar = sz - fd;
Mint ld = ter.front().second;
sparse_power_series<Mint> nrm;
nrm.reserve(ter.size());
for (auto [deg, cf] : ter) {
int nd = deg - fd;
if (nd >= tar)
break;
nrm.emplace_back(nd, cf / ld);
}
std::vector<Mint> uni =
sparse_fps_detail::unit_power(tar, nrm, Mint(1) / Mint(2));
std::vector<Mint> res(sz);
for (int idx = 0; idx < tar; idx++) {
res[shf + idx] = Mint(*rt) * uni[idx];
}
return res;
}
} // namespace noya
#endif // NOYA_SPARSE_FORMAL_POWER_SERIES_HPP
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <optional>
#include <utility>
#include <vector>
/// @complexity Time: O(nk), where k is the number of nonzero input terms.
/// Space: O(n + k).
/// @complexity Time: O(log^2 p).
/// Space: O(1).
/// @complexity Time: O(log^3 n) primality testing; Pollard-rho factorization is expected about O(n^(1/4)).
/// Space: O(log n) recursion and factors.
namespace noya {
namespace factorize_internal {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
return u64(u128(a) * b % mod);
}
inline u64 power_mod(u64 a, u64 exp, u64 mod) {
u64 res = 1;
while (exp > 0) {
if (exp & 1) {
res = multiply_mod(res, a, mod);
}
a = multiply_mod(a, a, mod);
exp >>= 1;
}
return res;
}
inline bool miller_rabin(u64 n) {
if (n < 2) {
return false;
}
for (u64 p :
std::array<u64, 12>{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
if (n % p == 0) {
return n == p;
}
}
int shf = __builtin_ctzll(n - 1);
u64 odd = (n - 1) >> shf;
for (u64 bas :
std::array<u64, 7>{2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
if (bas % n == 0) {
continue;
}
u64 val = power_mod(bas % n, odd, n);
if (val == 1 || val == n - 1) {
continue;
}
bool cmp = true;
for (int i = 1; i < shf; i++) {
val = multiply_mod(val, val, n);
if (val == n - 1) {
cmp = false;
break;
}
}
if (cmp) {
return false;
}
}
return true;
}
inline u64 splitmix64(u64 &st) {
u64 z = (st += 0x9e3779b97f4a7c15ULL);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
return z ^ (z >> 31);
}
inline u64 pollard_rho(u64 n) {
if (n % 2 == 0) {
return 2;
}
if (n % 3 == 0) {
return 3;
}
static u64 st = 0x123456789abcdef0ULL;
while (true) {
u64 y = splitmix64(st) % (n - 1) + 1;
u64 c = splitmix64(st) % (n - 1) + 1;
constexpr u64 blk = 128;
u64 g = 1;
u64 r = 1;
u64 q = 1;
u64 x = 0;
u64 sy = 0;
auto nxt = [&](u64 val) {
return u64((u128(multiply_mod(val, val, n)) + c) % n);
};
while (g == 1) {
x = y;
for (u64 i = 0; i < r; i++) {
y = nxt(y);
}
for (u64 off = 0; off < r && g == 1; off += blk) {
sy = y;
for (u64 i = 0; i < std::min(blk, r - off); i++) {
y = nxt(y);
u64 dif = x > y ? x - y : y - x;
q = multiply_mod(q, dif, n);
}
g = std::gcd(q, n);
}
r <<= 1;
}
if (g == n) {
do {
sy = nxt(sy);
u64 dif = x > sy ? x - sy : sy - x;
g = std::gcd(dif, n);
} while (g == 1);
}
if (g != n) {
return g;
}
}
}
inline void collect_factors(u64 n, std::vector<u64> &res) {
if (n == 1) {
return;
}
if (miller_rabin(n)) {
res.push_back(n);
return;
}
u64 fct = pollard_rho(n);
collect_factors(fct, res);
collect_factors(n / fct, res);
}
} // namespace factorize_internal
/// @brief Deterministic Miller-Rabin primality test for unsigned 64-bit
/// integers.
inline bool is_prime(std::uint64_t n) {
return factorize_internal::miller_rabin(n);
}
/// @brief Return the prime factors of n with multiplicity in increasing order.
inline std::vector<std::uint64_t> prime_factors(std::uint64_t n) {
assert(n >= 1);
std::vector<std::uint64_t> res;
factorize_internal::collect_factors(n, res);
std::sort(res.begin(), res.end());
return res;
}
/// @brief Return the prime factorization of n as (prime, exponent) pairs.
inline std::vector<std::pair<std::uint64_t, int>> factorize(std::uint64_t n) {
std::vector<std::pair<std::uint64_t, int>> res;
for (std::uint64_t p : prime_factors(n)) {
if (res.empty() || res.back().first != p) {
res.emplace_back(p, 1);
} else {
res.back().second++;
}
}
return res;
}
} // namespace noya
namespace noya {
/// @brief Compute the smaller square root modulo a prime, or nullopt if no
/// square root exists.
inline std::optional<std::uint64_t> mod_sqrt(std::uint64_t val,
std::uint64_t mod) {
assert(mod >= 2 && is_prime(mod));
val %= mod;
if (mod == 2 || val == 0) {
return val;
}
using factorize_internal::multiply_mod;
using factorize_internal::power_mod;
if (power_mod(val, (mod - 1) / 2, mod) != 1) {
return std::nullopt;
}
if (mod % 4 == 3) {
std::uint64_t rt = power_mod(val, (mod + 1) / 4, mod);
return std::min(rt, mod - rt);
}
std::uint64_t odd = mod - 1;
int exp = 0;
while ((odd & 1) == 0) {
odd >>= 1;
exp++;
}
std::uint64_t nqr = 2;
while (power_mod(nqr, (mod - 1) / 2, mod) != mod - 1) {
nqr++;
}
std::uint64_t rt = power_mod(val, (odd + 1) / 2, mod);
std::uint64_t rem = power_mod(val, odd, mod);
std::uint64_t ste = power_mod(nqr, odd, mod);
int rmn = exp;
while (rem != 1) {
std::uint64_t squ = rem;
int shf = 0;
while (squ != 1 && shf < rmn) {
squ = multiply_mod(squ, squ, mod);
shf++;
}
assert(shf < rmn);
std::uint64_t mul =
power_mod(ste, std::uint64_t(1) << (rmn - shf - 1), mod);
rt = multiply_mod(rt, mul, mod);
ste = multiply_mod(mul, mul, mod);
rem = multiply_mod(rem, ste, mod);
rmn = shf;
}
return std::min(rt, mod - rt);
}
} // namespace noya
namespace noya {
template <class Mint>
using sparse_power_series = std::vector<std::pair<int, Mint>>;
namespace sparse_fps_detail {
template <class Mint> std::vector<Mint> coefficient_inverses(int sz) {
std::vector<Mint> inv(sz);
if (sz > 1) {
inv[1] = Mint(1);
}
for (int vl = 2; vl < sz; vl++) {
inv[vl] = -Mint(Mint::mod() / vl) * inv[Mint::mod() % vl];
}
return inv;
}
template <class Mint>
void check_terms(int sz, const sparse_power_series<Mint> &ter) {
int pre = -1;
for (auto [deg, cf] : ter) {
assert(pre < deg && deg < sz && cf != Mint{});
pre = deg;
}
}
template <class Mint>
std::vector<Mint> unit_power(int sz, const sparse_power_series<Mint> &ter,
Mint exp) {
std::vector<Mint> res(sz);
if (sz == 0)
return res;
res[0] = Mint(1);
std::vector<Mint> inv = coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
vl += cf * res[deg - td] * (exp * Mint(td) - Mint(deg - td));
}
res[deg] = vl * inv[deg];
}
return res;
}
template <class Mint> Mint scalar_power(Mint vl, std::uint64_t exp) {
Mint res = Mint(1);
while (exp > 0) {
if (exp & 1)
res *= vl;
vl *= vl;
exp >>= 1;
}
return res;
}
} // namespace sparse_fps_detail
/// @brief Invert a sparse formal power series. From `f g = 1`, each new
/// coefficient of g is a dot product against the nonconstant terms of f.
template <class Mint>
std::vector<Mint> sparse_fps_inverse(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(!ter.empty() && ter.front().first == 0);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
Mint ic = Mint(1) / ter.front().second;
res[0] = ic;
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
vl += cf * res[deg - td];
}
res[deg] = -vl * ic;
}
return res;
}
/// @brief Compute log(f) for sparse f with f(0)=1. Comparing coefficients in
/// `f (log f)' = f'` gives one linear recurrence per output coefficient.
template <class Mint>
std::vector<Mint> sparse_fps_logarithm(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(!ter.empty() && ter.front().first == 0 &&
ter.front().second == Mint(1));
std::vector<Mint> res(sz);
std::vector<Mint> inv = sparse_fps_detail::coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td == 0)
continue;
if (td > deg)
break;
if (td == deg) {
vl += Mint(deg) * cf;
} else {
vl -= cf * Mint(deg - td) * res[deg - td];
}
}
res[deg] = vl * inv[deg];
}
return res;
}
/// @brief Compute exp(f) for sparse f with f(0)=0. The differential identity
/// `g' = f' g` gives the recurrence directly.
template <class Mint>
std::vector<Mint> sparse_fps_exponential(int sz,
const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
assert(ter.empty() || ter.front().first > 0);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
res[0] = Mint(1);
std::vector<Mint> inv = sparse_fps_detail::coefficient_inverses<Mint>(sz);
for (int deg = 1; deg < sz; deg++) {
Mint vl{};
for (auto [td, cf] : ter) {
if (td > deg)
break;
vl += Mint(td) * cf * res[deg - td];
}
res[deg] = vl * inv[deg];
}
return res;
}
/// @brief Raise a sparse series to a nonnegative integer exponent. The first
/// nonzero monomial determines the output shift and scale; after normalization
/// `f g' = exp f' g` yields an O(nk) recurrence.
template <class Mint>
std::vector<Mint> sparse_fps_power(int sz, const sparse_power_series<Mint> &ter,
std::uint64_t exp) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
std::vector<Mint> res(sz);
if (sz == 0)
return res;
if (exp == 0) {
res[0] = Mint(1);
return res;
}
if (ter.empty())
return res;
int fd = ter.front().first;
if (fd > 0 && exp > std::uint64_t((sz - 1) / fd)) {
return res;
}
int shf = int(std::uint64_t(fd) * exp);
int tar = sz - shf;
Mint ld = ter.front().second;
sparse_power_series<Mint> nrm;
nrm.reserve(ter.size());
for (auto [deg, cf] : ter) {
int nd = deg - fd;
if (nd >= tar)
break;
nrm.emplace_back(nd, cf / ld);
}
std::vector<Mint> uni = sparse_fps_detail::unit_power(tar, nrm, Mint(exp));
Mint scl = sparse_fps_detail::scalar_power(ld, exp);
for (int idx = 0; idx < tar; idx++) {
res[shf + idx] = scl * uni[idx];
}
return res;
}
/// @brief Return one square root of a sparse series when it exists. Its first
/// degree must be even and its leading coefficient a quadratic residue; after
/// removing both, the unit series is raised to the field exponent 1/2.
template <class Mint>
std::optional<std::vector<Mint>>
sparse_fps_square_root(int sz, const sparse_power_series<Mint> &ter) {
assert(sz >= 0);
sparse_fps_detail::check_terms(sz, ter);
if (ter.empty())
return std::vector<Mint>(sz);
int fd = ter.front().first;
if (fd % 2 != 0)
return std::nullopt;
auto rt = mod_sqrt(std::uint64_t(ter.front().second.val()),
std::uint64_t(Mint::mod()));
if (!rt)
return std::nullopt;
int shf = fd / 2;
int tar = sz - fd;
Mint ld = ter.front().second;
sparse_power_series<Mint> nrm;
nrm.reserve(ter.size());
for (auto [deg, cf] : ter) {
int nd = deg - fd;
if (nd >= tar)
break;
nrm.emplace_back(nd, cf / ld);
}
std::vector<Mint> uni =
sparse_fps_detail::unit_power(tar, nrm, Mint(1) / Mint(2));
std::vector<Mint> res(sz);
for (int idx = 0; idx < tar; idx++) {
res[shf + idx] = Mint(*rt) * uni[idx];
}
return res;
}
} // namespace noya