This documentation is automatically generated by competitive-verifier/competitive-verifier
#include "cp-algo/math/fft.hpp"
#ifndef CP_ALGO_MATH_FFT_HPP
#define CP_ALGO_MATH_FFT_HPP
#include "../number_theory/modint.hpp"
#include "cvector.hpp"
namespace cp_algo::math::fft {
template<modint_type base>
struct dft<base> {
int split;
cvector A, B;
dft(auto const& a, size_t n): A(n), B(n) {
split = int(std::sqrt(base::mod())) + 1;
cvector::exec_on_roots(2 * n, size(a), [&](size_t i, point rt) {
size_t ti = std::min(i, i - n);
auto rem = std::remainder(a[i].rem(), split);
auto quo = (ftype(a[i].rem()) - rem) / split;
A.set(ti, A.get(ti) + rem * rt);
B.set(ti, B.get(ti) + quo * rt);
});
if(n) {
A.fft();
B.fft();
}
}
void mul(auto &&C, auto const& D, auto &res, size_t k) {
assert(A.size() == C.size());
size_t n = A.size();
if(!n) {
res = {};
return;
}
for(size_t i = 0; i < n; i += flen) {
auto tmp = A.vget(i) * D.vget(i) + B.vget(i) * C.vget(i);
A.set(i, A.vget(i) * C.vget(i));
B.set(i, B.vget(i) * D.vget(i));
C.set(i, tmp);
}
A.ifft();
B.ifft();
C.ifft();
auto splitsplit = (base(split) * split).rem();
cvector::exec_on_roots(2 * n, std::min(n, k), [&](size_t i, point rt) {
rt = conj(rt);
auto Ai = A.get(i) * rt;
auto Bi = B.get(i) * rt;
auto Ci = C.get(i) * rt;
int64_t A0 = llround(real(Ai));
int64_t A1 = llround(real(Ci));
int64_t A2 = llround(real(Bi));
res[i] = A0 + A1 * split + A2 * splitsplit;
if(n + i >= k) {
return;
}
int64_t B0 = llround(imag(Ai));
int64_t B1 = llround(imag(Ci));
int64_t B2 = llround(imag(Bi));
res[n + i] = B0 + B1 * split + B2 * splitsplit;
});
}
void mul_inplace(auto &&B, auto& res, size_t k) {
mul(B.A, B.B, res, k);
}
void mul(auto const& B, auto& res, size_t k) {
mul(cvector(B.A), B.B, res, k);
}
std::vector<base> operator *= (dft &B) {
std::vector<base> res(2 * A.size());
mul_inplace(B, res, size(res));
return res;
}
std::vector<base> operator *= (dft const& B) {
std::vector<base> res(2 * A.size());
mul(B, res, size(res));
return res;
}
auto operator * (dft const& B) const {
return dft(*this) *= B;
}
point operator [](int i) const {return A.get(i);}
};
void mul_slow(auto &a, auto const& b, size_t k) {
if(empty(a) || empty(b)) {
a.clear();
} else {
size_t n = std::min(k, size(a));
size_t m = std::min(k, size(b));
a.resize(k);
for(int j = int(k - 1); j >= 0; j--) {
a[j] *= b[0];
for(int i = std::max(j - (int)n, 0) + 1; i < std::min(j + 1, (int)m); i++) {
a[j] += a[j - i] * b[i];
}
}
}
}
size_t com_size(size_t as, size_t bs) {
if(!as || !bs) {
return 0;
}
return std::max(flen, std::bit_ceil(as + bs - 1) / 2);
}
void mul_truncate(auto &a, auto const& b, size_t k) {
using base = std::decay_t<decltype(a[0])>;
if(std::min({k, size(a), size(b)}) < magic) {
mul_slow(a, b, k);
return;
}
auto n = std::max(flen, std::bit_ceil(
std::min(k, size(a)) + std::min(k, size(b)) - 1
) / 2);
a.resize(k);
auto A = dft<base>(a, n);
if(&a == &b) {
A.mul(A, a, k);
} else {
A.mul_inplace(dft<base>(b | std::views::take(k), n), a, k);
}
}
void mul(auto &a, auto const& b) {
if(size(a)) {
mul_truncate(a, b, size(a) + size(b) - 1);
}
}
}
#endif // CP_ALGO_MATH_FFT_HPP
#line 1 "cp-algo/math/fft.hpp"
#line 1 "cp-algo/number_theory/modint.hpp"
#line 1 "cp-algo/math/common.hpp"
#include <functional>
#include <cstdint>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
const int maxn = CP_ALGO_MAXN;
#else
const int maxn = 1 << 19;
#endif
const int magic = 64; // threshold for sizes to run the naive algo
auto bpow(auto const& x, auto n, auto const& one, auto op) {
if(n == 0) {
return one;
} else {
auto t = bpow(x, n / 2, one, op);
t = op(t, t);
if(n % 2) {
t = op(t, x);
}
return t;
}
}
auto bpow(auto x, auto n, auto ans) {
return bpow(x, n, ans, std::multiplies{});
}
template<typename T>
T bpow(T const& x, auto n) {
return bpow(x, n, T(1));
}
}
#line 4 "cp-algo/number_theory/modint.hpp"
#include <iostream>
#include <cassert>
namespace cp_algo::math {
inline constexpr auto inv2(auto x) {
assert(x % 2);
std::make_unsigned_t<decltype(x)> y = 1;
while(y * x != 1) {
y *= 2 - x * y;
}
return y;
}
template<typename modint, typename _Int>
struct modint_base {
using Int = _Int;
using UInt = std::make_unsigned_t<Int>;
static constexpr size_t bits = sizeof(Int) * 8;
using Int2 = std::conditional_t<bits <= 32, int64_t, __int128_t>;
using UInt2 = std::conditional_t<bits <= 32, uint64_t, __uint128_t>;
static Int mod() {
return modint::mod();
}
static UInt imod() {
return modint::imod();
}
static UInt2 pw128() {
return modint::pw128();
}
static UInt m_reduce(UInt2 ab) {
if(mod() % 2 == 0) [[unlikely]] {
return UInt(ab % mod());
} else {
UInt2 m = (UInt)ab * imod();
return UInt((ab + m * mod()) >> bits);
}
}
static UInt m_transform(UInt a) {
if(mod() % 2 == 0) [[unlikely]] {
return a;
} else {
return m_reduce(a * pw128());
}
}
modint_base(): r(0) {}
modint_base(Int2 rr): r(UInt(rr % mod())) {
r = std::min(r, r + mod());
r = m_transform(r);
}
modint inv() const {
return bpow(to_modint(), mod() - 2);
}
modint operator - () const {
modint neg;
neg.r = std::min(-r, 2 * mod() - r);
return neg;
}
modint& operator /= (const modint &t) {
return to_modint() *= t.inv();
}
modint& operator *= (const modint &t) {
r = m_reduce((UInt2)r * t.r);
return to_modint();
}
modint& operator += (const modint &t) {
r += t.r; r = std::min(r, r - 2 * mod());
return to_modint();
}
modint& operator -= (const modint &t) {
r -= t.r; r = std::min(r, r + 2 * mod());
return to_modint();
}
modint operator + (const modint &t) const {return modint(to_modint()) += t;}
modint operator - (const modint &t) const {return modint(to_modint()) -= t;}
modint operator * (const modint &t) const {return modint(to_modint()) *= t;}
modint operator / (const modint &t) const {return modint(to_modint()) /= t;}
// Why <=> doesn't work?..
auto operator == (const modint_base &t) const {return getr() == t.getr();}
auto operator != (const modint_base &t) const {return getr() != t.getr();}
auto operator <= (const modint_base &t) const {return getr() <= t.getr();}
auto operator >= (const modint_base &t) const {return getr() >= t.getr();}
auto operator < (const modint_base &t) const {return getr() < t.getr();}
auto operator > (const modint_base &t) const {return getr() > t.getr();}
Int rem() const {
UInt R = getr();
return 2 * R > (UInt)mod() ? R - mod() : R;
}
// Only use if you really know what you're doing!
UInt modmod() const {return (UInt)8 * mod() * mod();};
void add_unsafe(UInt t) {r += t;}
void pseudonormalize() {r = std::min(r, r - modmod());}
modint const& normalize() {
if(r >= (UInt)mod()) {
r %= mod();
}
return to_modint();
}
void setr(UInt rr) {r = m_transform(rr);}
UInt getr() const {
UInt res = m_reduce(r);
return std::min(res, res - mod());
}
void setr_direct(UInt rr) {r = rr;}
UInt getr_direct() const {return r;}
private:
UInt r;
modint& to_modint() {return static_cast<modint&>(*this);}
modint const& to_modint() const {return static_cast<modint const&>(*this);}
};
template<typename modint>
concept modint_type = std::is_base_of_v<modint_base<modint, typename modint::Int>, modint>;
template<modint_type modint>
std::istream& operator >> (std::istream &in, modint &x) {
typename modint::UInt r;
auto &res = in >> r;
x.setr(r);
return res;
}
template<modint_type modint>
std::ostream& operator << (std::ostream &out, modint const& x) {
return out << x.getr();
}
template<auto m>
struct modint: modint_base<modint<m>, decltype(m)> {
using Base = modint_base<modint<m>, decltype(m)>;
using Base::Base;
static constexpr Base::UInt im = m % 2 ? inv2(-m) : 0;
static constexpr Base::UInt r2 = (typename Base::UInt2)(-1) % m + 1;
static constexpr Base::Int mod() {return m;}
static constexpr Base::UInt imod() {return im;}
static constexpr Base::UInt2 pw128() {return r2;}
};
template<typename Int = int64_t>
struct dynamic_modint: modint_base<dynamic_modint<Int>, Int> {
using Base = modint_base<dynamic_modint<Int>, Int>;
using Base::Base;
static Int mod() {return m;}
static Base::UInt imod() {return im;}
static Base::UInt2 pw128() {return r2;}
static void switch_mod(Int nm) {
m = nm;
im = m % 2 ? inv2(-m) : 0;
r2 = static_cast<Base::UInt>(static_cast<Base::UInt2>(-1) % m + 1);
}
// Wrapper for temp switching
auto static with_mod(Int tmp, auto callback) {
struct scoped {
Int prev = mod();
~scoped() {switch_mod(prev);}
} _;
switch_mod(tmp);
return callback();
}
private:
static thread_local Int m;
static thread_local Base::UInt im, r2;
};
template<typename Int>
Int thread_local dynamic_modint<Int>::m = 1;
template<typename Int>
dynamic_modint<Int>::Base::UInt thread_local dynamic_modint<Int>::im = -1;
template<typename Int>
dynamic_modint<Int>::Base::UInt thread_local dynamic_modint<Int>::r2 = 0;
}
#line 1 "cp-algo/math/cvector.hpp"
#include <algorithm>
#line 5 "cp-algo/math/cvector.hpp"
#include <complex>
#include <vector>
#include <ranges>
namespace cp_algo::math::fft {
using ftype = double;
static constexpr size_t bytes = 32;
static constexpr size_t flen = bytes / sizeof(ftype);
using point = std::complex<ftype>;
using vftype [[gnu::vector_size(bytes)]] = ftype;
using vpoint = std::complex<vftype>;
#define WITH_IV(...) \
[&]<size_t ... i>(std::index_sequence<i...>) { \
return __VA_ARGS__; \
}(std::make_index_sequence<flen>());
template<typename ft>
constexpr ft to_ft(auto x) {
return ft{} + x;
}
template<typename pt>
constexpr pt to_pt(point r) {
using ft = std::conditional_t<std::is_same_v<point, pt>, ftype, vftype>;
return {to_ft<ft>(r.real()), to_ft<ft>(r.imag())};
}
struct cvector {
static constexpr size_t pre_roots = 1 << 17;
std::vector<vftype> x, y;
cvector(size_t n) {
n = std::max(flen, std::bit_ceil(n));
x.resize(n / flen);
y.resize(n / flen);
}
template<class pt = point>
void set(size_t k, pt t) {
if constexpr(std::is_same_v<pt, point>) {
x[k / flen][k % flen] = real(t);
y[k / flen][k % flen] = imag(t);
} else {
x[k / flen] = real(t);
y[k / flen] = imag(t);
}
}
template<class pt = point>
pt get(size_t k) const {
if constexpr(std::is_same_v<pt, point>) {
return {x[k / flen][k % flen], y[k / flen][k % flen]};
} else {
return {x[k / flen], y[k / flen]};
}
}
vpoint vget(size_t k) const {
return get<vpoint>(k);
}
size_t size() const {
return flen * std::size(x);
}
void dot(cvector const& t) {
size_t n = size();
for(size_t k = 0; k < n; k += flen) {
set(k, get<vpoint>(k) * t.get<vpoint>(k));
}
}
static const cvector roots;
template<class pt = point>
static pt root(size_t n, size_t k) {
if(n < pre_roots) {
return roots.get<pt>(n + k);
} else {
auto arg = std::numbers::pi / ftype(n);
if constexpr(std::is_same_v<pt, point>) {
return {cos(ftype(k) * arg), sin(ftype(k) * arg)};
} else {
return WITH_IV(pt{vftype{cos(ftype(k + i) * arg)...},
vftype{sin(ftype(k + i) * arg)...}});
}
}
}
template<class pt = point>
static void exec_on_roots(size_t n, size_t m, auto &&callback) {
size_t step = sizeof(pt) / sizeof(point);
pt cur;
pt arg = to_pt<pt>(root<point>(n, step));
for(size_t i = 0; i < m; i += step) {
if(i % 64 == 0 || n < pre_roots) {
cur = root<pt>(n, i);
} else {
cur *= arg;
}
callback(i, cur);
}
}
void ifft() {
size_t n = size();
for(size_t i = 1; i < n; i *= 2) {
for(size_t j = 0; j < n; j += 2 * i) {
auto butterfly = [&]<class pt>(size_t k, pt rt) {
k += j;
auto t = get<pt>(k + i) * conj(rt);
set(k + i, get<pt>(k) - t);
set(k, get<pt>(k) + t);
};
if(2 * i <= flen) {
exec_on_roots(i, i, butterfly);
} else {
exec_on_roots<vpoint>(i, i, butterfly);
}
}
}
for(size_t k = 0; k < n; k += flen) {
set(k, get<vpoint>(k) /= to_pt<vpoint>(ftype(n)));
}
}
void fft() {
size_t n = size();
for(size_t i = n / 2; i >= 1; i /= 2) {
for(size_t j = 0; j < n; j += 2 * i) {
auto butterfly = [&]<class pt>(size_t k, pt rt) {
k += j;
auto A = get<pt>(k) + get<pt>(k + i);
auto B = get<pt>(k) - get<pt>(k + i);
set(k, A);
set(k + i, B * rt);
};
if(2 * i <= flen) {
exec_on_roots(i, i, butterfly);
} else {
exec_on_roots<vpoint>(i, i, butterfly);
}
}
}
}
};
const cvector cvector::roots = []() {
cvector res(pre_roots);
for(size_t n = 1; n < res.size(); n *= 2) {
auto base = std::polar(1., std::numbers::pi / ftype(n));
point cur = 1;
for(size_t k = 0; k < n; k++) {
if((k & 15) == 0) {
cur = std::polar(1., std::numbers::pi * ftype(k) / ftype(n));
}
res.set(n + k, cur);
cur *= base;
}
}
return res;
}();
template<typename base>
struct dft {
cvector A;
dft(std::vector<base> const& a, size_t n): A(n) {
for(size_t i = 0; i < std::min(n, a.size()); i++) {
A.set(i, a[i]);
}
if(n) {
A.fft();
}
}
std::vector<base> operator *= (dft const& B) {
assert(A.size() == B.A.size());
size_t n = A.size();
if(!n) {
return std::vector<base>();
}
A.dot(B.A);
A.ifft();
std::vector<base> res(n);
for(size_t k = 0; k < n; k++) {
res[k] = A.get(k);
}
return res;
}
auto operator * (dft const& B) const {
return dft(*this) *= B;
}
point operator [](int i) const {return A.get(i);}
};
}
#line 5 "cp-algo/math/fft.hpp"
namespace cp_algo::math::fft {
template<modint_type base>
struct dft<base> {
int split;
cvector A, B;
dft(auto const& a, size_t n): A(n), B(n) {
split = int(std::sqrt(base::mod())) + 1;
cvector::exec_on_roots(2 * n, size(a), [&](size_t i, point rt) {
size_t ti = std::min(i, i - n);
auto rem = std::remainder(a[i].rem(), split);
auto quo = (ftype(a[i].rem()) - rem) / split;
A.set(ti, A.get(ti) + rem * rt);
B.set(ti, B.get(ti) + quo * rt);
});
if(n) {
A.fft();
B.fft();
}
}
void mul(auto &&C, auto const& D, auto &res, size_t k) {
assert(A.size() == C.size());
size_t n = A.size();
if(!n) {
res = {};
return;
}
for(size_t i = 0; i < n; i += flen) {
auto tmp = A.vget(i) * D.vget(i) + B.vget(i) * C.vget(i);
A.set(i, A.vget(i) * C.vget(i));
B.set(i, B.vget(i) * D.vget(i));
C.set(i, tmp);
}
A.ifft();
B.ifft();
C.ifft();
auto splitsplit = (base(split) * split).rem();
cvector::exec_on_roots(2 * n, std::min(n, k), [&](size_t i, point rt) {
rt = conj(rt);
auto Ai = A.get(i) * rt;
auto Bi = B.get(i) * rt;
auto Ci = C.get(i) * rt;
int64_t A0 = llround(real(Ai));
int64_t A1 = llround(real(Ci));
int64_t A2 = llround(real(Bi));
res[i] = A0 + A1 * split + A2 * splitsplit;
if(n + i >= k) {
return;
}
int64_t B0 = llround(imag(Ai));
int64_t B1 = llround(imag(Ci));
int64_t B2 = llround(imag(Bi));
res[n + i] = B0 + B1 * split + B2 * splitsplit;
});
}
void mul_inplace(auto &&B, auto& res, size_t k) {
mul(B.A, B.B, res, k);
}
void mul(auto const& B, auto& res, size_t k) {
mul(cvector(B.A), B.B, res, k);
}
std::vector<base> operator *= (dft &B) {
std::vector<base> res(2 * A.size());
mul_inplace(B, res, size(res));
return res;
}
std::vector<base> operator *= (dft const& B) {
std::vector<base> res(2 * A.size());
mul(B, res, size(res));
return res;
}
auto operator * (dft const& B) const {
return dft(*this) *= B;
}
point operator [](int i) const {return A.get(i);}
};
void mul_slow(auto &a, auto const& b, size_t k) {
if(empty(a) || empty(b)) {
a.clear();
} else {
size_t n = std::min(k, size(a));
size_t m = std::min(k, size(b));
a.resize(k);
for(int j = int(k - 1); j >= 0; j--) {
a[j] *= b[0];
for(int i = std::max(j - (int)n, 0) + 1; i < std::min(j + 1, (int)m); i++) {
a[j] += a[j - i] * b[i];
}
}
}
}
size_t com_size(size_t as, size_t bs) {
if(!as || !bs) {
return 0;
}
return std::max(flen, std::bit_ceil(as + bs - 1) / 2);
}
void mul_truncate(auto &a, auto const& b, size_t k) {
using base = std::decay_t<decltype(a[0])>;
if(std::min({k, size(a), size(b)}) < magic) {
mul_slow(a, b, k);
return;
}
auto n = std::max(flen, std::bit_ceil(
std::min(k, size(a)) + std::min(k, size(b)) - 1
) / 2);
a.resize(k);
auto A = dft<base>(a, n);
if(&a == &b) {
A.mul(A, a, k);
} else {
A.mul_inplace(dft<base>(b | std::views::take(k), n), a, k);
}
}
void mul(auto &a, auto const& b) {
if(size(a)) {
mul_truncate(a, b, size(a) + size(b) - 1);
}
}
}