matsutaku-library
This documentation is automatically generated by competitive-verifier/competitive-verifier
test/yosupo/polynomial/log_of_formal_power_series_sparse.test.cpp
- View this file on GitHub
- Last update: 2024-07-11 22:21:57+09:00
- Problem: https://judge.yosupo.jp/problem/log_of_formal_power_series_sparse
Code
#define PROBLEM "https://judge.yosupo.jp/problem/log_of_formal_power_series_sparse"
#include "include/mtl/fps.hpp"
#include <bits/stdc++.h>
using namespace std;
int main() {
int n,k; cin>>n>>k;
vector<int> A(n);
for (int i = 0; i < k; i++) {
int j,a; cin>>j>>a;
A[j] = a;
}
Fps<> f(A.begin(), A.end());
auto g = f.log_sparse();
for (int i = 0; i < n; i++)
cout << g[i] << ' ';
cout << endl;
}#line 1 "test/yosupo/polynomial/log_of_formal_power_series_sparse.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/log_of_formal_power_series_sparse"
#line 2 "include/mtl/bit_manip.hpp"
#include <cstdint>
#include <cassert>
#if __cplusplus >= 202002L
#ifndef MTL_CPP20
#define MTL_CPP20
#endif
#include <bit>
#endif
namespace bm {
/// Count 1s for each 8 bits
inline constexpr uint64_t popcnt_e8(uint64_t x) {
x = (x & 0x5555555555555555) + ((x>>1) & 0x5555555555555555);
x = (x & 0x3333333333333333) + ((x>>2) & 0x3333333333333333);
x = (x & 0x0F0F0F0F0F0F0F0F) + ((x>>4) & 0x0F0F0F0F0F0F0F0F);
return x;
}
/// Count 1s
inline constexpr unsigned popcnt(uint64_t x) {
#ifdef MTL_CPP20
return std::popcount(x);
#else
return (popcnt_e8(x) * 0x0101010101010101) >> 56;
#endif
}
/// Alias to mtl::popcnt(x)
constexpr unsigned popcount(uint64_t x) {
return popcnt(x);
}
/// Count trailing 0s. s.t. *11011000 -> 3
inline constexpr unsigned ctz(uint64_t x) {
#ifdef MTL_CPP20
return std::countr_zero(x);
#else
return popcnt((x & (-x)) - 1);
#endif
}
/// Alias to mtl::ctz(x)
constexpr unsigned countr_zero(uint64_t x) {
return ctz(x);
}
/// Count trailing 1s. s.t. *11011011 -> 2
inline constexpr unsigned cto(uint64_t x) {
#ifdef MTL_CPP20
return std::countr_one(x);
#else
return ctz(~x);
#endif
}
/// Alias to mtl::cto(x)
constexpr unsigned countr_one(uint64_t x) {
return cto(x);
}
inline constexpr unsigned ctz8(uint8_t x) {
return x == 0 ? 8 : popcnt_e8((x & (-x)) - 1);
}
/// [00..0](8bit) -> 0, [**..*](not only 0) -> 1
inline constexpr uint8_t summary(uint64_t x) {
constexpr uint64_t hmask = 0x8080808080808080ull;
constexpr uint64_t lmask = 0x7F7F7F7F7F7F7F7Full;
auto a = x & hmask;
auto b = x & lmask;
b = hmask - b;
b = ~b;
auto c = (a | b) & hmask;
c *= 0x0002040810204081ull;
return uint8_t(c >> 56);
}
/// Extract target area of bits
inline constexpr uint64_t bextr(uint64_t x, unsigned start, unsigned len) {
uint64_t mask = len < 64 ? (1ull<<len)-1 : 0xFFFFFFFFFFFFFFFFull;
return (x >> start) & mask;
}
/// 00101101 -> 00111111 -count_1s-> 6
inline constexpr unsigned log2p1(uint8_t x) {
if (x & 0x80)
return 8;
uint64_t p = uint64_t(x) * 0x0101010101010101ull;
p -= 0x8040201008040201ull;
p = ~p & 0x8080808080808080ull;
p = (p >> 7) * 0x0101010101010101ull;
p >>= 56;
return p;
}
/// 00101100 -mask_mssb-> 00100000 -to_index-> 5
inline constexpr unsigned mssb8(uint8_t x) {
assert(x != 0);
return log2p1(x) - 1;
}
/// 00101100 -mask_lssb-> 00000100 -to_index-> 2
inline constexpr unsigned lssb8(uint8_t x) {
assert(x != 0);
return popcnt_e8((x & -x) - 1);
}
/// Count leading 0s. 00001011... -> 4
inline constexpr unsigned clz(uint64_t x) {
#ifdef MTL_CPP20
return std::countl_zero(x);
#else
if (x == 0)
return 64;
auto i = mssb8(summary(x));
auto j = mssb8(bextr(x, 8 * i, 8));
return 63 - (8 * i + j);
#endif
}
/// Alias to mtl::clz(x)
constexpr unsigned countl_zero(uint64_t x) {
return clz(x);
}
/// Count leading 1s. 11110100... -> 4
inline constexpr unsigned clo(uint64_t x) {
#ifdef MTL_CPP20
return std::countl_one(x);
#else
return clz(~x);
#endif
}
/// Alias to mtl::clo(x)
constexpr unsigned countl_one(uint64_t x) {
return clo(x);
}
inline constexpr unsigned clz8(uint8_t x) {
return x == 0 ? 8 : 7 - mssb8(x);
}
inline constexpr uint64_t bit_reverse(uint64_t x) {
x = ((x & 0x00000000FFFFFFFF) << 32) | ((x & 0xFFFFFFFF00000000) >> 32);
x = ((x & 0x0000FFFF0000FFFF) << 16) | ((x & 0xFFFF0000FFFF0000) >> 16);
x = ((x & 0x00FF00FF00FF00FF) << 8) | ((x & 0xFF00FF00FF00FF00) >> 8);
x = ((x & 0x0F0F0F0F0F0F0F0F) << 4) | ((x & 0xF0F0F0F0F0F0F0F0) >> 4);
x = ((x & 0x3333333333333333) << 2) | ((x & 0xCCCCCCCCCCCCCCCC) >> 2);
x = ((x & 0x5555555555555555) << 1) | ((x & 0xAAAAAAAAAAAAAAAA) >> 1);
return x;
}
/// Check if x is power of 2. 00100000 -> true, 00100001 -> false
constexpr bool has_single_bit(uint64_t x) noexcept {
#ifdef MTL_CPP20
return std::has_single_bit(x);
#else
return x != 0 && (x & (x - 1)) == 0;
#endif
}
/// Bit width needs to represent x. 00110110 -> 6
constexpr int bit_width(uint64_t x) noexcept {
#ifdef MTL_CPP20
return std::bit_width(x);
#else
return 64 - clz(x);
#endif
}
/// Ceil power of 2. 00110110 -> 01000000
constexpr uint64_t bit_ceil(uint64_t x) {
#ifdef MTL_CPP20
return std::bit_ceil(x);
#else
if (x == 0) return 1;
return 1ull << bit_width(x - 1);
#endif
}
/// Floor power of 2. 00110110 -> 00100000
constexpr uint64_t bit_floor(uint64_t x) {
#ifdef MTL_CPP20
return std::bit_floor(x);
#else
if (x == 0) return 0;
return 1ull << (bit_width(x) - 1);
#endif
}
} // namespace bm
#line 3 "include/mtl/modular.hpp"
#include <iostream>
#line 5 "include/mtl/modular.hpp"
template <int MOD>
class Modular {
private:
unsigned int val_;
public:
static constexpr unsigned int mod() { return MOD; }
template<class T>
static constexpr unsigned int safe_mod(T v) {
auto x = (long long)(v%(long long)mod());
if (x < 0) x += mod();
return (unsigned int) x;
}
constexpr Modular() : val_(0) {}
template<class T,
std::enable_if_t<
std::is_integral<T>::value && std::is_unsigned<T>::value
> * = nullptr>
constexpr Modular(T v) : val_(v%mod()) {}
template<class T,
std::enable_if_t<
std::is_integral<T>::value && !std::is_unsigned<T>::value
> * = nullptr>
constexpr Modular(T v) : val_(safe_mod(v)) {}
constexpr unsigned int val() const { return val_; }
constexpr Modular& operator+=(Modular x) {
val_ += x.val();
if (val_ >= mod()) val_ -= mod();
return *this;
}
constexpr Modular operator-() const { return {mod() - val_}; }
constexpr Modular& operator-=(Modular x) {
val_ += mod() - x.val();
if (val_ >= mod()) val_ -= mod();
return *this;
}
constexpr Modular& operator*=(Modular x) {
auto v = (long long) val_ * x.val();
if (v >= mod()) v %= mod();
val_ = v;
return *this;
}
constexpr Modular pow(long long p) const {
assert(p >= 0);
Modular t = 1;
Modular u = *this;
while (p) {
if (p & 1)
t *= u;
u *= u;
p >>= 1;
}
return t;
}
friend constexpr Modular pow(Modular x, long long p) {
return x.pow(p);
}
constexpr Modular inv() const { return pow(mod()-2); }
constexpr Modular& operator/=(Modular x) { return *this *= x.inv(); }
constexpr Modular operator+(Modular x) const { return Modular(*this) += x; }
constexpr Modular operator-(Modular x) const { return Modular(*this) -= x; }
constexpr Modular operator*(Modular x) const { return Modular(*this) *= x; }
constexpr Modular operator/(Modular x) const { return Modular(*this) /= x; }
constexpr Modular& operator++() { return *this += 1; }
constexpr Modular operator++(int) { Modular c = *this; ++(*this); return c; }
constexpr Modular& operator--() { return *this -= 1; }
constexpr Modular operator--(int) { Modular c = *this; --(*this); return c; }
constexpr bool operator==(Modular x) const { return val() == x.val(); }
constexpr bool operator!=(Modular x) const { return val() != x.val(); }
constexpr bool is_square() const {
return pow((mod()-1)/2) == 1;
}
/**
* Return x s.t. x * x = a mod p
* reference: https://zenn.dev/peria/articles/c6afc72b6b003c
*/
constexpr Modular sqrt() const {
if (!is_square())
throw std::runtime_error("not square");
auto mod_eight = mod() % 8;
if (mod_eight == 3 || mod_eight == 7) {
return pow((mod()+1)/4);
} else if (mod_eight == 5) {
auto x = pow((mod()+3)/8);
if (x * x != *this)
x *= Modular(2).pow((mod()-1)/4);
return x;
} else {
Modular d = 2;
while (d.is_square())
d += 1;
auto t = mod()-1;
int s = bm::ctz(t);
t >>= s;
auto a = pow(t);
auto D = d.pow(t);
int m = 0;
Modular dt = 1;
Modular du = D;
for (int i = 0; i < s; i++) {
if ((a*dt).pow(1u<<(s-1-i)) == -1) {
m |= 1u << i;
dt *= du;
}
du *= du;
}
return pow((t+1)/2) * D.pow(m/2);
}
}
friend std::ostream& operator<<(std::ostream& os, const Modular& x) {
return os << x.val();
}
friend std::istream& operator>>(std::istream& is, Modular& x) {
return is >> x.val_;
}
};
using Modular998244353 = Modular<998244353>;
using Modular1000000007 = Modular<(int)1e9+7>;
template<int Id=0>
class DynamicModular {
private:
static unsigned int mod_;
unsigned int val_;
public:
static unsigned int mod() { return mod_; }
static void set_mod(unsigned int m) { mod_ = m; }
template<class T>
static constexpr unsigned int safe_mod(T v) {
auto x = (long long)(v%(long long)mod());
if (x < 0) x += mod();
return (unsigned int) x;
}
constexpr DynamicModular() : val_(0) {}
template<class T,
std::enable_if_t<
std::is_integral<T>::value && std::is_unsigned<T>::value
> * = nullptr>
constexpr DynamicModular(T v) : val_(v%mod()) {}
template<class T,
std::enable_if_t<
std::is_integral<T>::value && !std::is_unsigned<T>::value
> * = nullptr>
constexpr DynamicModular(T v) : val_(safe_mod(v)) {}
constexpr unsigned int val() const { return val_; }
constexpr DynamicModular& operator+=(DynamicModular x) {
val_ += x.val();
if (val_ >= mod()) val_ -= mod();
return *this;
}
constexpr DynamicModular operator-() const { return {mod() - val_}; }
constexpr DynamicModular& operator-=(DynamicModular x) {
val_ += mod() - x.val();
if (val_ >= mod()) val_ -= mod();
return *this;
}
constexpr DynamicModular& operator*=(DynamicModular x) {
auto v = (long long) val_ * x.val();
if (v >= mod()) v %= mod();
val_ = v;
return *this;
}
constexpr DynamicModular pow(long long p) const {
assert(p >= 0);
DynamicModular t = 1;
DynamicModular u = *this;
while (p) {
if (p & 1)
t *= u;
u *= u;
p >>= 1;
}
return t;
}
friend constexpr DynamicModular pow(DynamicModular x, long long p) {
return x.pow(p);
}
// TODO: implement when mod is not prime
constexpr DynamicModular inv() const { return pow(mod()-2); }
constexpr DynamicModular& operator/=(DynamicModular x) { return *this *= x.inv(); }
constexpr DynamicModular operator+(DynamicModular x) const { return DynamicModular(*this) += x; }
constexpr DynamicModular operator-(DynamicModular x) const { return DynamicModular(*this) -= x; }
constexpr DynamicModular operator*(DynamicModular x) const { return DynamicModular(*this) *= x; }
constexpr DynamicModular operator/(DynamicModular x) const { return DynamicModular(*this) /= x; }
constexpr DynamicModular& operator++() { return *this += 1; }
constexpr DynamicModular operator++(int) { DynamicModular c = *this; ++(*this); return c; }
constexpr DynamicModular& operator--() { return *this -= 1; }
constexpr DynamicModular operator--(int) { DynamicModular c = *this; --(*this); return c; }
constexpr bool operator==(DynamicModular x) const { return val() == x.val(); }
constexpr bool operator!=(DynamicModular x) const { return val() != x.val(); }
constexpr bool is_square() const {
return val() == 0 or pow((mod()-1)/2) == 1;
}
/**
* Return x s.t. x * x = a mod p
* reference: https://zenn.dev/peria/articles/c6afc72b6b003c
*/
constexpr DynamicModular sqrt() const {
// assert mod is prime
if (!is_square())
throw std::runtime_error("not square");
if (val() < 2)
return val();
auto mod_eight = mod() % 8;
if (mod_eight == 3 || mod_eight == 7) {
return pow((mod()+1)/4);
} else if (mod_eight == 5) {
auto x = pow((mod()+3)/8);
if (x * x != *this)
x *= DynamicModular(2).pow((mod()-1)/4);
return x;
} else {
DynamicModular d = 2;
while (d.is_square())
++d;
auto t = mod()-1;
int s = bm::ctz(t);
t >>= s;
auto a = pow(t);
auto D = d.pow(t);
int m = 0;
DynamicModular dt = 1;
DynamicModular du = D;
for (int i = 0; i < s; i++) {
if ((a*dt).pow(1u<<(s-1-i)) == -1) {
m |= 1u << i;
dt *= du;
}
du *= du;
}
return pow((t+1)/2) * D.pow(m/2);
}
}
friend std::ostream& operator<<(std::ostream& os, const DynamicModular& x) {
return os << x.val();
}
friend std::istream& operator>>(std::istream& is, DynamicModular& x) {
return is >> x.val_;
}
};
template<int Id>
unsigned int DynamicModular<Id>::mod_;
#include <vector>
template<class ModInt>
struct ModularUtil {
static constexpr int mod = ModInt::mod();
static struct inv_table {
std::vector<ModInt> tb{0,1};
inv_table() : tb({0,1}) {}
} inv_;
void set_inv(int n) {
int m = inv_.tb.size();
if (m > n) return;
inv_.tb.resize(n+1);
for (int i = m; i < n+1; i++)
inv_.tb[i] = -inv_.tb[mod % i] * (mod / i);
}
ModInt& inv(int i) {
set_inv(i);
return inv_.tb[i];
}
};
template<class ModInt>
typename ModularUtil<ModInt>::inv_table ModularUtil<ModInt>::inv_;
#include <array>
namespace math {
constexpr int mod_pow_constexpr(int x, int p, int m) {
long long t = 1;
long long u = x;
while (p) {
if (p & 1) {
t *= u;
t %= m;
}
u *= u;
u %= m;
p >>= 1;
}
return (int) t;
}
constexpr int primitive_root_constexpr(int m) {
if (m == 2) return 1;
if (m == 167772161) return 3;
if (m == 469762049) return 3;
if (m == 754974721) return 11;
if (m == 880803841) return 26;
if (m == 998244353) return 3;
std::array<int, 20> divs{};
int cnt = 0;
int x = m-1;
if (x % 2 == 0) {
divs[cnt++] = 2;
x >>= bm::ctz(x);
}
for (int d = 3; d*d <= x; d += 2) {
if (x % d == 0) {
divs[cnt++] = d;
while (x % d == 0)
x /= d;
}
}
if (x > 1) divs[cnt++] = x;
for (int g = 2; g < m; g++) {
bool ok = true;
for (int i = 0; i < cnt; i++) {
if (mod_pow_constexpr(g, (m-1) / divs[i], m) == 1) {
ok = false;
break;
}
}
if (ok) return g;
}
return -1;
}
template<int m>
constexpr int primitive_root = primitive_root_constexpr(m);
}
#line 6 "include/mtl/ntt.hpp"
namespace math {
constexpr int ceil_pow2(unsigned long long x) {
return x == 0 ? 0 : 64-bm::clz(x-1);
}
}
#ifndef MTL_ARRAY_SET_CONSTEXPR
#if __cplusplus >= 201703L
#define MTL_ARRAY_SET_CONSTEXPR constexpr
#else
#define MTL_ARRAY_SET_CONSTEXPR
#endif
#endif
namespace _ntt {
template<int mod>
struct ntt_info {
using mint = Modular<mod>;
static constexpr int primitive_root() {
return math::primitive_root<mod>;
}
static constexpr int log_n() {
return bm::ctz(mod-1);
}
static MTL_ARRAY_SET_CONSTEXPR
std::array<Modular<mod>, log_n()> coeff(bool forward) {
mint r = mint(primitive_root()).pow((mod-1) >> log_n());
std::array<mint, log_n()> coeff{};
mint iw = forward ? r.inv() : r;
for (int i = log_n()-1; i >= 0; i--) {
coeff[i] = iw;
iw *= iw;
}
return coeff;
}
};
template<class T>
void _iterative_bit_reverse(std::vector<T>& f) {
int log_n = math::ceil_pow2(f.size());
int n = 1 << log_n;
for (int i = 0, j = 1; j < n-1; j++) {
for (int k = n >> 1; k > (i ^= k); k >>= 1);
if (i < j) std::swap(f[i], f[j]);
}
}
template <bool Forward, bool BitReverse, int mod>
void _fft_impl(std::vector<Modular<mod>>& f) {
using mint = Modular<mod>;
int log_n = math::ceil_pow2(f.size());
using info = ntt_info<mod>;
assert(info::log_n() >= log_n);
int n = 1 << log_n;
f.resize(n, 0);
if constexpr (!Forward and BitReverse)
_iterative_bit_reverse(f);
// Cooley-Tukey FFT
static MTL_ARRAY_SET_CONSTEXPR auto coeff = info::coeff(Forward);
if constexpr (Forward) {
for (int log_m = log_n-1; log_m >= 0; log_m -= 2) {
int m = 1<<log_m;
mint w0 = coeff[log_m];
if (log_m == 0) {
for (int chunk = 0; chunk < n; chunk += 2*m) {
mint w = 1;
for (int i = 0; i < m; i++) {
auto p = chunk + i;
auto a = f[p + 0];
auto b = f[p + m];
f[p + 0] = (a + b);
f[p + m] = (a - b) * w;
w *= w0;
}
}
} else { // 4-base
auto w1 = coeff[log_m-1];
auto wh = w0.pow(m/2).val();
for (int chunk = 0; chunk < n; chunk += 2*m) {
mint w = 1, x = 1;
for (int i = 0; i < m/2; i++) {
auto p = chunk + i;
auto ia = p + 0*m/2;
auto ib = p + 1*m/2;
auto ic = p + 2*m/2;
auto id = p + 3*m/2;
long long a = f[ia].val(), b = f[ib].val(), c = f[ic].val(), d = f[id].val();
auto s = a + c, t = b + d, u = (a - c + mod) * w.val() % mod, v = (b - d + mod) * w.val() % mod * wh % mod;
f[ia] = (s + t);
f[ib] = (s - t) * x.val();
f[ic] = (u + v);
f[id] = (u - v) * x.val();
w *= w0;
x *= w1;
}
}
}
}
} else {
int _log_m = 0;
if (log_n % 2 == 1) {
for (int chunk = 0; chunk < n; chunk += 2) {
auto p = chunk;
auto a = f[p + 0];
auto b = f[p + 1];
f[p + 0] = a + b;
f[p + 1] = a - b;
}
_log_m = 1;
}
// 4-base
for (int log_m = _log_m; log_m < log_n; log_m += 2) {
int m = 1<<(log_m+1);
mint w0 = coeff[log_m];
auto w1 = coeff[log_m+1];
auto wh = w1.pow(m/2).val();
for (int chunk = 0; chunk < n; chunk += 2*m) {
mint w = 1, x = 1;
for (int i = 0; i < m/2; i++) {
auto p = chunk + i;
auto ia = p + 0*m/2;
auto ib = p + 1*m/2;
auto ic = p + 2*m/2;
auto id = p + 3*m/2;
long long a = f[ia].val(), b = (long long)f[ib].val() * w.val() % mod,
c = f[ic].val(), d = (long long)f[id].val() * w.val() % mod;
auto s = a + b, t = a - b + mod, u = (c + d) * x.val() % mod, v = (c - d + mod) * x.val() % mod * wh % mod;
f[ia] = s + u;
f[ib] = t + v;
f[ic] = s - u;
f[id] = t - v;
w *= w0;
x *= w1;
}
}
}
}
if constexpr (Forward and BitReverse)
_iterative_bit_reverse(f);
if constexpr (!Forward) {
mint inv = mint(n).inv();
for (int i = 0; i < n; i++) f[i] *= inv;
}
}
template <int mod>
void _fft(std::vector<Modular<mod>>& f) {
_fft_impl<true, true>(f);
}
template <int mod>
void _ifft(std::vector<Modular<mod>>& f) {
_fft_impl<false, true>(f);
}
template <int mod>
void _convolution_fft(std::vector<Modular<mod>>& f) {
_fft_impl<true, false>(f);
}
template <int mod>
void _convolution_ifft(std::vector<Modular<mod>>& f) {
_fft_impl<false, false>(f);
}
} // namespace _ntt
template<class mint>
void ntt_inline(std::vector<mint>& f) {
_ntt::_fft(f);
}
template<class mint>
std::vector<mint> ntt(std::vector<mint> f) {
_ntt::_fft(f);
return f;
}
template<class mint>
void intt_inline(std::vector<mint>& f) {
_ntt::_ifft(f);
}
template<class mint>
std::vector<mint> intt(std::vector<mint> f) {
_ntt::_ifft(f);
return f;
}
template<class mint>
void ntt_convolution_inline(std::vector<mint>& f) {
_ntt::_convolution_fft(f);
}
template<class mint>
std::vector<mint> ntt_convolution(std::vector<mint> f) {
_ntt::_convolution_fft(f);
return f;
}
template<class mint>
void intt_convolution_inline(std::vector<mint>& f) {
_ntt::_convolution_ifft(f);
}
template<class mint>
std::vector<mint> intt_convolution(std::vector<mint> f) {
_ntt::_convolution_ifft(f);
return f;
}
#line 3 "include/mtl/fft.hpp"
#include <complex>
class Fft {
public:
using complex_t = std::complex<double>;
using value_type = complex_t;
private:
size_t n_;
public:
explicit Fft(size_t size) {
n_ = 1;
while (n_ < size)
n_ *= 2;
};
size_t n() const { return n_; }
void fft_inline(std::vector<complex_t>& f) const {
_fft(f);
}
std::vector<complex_t> fft(const std::vector<complex_t>& f) const {
auto ff = f;
ff.resize(n_, 0);
_fft(ff);
return ff;
}
void ifft_inline(std::vector<complex_t>& f) const {
_ifft(f);
}
std::vector<complex_t> ifft(const std::vector<complex_t>& f) const {
auto ff = f;
ff.resize(n_, 0);
_ifft(ff);
return ff;
}
private:
template <bool Forward>
void _fft_impl(std::vector<complex_t>& f) const {
// iterative bit reversal
for (size_t i = 0, j = 1; j < n_-1; j++) {
for (size_t k = n_ >> 1; k > (i ^= k); k >>= 1);
if (i < j) std::swap(f[i], f[j]);
}
// Cooley-Tukey FFT
for (size_t m = 1; m < n_; m *= 2) {
double theta0 = Forward ? M_PI/m : -M_PI/m;
for (size_t chunk = 0; chunk < n_; chunk += 2*m) {
for (size_t i = 0; i < m; i++) {
complex_t w = {cos(theta0*i), -sin(theta0*i)};
auto p = chunk + i;
auto a = f[p + 0];
auto b = f[p + m] * w;
f[p + 0] = a + b;
f[p + m] = a - b;
}
}
}
}
void _fft(std::vector<complex_t>& f) const {
_fft_impl<true>(f);
}
void _ifft(std::vector<complex_t>& f) const {
_fft_impl<false>(f);
for (auto& x : f) {
x /= n_;
}
}
};
#line 4 "include/mtl/convolution.hpp"
const int CONVOLUTION_NAIVE_THRESHOLD = 60;
template<class T>
std::vector<T> convolution_naive(const std::vector<T>& f, const std::vector<T>& g) {
auto n = f.size(), m = g.size();
if (n == 0 or m == 0)
return {};
std::vector<T> h(n+m-1);
if (n < m) {
for (size_t j = 0; j < m; j++)
for (size_t i = 0; i < n; i++)
h[i+j] += f[i] * g[j];
} else {
for (size_t i = 0; i < n; i++)
for (size_t j = 0; j < m; j++)
h[i+j] += f[i] * g[j];
}
return h;
}
template<class T>
std::vector<T> convolution_fft(std::vector<T> f, std::vector<T> g) {
auto n = f.size(), m = g.size();
size_t len = n + m - 1;
size_t l = 1 << math::ceil_pow2(len);
if (n+m-3 < l/2) {
auto fb = f.back(), gb = g.back();
f.pop_back(); g.pop_back();
std::vector<T> h(n+m-1, 0);
h[n+m-2] = fb * gb;
for (size_t i = 0; i < n; i++)
h[i + m-1] += f[i] * gb;
for (size_t i = 0; i < m; i++)
h[i + n-1] += g[i] * fb;
auto fg = convolution(std::move(f), std::move(g));
for (size_t i = 0; i < n+m-3; i++)
h[i] += fg[i];
return h;
}
Fft fft(len);
f.resize(fft.n(), 0);
g.resize(fft.n(), 0);
auto same = f == g;
fft.fft_inline(f);
if (same) {
for (size_t i = 0; i < fft.n(); i++)
f[i] *= f[i];
} else {
fft.fft_inline(g);
for (size_t i = 0; i < fft.n(); i++)
f[i] *= g[i];
}
fft.ifft_inline(f);
std::vector<T> res(len);
for (size_t i = 0; i < len; i++)
res[i] = f[i].real();
return res;
}
template<int Mod>
std::vector<Modular<Mod>> convolution_ntt(std::vector<Modular<Mod>> f, std::vector<Modular<Mod>> g) {
auto n = f.size(), m = g.size();
size_t l = 1 << math::ceil_pow2(n + m - 1);
if (n+m-3 < l/2) {
auto fb = f.back(), gb = g.back();
f.pop_back(); g.pop_back();
std::vector<Modular<Mod>> h(n+m-1, 0);
h[n+m-2] = fb * gb;
for (size_t i = 0; i < n-1; i++)
h[i + m-1] += f[i] * gb;
for (size_t i = 0; i < m-1; i++)
h[i + n-1] += g[i] * fb;
auto fg = convolution(std::move(f), std::move(g));
for (size_t i = 0; i < fg.size(); i++)
h[i] += fg[i];
return h;
}
f.resize(l, 0);
g.resize(l, 0);
bool same = f == g;
ntt_convolution_inline(f);
if (same) {
for (size_t i = 0; i < l; i++)
f[i] *= f[i];
} else {
ntt_convolution_inline(g);
for (size_t i = 0; i < l; i++)
f[i] *= g[i];
}
intt_convolution_inline(f);
return f;
}
template<int Mod>
std::vector<Modular<Mod>> convolution_fft(std::vector<Modular<Mod>> f, std::vector<Modular<Mod>> g) {
return convolution_ntt(std::move(f), std::move(g));
}
std::vector<int> convolution_fft(const std::vector<int>& a, const std::vector<int>& b) {
auto n = a.size(), m = b.size();
std::vector<Modular998244353> f(n), g(m);
for (size_t i = 0; i < n; i++)
f[i] = a[i];
for (size_t i = 0; i < m; i++)
g[i] = b[i];
auto h = convolution_fft(std::move(f), std::move(g));
std::vector<int> res(n+m-1);
for (size_t i = 0; i < n+m-1; i++)
res[i] = h[i].val();
return res;
}
template<class T>
std::vector<T> convolution(const std::vector<T>& f, const std::vector<T>& g) {
auto n = f.size(), m = g.size();
if (n == 0 or m == 0)
return {};
if (std::min(n, m) <= CONVOLUTION_NAIVE_THRESHOLD)
return convolution_naive(f, g);
return convolution_fft(f, g);
}
template<class T>
std::vector<T> convolution(std::vector<T>&& f, std::vector<T>&& g) {
auto n = f.size(), m = g.size();
if (n == 0 or m == 0)
return {};
if (std::min(n, m) <= CONVOLUTION_NAIVE_THRESHOLD)
return convolution_naive(f, g);
return convolution_fft(std::move(f), std::move(g));
}
template<class mint>
std::vector<mint> convolution_garner(const std::vector<mint>& f,
const std::vector<mint>& g) {
auto n = f.size();
auto m = g.size();
if (n == 0 or m == 0)
return {};
constexpr long long nttprimes[] = {754974721, 167772161, 469762049};
using mint0 = Modular<754974721>;
using mint1 = Modular<167772161>;
using mint2 = Modular<469762049>;
std::vector<mint0> a0(n), b0(m);
std::vector<mint1> a1(n), b1(m);
std::vector<mint2> a2(n), b2(m);
for(size_t i = 0; i < n; i++) a0[i] = f[i].val(), a1[i] = f[i].val(), a2[i] = f[i].val();
for(size_t i = 0; i < m; i++) b0[i] = g[i].val(), b1[i] = g[i].val(), b2[i] = g[i].val();
auto c0 = convolution_ntt(a0, b0);
auto c1 = convolution_ntt(a1, b1);
auto c2 = convolution_ntt(a2, b2);
constexpr long long m01 = (long long) nttprimes[0] * nttprimes[1];
constexpr long long m0_inv_m1 = mint1(nttprimes[0]).inv().val();
constexpr long long m01_inv_m2 = mint2(m01).inv().val();
const int mod = mint::mod();
auto garner = [&](mint0 x0, mint1 x1, mint2 x2) -> mint {
int r0 = x0.val(), r1 = x1.val(), r2 = x2.val();
int v1 = (m0_inv_m1 * (r1 + nttprimes[1] - r0)) % nttprimes[1];
auto v2 = (mint2(r2) - r0 - mint2(nttprimes[0]) * v1) * mint2(m01_inv_m2);
return mint(r0 + 1LL * nttprimes[0] * v1 + m01 % mod * v2.val());
};
std::vector<mint> c(c0.size());
for(size_t i = 0; i < c.size(); i++) c[i] = garner(c0[i], c1[i], c2[i]);
return c;
}
template<int Mod>
std::vector<Modular<Mod>> convolution_large(const std::vector<Modular<Mod>>& f,
const std::vector<Modular<Mod>>& g) {
auto n = f.size();
auto m = g.size();
if (n == 0 or m == 0)
return {};
constexpr auto h = bm::ctz(Mod-1);
if (n+m-1 < 1u<<h)
return convolution(f ,g);
size_t l = h-1;
int cs = ((n-1)>>l)+1;
int ds = ((m-1)>>l)+1;
using result_type = std::vector<Modular<Mod>>;
std::vector<result_type> c(cs, result_type(2<<l));
std::vector<result_type> d(ds, result_type(2<<l));
size_t mask = (1<<l)-1;
for (size_t i = 0; i < n; i++)
c[i>>l][i&mask] = f[i];
for (size_t i = 0; i < m; i++)
d[i>>l][i&mask] = g[i];
for (auto& v: c) ntt_convolution_inline(v);
for (auto& v: d) ntt_convolution_inline(v);
result_type res(n+m-1);
result_type e(2<<l);
for (int i = 0; i < cs+ds-1; i++) {
e.assign(2<<l, 0);
for (int j = std::max(i-ds+1, 0); j < std::min(cs, i+1); j++) {
int k = i-j;
assert(k >= 0 and k < ds);
for (int x = 0; x < 2<<l; x++)
e[x] += c[j][x] * d[k][x];
}
intt_convolution_inline(e);
auto len = std::min((size_t)2<<l, n+m-1-(i<<l));
for (size_t x = 0; x < len; x++)
res[(i<<l) + x] += e[x];
}
return res;
}
size_t complexity_of_convolution(size_t fsize, size_t gsize) {
int lg = math::ceil_pow2(fsize + gsize - 1);
return lg << lg;
}
#line 5 "include/mtl/fps.hpp"
#include <initializer_list>
#line 7 "include/mtl/fps.hpp"
#include <algorithm>
template<int Mod>
class SparseFps;
template<int Mod = 998244353>
class Fps : public std::vector<Modular<Mod>> {
using _base = std::vector<Modular<Mod>>;
public:
using mint = Modular<Mod>;
using f = std::vector<mint>;
explicit Fps(size_t n=0, mint v=0) : f(n, v) {}
template<typename It> Fps(It begin, It end) : f(begin, end) {}
Fps(std::initializer_list<mint> l) : f(l) {}
explicit Fps(const std::vector<mint>& l) : f(l) {}
explicit Fps(std::vector<mint>&& l) : f(std::move(l)) {}
Fps& normalize() {
int n = (int) _base::size();
while (n-1 >= 0 and (*this)[n-1] == 0) n--;
_base::resize(n);
return *this;
}
Fps& inline_pre(size_t m) {
_base::resize(m, 0);
return *this;
}
/// f mod x^m
Fps pre(size_t m) const {
if (m <= _base::size())
return Fps(_base::begin(), _base::begin()+m);
else
return Fps(*this).inline_pre(m);
}
bool operator==(const Fps& rhs) const {
for (size_t i = 0, n = std::max(_base::size(), rhs.size()); i < n; i++) {
auto l = i < _base::size() ? (*this)[i] : 0;
auto r = i < rhs.size() ? rhs[i] : 0;
if (l != r) return false;
}
return true;
}
bool operator!=(const Fps& rhs) const {
return !(*this == rhs);
}
Fps& operator+=(mint x) {
if (_base::empty())
_base::resize(1);
(*this)[0] += x;
return *this;
}
Fps& operator-=(mint x) { return *this += -x; }
Fps& operator*=(mint x) {
for (auto& v : *this) v *= x;
return *this;
}
Fps& operator/=(mint x) { return *this *= x.inv(); }
Fps operator+(mint x) const { return Fps(*this) += x; }
Fps operator-(mint x) const { return Fps(*this) -= x; }
Fps operator*(mint x) const { return Fps(*this) *= x; }
Fps operator/(mint x) const { return Fps(*this) /= x; }
Fps operator-() const {
auto g = *this;
for (auto& v : g) v = -v;
return g;
}
friend Fps operator+(mint x, const Fps& y) { return y+x; }
friend Fps operator-(mint x, const Fps& y) { return -y+x; }
friend Fps operator*(mint x, const Fps& y) { return y*x; }
friend Fps operator/(mint x, const Fps& y) { return y.inv()*x; }
Fps& operator+=(const Fps& r) {
if (r.size() > _base::size())
_base::resize(r.size());
for (size_t i = 0; i < r.size(); i++)
(*this)[i] += r[i];
return *this;
}
Fps& operator-=(const Fps& r) {
if (r.size() > _base::size())
_base::resize(r.size());
for (size_t i = 0; i < r.size(); i++)
(*this)[i] -= r[i];
return *this;
}
Fps& dot(const Fps& r) {
if (_base::size() > r.size())
_base::resize(r.size());
for (size_t i = 0; i < _base::size(); i++)
(*this)[i] *= r[i];
return *this;
}
size_t count_terms(int n = -1) const {
if (n == -1) n = (int)_base::size();
n = std::min(n, (int)_base::size());
return std::count_if(_base::begin(), _base::begin()+n, [](mint x) { return x != 0; });
}
using sparse_type = std::vector<std::pair<size_t, mint>>;
sparse_type term_ties() const {
return term_ties(0, _base::size());
}
sparse_type term_ties(size_t n) const {
return term_ties(0, n);
}
sparse_type term_ties(size_t front, size_t back) const {
sparse_type ret;
for (size_t i = front; i < back; i++)
if ((*this)[i] != 0)
ret.emplace_back(i, (*this)[i]);
return ret;
}
private:
template<class F>
Fps& _mul_set_dense(F&& r) {
return *this = Fps(convolution_ntt(std::move(*this), std::forward<F>(r)));
}
/**
* Complexity: O(NK)
* where N is count of non-zero terms of self and K is count of non-zero terms of r
*/
Fps& _mul_set_sparse(const Fps& r) {
auto ri = r.term_ties();
if (ri.empty()) return *this = Fps();
Fps ret(_base::size() + ri.back().first);
for (size_t i = 0; i < _base::size(); i++) if ((*this)[i] != 0) for (auto j:ri) {
ret[i + j.first] += (*this)[i] * j.second;
}
return *this = ret;
}
template<class F>
Fps& _mul_set(F&& r) {
return
r.count_terms() < 100 ?
_mul_set_sparse(std::forward<F>(r)) :
_mul_set_dense(std::forward<F>(r));
}
public:
Fps& operator*=(const Fps& r) {
return _mul_set(r);
}
Fps& operator*=(Fps&& r) {
return _mul_set(std::move(r));
}
Fps inv_dense(int n = -1) const {
assert(!_base::empty() and (*this)[0] != 0);
if (n == -1) n = (int) _base::size();
if (n == 0) return Fps();
// Newton descent
// find g, s.t. F(g) = a
// g_{n+1} = g_n - (F(g_n) - a) / F'(g_n)
// find g, s.t. g^{-1} = f
// g_{n+1} = g_n - (g_n^{-1} - f) / -g_n^{-2}
// = 2g_n - g_n^2 f
// g_0 = f_0^{-1}
Fps g,fm,fgg;
g.reserve(n); fm.reserve(n); fgg.reserve(n*2-1);
g.push_back((*this)[0].inv());
fm.push_back((*this)[0]);
for (int m = 1; m < n; m <<= 1) {
int nm = std::min(m*2, n);
fm.resize(nm);
for (int i = m; i < std::min(nm, (int)_base::size()); i++)
fm[i] = (*this)[i];
fgg = g;
fgg *= g;
fgg *= fm;
fgg.inline_pre(nm);
g.resize(nm);
for (int i = m; i < nm; i++)
g[i] = -fgg[i];
}
// assert((pre(n)*g).pre(n) == Fps{1});
return g;
}
Fps inv_sparse(int n = -1) const;
Fps inv(int n = -1) const {
if (n == -1) n = (int) _base::size();
return count_terms(n) < 100 ? inv_sparse(n) : inv_dense(n);
}
Fps& operator/=(const Fps& r) {
return *this *= r.inv();
}
// f x^m = sum f_i x^{i+m}
Fps operator<<(size_t m) const {
Fps ret(m + _base::size());
std::copy(_base::begin(), _base::end(), ret.begin()+m);
return ret;
}
// sum f_{i+m} x^i
Fps operator>>(size_t m) const {
if (m >= _base::size()) return Fps();
return Fps(_base::begin()+m, _base::end());
}
Fps& operator<<=(size_t m) {
size_t s = _base::size();
_base::resize(m + s);
for (size_t i = 0; i < s; i++)
(*this)[m+s-1-i] = (*this)[s-1-i];
std::fill(_base::begin(), _base::begin()+m, 0);
return *this;
}
Fps& operator>>=(size_t m) {
if (m >= _base::size()) return *this = Fps();
_base::erase(_base::begin(), _base::begin()+m);
return *this;
}
Fps& inline_diff() {
if (_base::empty()) return *this;
for (size_t i = 1; i < _base::size(); i++)
(*this)[i-1] = (*this)[i] * i;
_base::pop_back();
return *this;
}
Fps diff() const {
return Fps(*this).inline_diff();
}
Fps& inline_inte() {
if (_base::empty()) return *this;
_base::push_back(0);
ModularUtil<mint> mu;
mu.set_inv(_base::size());
for (int i = (int) _base::size()-1; i > 0; i--)
(*this)[i] = (*this)[i-1] * mu.inv(i);
(*this)[0] = 0;
return *this;
}
Fps inte() const {
return Fps(*this).inline_inte();
}
Fps log_dense(int n = -1) const {
assert(!_base::empty() and _base::operator[](0) == 1);
if (n == -1) n = (int) _base::size();
// integral(f' / f)
return (diff() * inv_dense(n-1)).inline_pre(n-1).inline_inte().inline_pre(n);
}
Fps log_sparse(int n = -1) const {
assert(!_base::empty() and _base::operator[](0) == 1);
if (n == -1) n = (int) _base::size();
// integral(f' / f)
return (diff() / SparseFps<Mod>(*this, n-1)).inline_pre(n-1).inline_inte().inline_pre(n);
}
/**
* define log (1-f) = -sum_n f^n / n
* satisfy log(f)' = f'/f, log(fg) = log f + log g
*/
Fps log(int n = -1) const {
assert(!_base::empty() and _base::operator[](0) == 1);
if (n == -1) n = (int) _base::size();
// integral(f' / f)
return (diff() * inv(n-1)).inline_pre(n-1).inline_inte().inline_pre(n);
}
Fps exp_dense(int n = -1) const;
Fps exp_sparse(int n = -1) const;
/**
* define exp f = sum_n f^n / n!
* satisfy (exp f)' = (exp f)f', exp(f + g) = exp(f)exp(g)
*/
Fps exp(int n = -1) const {
assert(_base::empty() or _base::operator[](0) == 0);
return count_terms() < 100 ? exp_sparse(n) : exp_dense(n);
}
private:
Fps _pow_1_dense(long long n) const;
Fps _pow_1_sparse(long long n) const;
Fps _pow_1(long long n) const {
assert(!_base::empty() and _base::operator[](0) == 1);
return count_terms() < 100 ? _pow_1_sparse(n) : _pow_1_dense(n);
}
public:
/**
* f^n = exp(n log f)
* f^{ab} = (f^a)^b, f^{a+b} = f^a f^b
*/
Fps pow(long long n) const;
private:
/**
* define f^{1/2}' = 1/2f^{-1/2} if [x^0]f = 1
* F = f^{1/2}
* fF' = 1/2 F f'
*/
Fps _sqrt_1() const {
assert(!_base::empty() and (*this)[0] == 1);
return _pow_1(mint(2).inv().val());
}
public:
Fps sqrt() const {
// f = sum_i a_i x^i = c_k x^k sum_i a_i/c_k x^{i-k} = c_k x^k g
// f^{1/2} = c_k^{1/2} x^{k/2} g^{1/2}
// s.t. c_k is square and k is even
size_t k = 0;
while (k < _base::size() and (*this)[k] == 0) k++;
if (k == _base::size()) return Fps();
if (k%2==1)
throw std::runtime_error("minimum degree must be even");
auto ck = (*this)[k];
if (!ck.is_square())
throw std::runtime_error("not square");
auto g = *this >> k;
g /= ck;
g = g._sqrt_1();
g *= ck.sqrt();
g <<= k/2;
return g;
}
bool is_square() const {
size_t k = 0;
while (k < _base::size() and (*this)[k] == 0) k++;
if (k == _base::size()) return true;
if (k%2==1)
return false;
if (!(*this)[k].is_square())
return false;
return true;
}
private:
template<class F>
Fps& _mod_set(F&& r) {
normalize();
Fps q = *this / r;
return *this -= q * std::forward<F>(r);
}
public:
Fps& operator%=(const Fps& r) {
return _mod_set(r);
}
Fps& operator%=(Fps&& r) {
return _mod_set(std::move(r));
}
Fps operator+(const Fps& r) const { return Fps(*this) += r; }
Fps operator+(Fps&& r) const { return Fps(*this) += std::move(r); }
Fps operator-(const Fps& r) const { return Fps(*this) -= r; }
Fps operator-(Fps&& r) const { return Fps(*this) -= std::move(r); }
Fps operator*(const Fps& r) const { return Fps(*this) *= r; }
Fps operator*(Fps&& r) const { return Fps(*this) *= std::move(r); }
Fps operator/(const Fps& r) const { return Fps(*this) /= r; }
Fps operator/(Fps&& r) const { return Fps(*this) /= std::move(r); }
Fps operator%(const Fps& r) const { return Fps(*this) %= r; }
Fps operator%(Fps&& r) const { return Fps(*this) %= std::move(r); }
};
#line 4 "include/mtl/sparse_fps.hpp"
#include <cstddef>
#line 6 "include/mtl/sparse_fps.hpp"
#include <utility>
#line 9 "include/mtl/sparse_fps.hpp"
template<int Mod = 998244353>
class SparseFps : public std::vector<std::pair<size_t, Modular<Mod>>> {
using _base = std::vector<std::pair<size_t, Modular<Mod>>>;
public:
using mint = Modular<Mod>;
using value_type = typename _base::value_type;
using fps_type = Fps<Mod>;
SparseFps() = default;
template<class InputIt>
SparseFps(InputIt first, InputIt last) {
for (auto it = first; it != last; ++it) {
if (it->second != mint(0)) {
this->emplace_back(*it);
}
}
}
SparseFps(std::initializer_list<value_type> init) : SparseFps(init.begin(), init.end()) {}
SparseFps(std::initializer_list<mint> init) {
for (auto it = init.begin(); it != init.end(); ++it) {
if (*it != mint(0)) {
auto k = it - init.begin();
this->emplace_back(k, *it);
}
}
}
SparseFps(const fps_type& rhs, size_t n = Mod) {
n = std::min(n, rhs.size());
for (size_t i = 0; i < n; ++i) {
if (rhs[i] != mint(0)) {
this->emplace_back(i, rhs[i]);
}
}
}
size_t deg() const {
return this->empty() ? 0 : this->back().first+1;
}
SparseFps& inline_pre(int n) {
for (size_t i = 0; i < this->size(); i++) {
if ((*this)[i].first >= n) {
this->resize(i);
break;
}
}
}
SparseFps pre(int n) const {
auto it = std::lower_bound(this->begin(), this->end(), std::make_pair(n, mint(0)));
return SparseFps(this->begin(), it);
}
template<class Value>
void push(Value&& x) {
if (x.second == mint(0)) return;
this->emplace_back(std::forward<Value>(x));
}
void emplace(size_t deg, mint val) {
push(std::make_pair(deg, val));
}
fps_type to_fps() const {
if (this->empty()) return fps_type();
fps_type res(deg(), 0);
for (auto& x : *this) {
res[x.first] = x.second;
}
return res;
}
operator fps_type() const {
return to_fps();
}
bool operator==(const SparseFps& rhs) const {
auto it1 = this->begin();
auto it2 = rhs.begin();
while (it1 != this->end() and it2 != rhs.end()) {
while (it1 != this->end() and it1->second == mint(0)) ++it1;
while (it2 != rhs.end() and it2->second == mint(0)) ++it2;
if (it1 == this->end() or it2 == rhs.end()) break;
if (*it1 != *it2) return false;
++it1; ++it2;
}
return it1 == this->end() and it2 == rhs.end();
}
SparseFps operator+(const SparseFps& rhs) const {
// degrees must be sorted
SparseFps res;
auto it1 = this->begin(), it2 = rhs.begin();
while (it1 != this->end() && it2 != rhs.end()) {
if (it1->first == it2->first) {
res.push_back(*it1 + *it2);
++it1; ++it2;
} else if (it1->first < it2->first) {
res.push_back(*it1);
++it1;
} else {
res.push_back(*it2);
++it2;
}
}
while (it1 != this->end()) {
res.push_back(*it1);
++it1;
}
while (it2 != rhs.end()) {
res.push_back(*it2);
++it2;
}
return res;
}
SparseFps operator-(const SparseFps& rhs) const {
// degrees must be sorted
SparseFps res;
auto it1 = this->begin(), it2 = rhs.begin();
while (it1 != this->end() && it2 != rhs.end()) {
if (it1->first == it2->first) {
res.push_back(*it1 - *it2);
++it1; ++it2;
} else if (it1->first < it2->first) {
res.push_back(*it1);
++it1;
} else {
res.push_back(-*it2);
++it2;
}
}
while (it1 != this->end()) {
res.push_back(*it1);
++it1;
}
while (it2 != rhs.end()) {
res.push_back(-*it2);
++it2;
}
return res;
}
SparseFps operator-() const {
SparseFps res;
for (auto& x : *this) {
res.push_back(x.first, -x.second);
}
return res;
}
SparseFps operator*(const SparseFps& rhs) const {
if (this->empty() or rhs.empty()) return SparseFps();
SparseFps res;
for (auto& x : *this) {
for (auto& y : rhs) {
res.push_back(x.first + y.first, x.second * y.second);
}
}
return res;
}
fps_type operator*(const fps_type& rhs) const {
if (this->empty() or rhs.empty()) return fps_type{};
fps_type res(deg() - 1 + rhs.size(), 0);
for (auto& x : *this) {
for (size_t i = 0; i < rhs.size(); ++i) {
res[i+x.first] += x.second * rhs[i];
}
}
return res;
}
friend fps_type operator*(const fps_type& lhs, const SparseFps& rhs) {
return rhs * lhs;
}
fps_type inv(int n) const {
assert((!this->empty() and this->front().first == 0 and this->front().second != mint(0)));
if (n == 0) return fps_type();
// fg = 1 => (n>0) sum_i f_i g_{n-i} = 0
// f_0 g_n = - sum_{i=1}^n f_i g_{n-i}
fps_type g(n);
auto ifz = this->front().second.inv();
g[0] = ifz;
for (size_t i = 1; i < (size_t) n; i++) {
mint s = 0;
for (size_t j = 1; j < this->size(); j++) {
auto& x = this->operator[](j);
if (x.first > i) break;
s += x.second * g[i-x.first];
}
g[i] = -s * ifz;
}
// assert((pre(n)*g).pre(n) == Fps{1});
return g;
}
friend fps_type operator/(const fps_type& lhs, const SparseFps& rhs) {
if (lhs == fps_type()) return fps_type();
assert(!rhs.empty() and rhs[0].first == 0 and rhs[0].second != mint(0));
if (rhs[0].second == mint(1)) {
return div_one(lhs, rhs);
} else {
auto icoeff = rhs[0].second.inv();
return div_one(lhs, rhs * icoeff) * icoeff;
}
}
private:
friend fps_type div_one(const fps_type& lhs, const SparseFps& rhs) {
assert(lhs != fps_type());
assert(!rhs.empty() and rhs[0].second == mint(1));
fps_type ret(lhs.size() + rhs.deg() - 1);
std::copy(lhs.begin(), lhs.end(), ret.begin());
for (size_t i = 1; i < ret.size(); i++) {
for (size_t j = 1; j < rhs.size(); j++) {
auto& x = rhs[j];
if (i < x.first) break;
ret[i] += ret[i-x.first] * -x.second;
}
}
return ret;
}
public:
SparseFps diff(int n = Mod) const {
SparseFps res;
for (auto& x : *this) {
if (x.first == 0) continue;
if (x.first >= n) break;
res.push_back(x.first - 1, x.second * x.first);
}
return res;
}
SparseFps& inline_diff() {
for (auto& x : *this) {
x.second *= x.first;
--x.first;
}
return *this;
}
SparseFps inte(int n = Mod) const {
SparseFps res;
for (auto& x : *this) {
if (x.first + 1 >= n) break;
res.push_back(x.first + 1, x.second / (x.first + 1));
}
return res;
}
SparseFps inline_inte() const {
for (auto& x : *this) {
x.second /= x.first + 1;
++x.first;
}
return *this;
}
fps_type exp(int deg) const;
fps_type log(int n = Mod) const {
assert(!_base::empty() and _base::operator[](0) == 1);
// integral(f' / f)
return (diff(n).to_fps().inline_pre(n-1) / pre(n-1)).inline_pre(n-1).inline_inte().inline_pre(n);
}
fps_type pow(long long x, int n) const;
SparseFps operator+(mint x) const {
if (!this->empty() and this->front().first == 0) {
auto res = *this;
res[0].second += x;
return res;
} else {
SparseFps res{x};
res.insert(res.end(), this->begin(), this->end());
return res;
}
}
SparseFps operator-(mint x) const {
return *this + (-x);
}
SparseFps operator*(mint x) const {
assert(x != mint(0));
auto res = *this;
for (auto& y : res) {
y.second *= x;
}
return res;
}
SparseFps operator/(mint x) const {
assert(x != mint(0));
return *this * x.inv();
}
SparseFps& operator*=(mint x) {
assert(x != mint(0));
for (auto& y : *this) {
y.second *= x;
}
return *this;
}
SparseFps& operator/=(mint x) {
assert(x != mint(0));
return *this *= x.inv();
}
};
#line 347 "include/mtl/fps.hpp"
template<int Mod>
Fps<Mod> Fps<Mod>::inv_sparse(int n) const {
assert(!_base::empty() and (*this)[0] != 0);
if (n == -1) n = (int) _base::size();
if (n == 0) return Fps();
return SparseFps<Mod>(*this, n).inv(n);
}
#line 3 "test/yosupo/polynomial/log_of_formal_power_series_sparse.test.cpp"
#include <bits/stdc++.h>
using namespace std;
int main() {
int n,k; cin>>n>>k;
vector<int> A(n);
for (int i = 0; i < k; i++) {
int j,a; cin>>j>>a;
A[j] = a;
}
Fps<> f(A.begin(), A.end());
auto g = f.log_sparse();
for (int i = 0; i < n; i++)
cout << g[i] << ' ';
cout << endl;
}
Test cases
| Env | Name | Status | Elapsed | Memory |
|---|---|---|---|---|
| g++ | example_00 |
AC |
6 ms | 3 MB |
| g++ | example_01 |
AC |
5 ms | 3 MB |
| g++ | max_random_00 |
AC |
113 ms | 30 MB |
| g++ | max_random_01 |
AC |
113 ms | 30 MB |
| g++ | max_random_02 |
AC |
110 ms | 30 MB |
| g++ | max_random_03 |
AC |
109 ms | 30 MB |
| g++ | max_random_04 |
AC |
107 ms | 30 MB |
| g++ | min_K_00 |
AC |
33 ms | 12 MB |
| g++ | min_K_01 |
AC |
39 ms | 15 MB |
| g++ | random_00 |
AC |
37 ms | 12 MB |
| g++ | random_01 |
AC |
42 ms | 16 MB |
| g++ | random_02 |
AC |
54 ms | 19 MB |
| g++ | random_03 |
AC |
37 ms | 14 MB |
| g++ | random_04 |
AC |
65 ms | 25 MB |
| g++ | small_N_00 |
AC |
6 ms | 3 MB |
| g++ | small_N_01 |
AC |
5 ms | 3 MB |
| g++ | small_N_02 |
AC |
5 ms | 3 MB |
| g++ | small_N_03 |
AC |
5 ms | 3 MB |
| g++ | small_N_04 |
AC |
5 ms | 3 MB |
| g++ | small_dense_00 |
AC |
103 ms | 27 MB |
| g++ | small_dense_01 |
AC |
124 ms | 27 MB |
| g++ | small_dense_02 |
AC |
107 ms | 27 MB |
| g++ | small_dense_03 |
AC |
74 ms | 28 MB |
| g++ | small_dense_04 |
AC |
110 ms | 27 MB |
| clang++ | example_00 |
AC |
6 ms | 3 MB |
| clang++ | example_01 |
AC |
5 ms | 3 MB |
| clang++ | max_random_00 |
AC |
95 ms | 30 MB |
| clang++ | max_random_01 |
AC |
96 ms | 30 MB |
| clang++ | max_random_02 |
AC |
96 ms | 31 MB |
| clang++ | max_random_03 |
AC |
95 ms | 30 MB |
| clang++ | max_random_04 |
AC |
95 ms | 30 MB |
| clang++ | min_K_00 |
AC |
31 ms | 12 MB |
| clang++ | min_K_01 |
AC |
35 ms | 15 MB |
| clang++ | random_00 |
AC |
31 ms | 14 MB |
| clang++ | random_01 |
AC |
40 ms | 16 MB |
| clang++ | random_02 |
AC |
51 ms | 19 MB |
| clang++ | random_03 |
AC |
35 ms | 14 MB |
| clang++ | random_04 |
AC |
58 ms | 24 MB |
| clang++ | small_N_00 |
AC |
6 ms | 3 MB |
| clang++ | small_N_01 |
AC |
5 ms | 3 MB |
| clang++ | small_N_02 |
AC |
5 ms | 3 MB |
| clang++ | small_N_03 |
AC |
5 ms | 3 MB |
| clang++ | small_N_04 |
AC |
5 ms | 3 MB |
| clang++ | small_dense_00 |
AC |
112 ms | 27 MB |
| clang++ | small_dense_01 |
AC |
118 ms | 27 MB |
| clang++ | small_dense_02 |
AC |
99 ms | 27 MB |
| clang++ | small_dense_03 |
AC |
61 ms | 28 MB |
| clang++ | small_dense_04 |
AC |
104 ms | 27 MB |