先把代码发上来
#include<iostream>
#include<algorithm>
#include<functional>
using i32 = int;
using u32 = unsigned;
using i64 = long long;
using u64 = unsigned long long;
#define IL __inline__ __attribute__((always_inline))
#define RC(T, x) reinterpret_cast<T>(x)
/*
一个高性能的多项式实现
作者为 QedDust413 & killer_joke
*/
namespace Setting{
constexpr u32 mod = 998244353;
constexpr size_t sta_l_MB = 64;
constexpr int Detail = 1;
}
namespace __qed{
constexpr int bit_up(int x){return 1 << (32 - __builtin_clz(x));}
constexpr int crz_32(int x){return __builtin_ctz(x);}
constexpr int cro_32(int x){return __builtin_ctz(~x);}
constexpr int lg2u(int x){return x == 1 ? 0 : (32 - __builtin_clz(x - 1));}
constexpr int lg2d(int x){return std::__lg(x);}
constexpr bool ispow2(u64 x){return (x != 0) && ((x & (x - 1)) == 0);}
constexpr u32 primitive_root(u32 M){
u32 qed = 0, n = M - 1, d[11] = {};
for (int i = 2; i * i <= n; i++){
if (n % i == 0){
d[qed++] = i;
do{n /= i;}while (n % i == 0);
}
}
if (n > 1){d[qed++] = n;}
for (u32 g = 2, r = 0; ; ++g){
for (u32 i = 0; i < qed; ++i){
u32 b = (M - 1) / d[i], a = g;
for(r = 1; b; b >>= 1, a = u64(a) * a % M){
if(b & 1){r = u64(r) * a % M;}
}
if (r == 1){break;}
}
if(r != 1){return g;}
}
}
}
namespace Stat_Info{
using Setting::Detail;
i64 ntt_size;
i64 fill0_size, copy_size, rev_size;
size_t max_cost_sta_l;
constexpr const char* Author = "QedDust413 & killer_joke";
constexpr const char* Thanks = "negiizhao, chaihf, rogeryoungh, EntropyIncreaser, hly1204, Min_25, Qdedawsd2233, bh1234666, yosupo, Pulsating_Dust, KKKKa, qdc, and more.";
template<typename Tf = std::ostream>void report(Tf &outf = std::clog){
if constexpr(Detail <= 0){outf << "Statistics are turned off.\n";}
if constexpr(Detail > 0){outf << "ntt_size:" << ntt_size << "\n";}
if constexpr(Detail > 1){outf << "fill0_size:" << fill0_size << " copy_size:" << copy_size << " rev_size:" << rev_size << "\nmax_cost_sta:" << (double(max_cost_sta_l) / double(1 << 20)) << "\n";}
outf.flush();
}
}
namespace mem_helper{
using namespace __qed;
template<size_t align>void *to_align(void *mem){
static_assert(ispow2(align));
return RC(void*, (RC(u64, mem) + (align - 1)) & (-align));
}
template<size_t align>bool is_align(const void *mem){
static_assert(ispow2(align));
return (RC(u64, mem) & (align - 1)) == 0;
}
char _mem[Setting::sta_l_MB << 20], *_now = _mem;
struct pre_aloc{
char *t;
pre_aloc():t(_now){}
~pre_aloc(){_now = t;}
};
template<typename T, size_t align>T *AlocS(size_t n){
T *r = RC(T*, to_align<align>(_now));
_now = RC(char*, r + n);
if constexpr(Stat_Info::Detail > 1){
Stat_Info::max_cost_sta_l = std::max<size_t>(_now - _mem, Stat_Info::max_cost_sta_l);
}
return r;
}
template<size_t align>void *_alocH(size_t n){
constexpr size_t ptr_size = sizeof(void*), Extra = (align + ptr_size);
void *base = std::malloc(n + Extra), *p = to_align<align>(RC(char*, base) + ptr_size);
return RC(void**, p)[-1] = base, p;
}
template<typename T, size_t align>T *AlocH(size_t n){
return RC(T*, _alocH<align>(n * sizeof(T)));
}
void Freed(void* p){
std::free(RC(void**, p)[-1]);
}
template<typename T, size_t align>struct Trivial_Aligned_Allocator{
static_assert(std::is_trivial<T>::value);
typedef T value_type;
T *allocate(std::size_t n) {return AlocH<T, align>(n);}
template <class T1>struct rebind {using other = Trivial_Aligned_Allocator<T1, align>;};
void deallocate(T *p, std::size_t n) {Freed(p);}
};
}
namespace Montgo{
//Author : QedDust413
struct Mont32{
u32 Mod, Mod2, Inv, NInv, R2;
constexpr Mont32(u32 n):Mod(n), Mod2(n << 1), Inv(1), NInv(), R2((-u64(n)) % n){
for (int i = 0; i < 5; ++i){Inv *= 2 - n * Inv;}NInv = -Inv;
}
constexpr IL u32 reduce (u64 x)const{return (x + u64(u32(x) * NInv) * Mod) >> 32;}
constexpr IL u32 reduce_s (u64 x)const{
u32 r = (x >> 32) - ((u64(u32(x) * Inv) * Mod) >> 32);
return r >> 31 ? r + Mod : r;
}
constexpr IL u32 mul(u32 x, u32 y)const{return reduce(u64(x) * y);}
constexpr IL u32 mul_s(u32 x, u32 y)const{return reduce_s(u64(x) * y);}
constexpr IL u32 In(u32 x)const{return mul(x, R2);}
constexpr IL u32 In_s(u32 x)const{return mul_s(x, R2);}
constexpr IL u32 Out(u32 x)const{
u32 r = (x + (u64(x * NInv) * Mod)) >> 32;
return r < Mod ? r : r - Mod;
}
};
}
namespace field_Z{
//Author : QedDust413
using Setting::mod;
constexpr Montgo::Mont32 Space(mod);
constexpr u32 mod2 = Space.Mod2;
constexpr IL u32 shrink(u32 x){return x >= mod ? x - mod : x;}
constexpr IL u32 dilate2(u32 x){return x >> 31 ? x + mod2 : x;}
using Z = u32;
constexpr bool isgood(Z x){return x < mod2;}
constexpr IL Z InZ(u32 x){return Space.In(x);}
constexpr IL Z InZs(u32 x){return Space.In_s(x);}
constexpr Z zero_Z(0), one_Z(InZs(1)), not_exist_Z(-1);
constexpr IL u32 OutZ(Z x){return Space.Out(x);}
constexpr bool equalZ(Z x, Z y){return shrink(x) == shrink(y);}
namespace calc{
constexpr IL Z addZ(Z x, Z y){return dilate2(x + y - mod2);}
constexpr IL Z subZ(Z x, Z y){return dilate2(x - y);}
constexpr IL Z mulZ(Z x, Z y){return Space.mul(x, y);}
constexpr Z powZ(Z a, u32 b, Z r = one_Z){
for(; b; a = mulZ(a, a), b >>= 1){if(b & 1){r = mulZ(r, a);}}
return r;
}
constexpr IL Z invZ(Z x){return powZ(x, mod - 2);}
constexpr IL Z divZ(Z x, Z y){return powZ(y, mod - 2, x);}
template<bool strict = true>constexpr IL Z mulZs(Z x, Z y){
if constexpr(strict){return Space.mul_s(x, y);}
return mulZ(x, y);
}
constexpr IL Z negZ(Z x){return mod2 - x;}
namespace extend_{
constexpr Z absZ(Z x){
u32 r = OutZ(x);
return InZs(std::min(r, mod - r));
}
constexpr Z legendre(Z x){return powZ(x, (mod - 1) >> 1);}
constexpr bool isQR(Z x){return equalZ(legendre(x), one_Z);}
constexpr Z sqrtZ(Z x){
if(!isQR(x)){return not_exist_Z;}
Z a(1), I_mul(0);
while(isQR(I_mul = subZ(mulZ(a, a), x))){++a;}
struct dZ{
Z r, i;
constexpr void Mul(dZ d, Z I_){*this = {addZ(mulZ(r, d.r), mulZ(mulZ(I_, i), d.i)), addZ(mulZ(r, d.i), mulZ(i, d.r))};}
};
dZ s = {a, one_Z}, r = {one_Z, zero_Z};
for(u32 b = (mod + 1) >> 1; b; b >>= 1, s.Mul(s, I_mul)){if(b & 1){r.Mul(s, I_mul);}}
return absZ(r.r);
}
}
}
template<int fixes, bool strict = true>constexpr Z trans(Z x){
constexpr Z o = fixes > 0 ? calc::powZ(Space.R2, fixes) : calc::powZ(1, -fixes);
return calc::mulZs<strict>(x, o);
}
template<bool strict = true>constexpr Z trans(Z x, int fixes){
return calc::mulZs<strict>(x, fixes > 0 ? calc::powZ(Space.R2, fixes) : calc::powZ(1, -fixes));
}
namespace Const{
constexpr Z _half = InZs((mod + 1) >> 1), _neghalf = InZs((mod - 1) >> 1), neg_one = InZs(mod - 1);
constexpr Z img = calc::extend_::sqrtZ(neg_one), imgI = mod - img, _imghalf = calc::mulZs(_half, img);
}
}
#pragma GCC target("avx2")
#include <immintrin.h>
#define VEC(sz, T) __attribute((vector_size(sz))) T
namespace SIMD{
using i32x8 = VEC(32, i32);
using u32x8 = VEC(32, u32);
using i64x4 = VEC(32, i64);
using u64x4 = VEC(32, u64);
using I256 = __m256i;
constexpr IL u32x8 setu32x8(u32 x){return (u32x8){x, x, x, x, x, x, x, x};}
template<int typ>IL u32x8 shuffle(const u32x8 &x){return RC(u32x8, _mm256_shuffle_epi32(RC(I256, x), typ));}
template<int typ>IL u32x8 blend(const u32x8 &x, const u32x8 &y){return RC(u32x8, _mm256_blend_epi32(RC(I256, x), RC(I256, y), typ));}
IL I256 swaplohi128(const I256 &x){return _mm256_permute2x128_si256(x, x, 1);}
IL u32x8& x8(u32 *data){return *RC(u32x8* ,data);}
IL const u32x8& x8(const u32 *data){return *RC(const u32x8*, data);}
IL I256 loadu(const void* data){return _mm256_loadu_si256(RC(const __m256i_u*, data));}
IL void storeu(const I256 &x, void* data){return _mm256_storeu_si256(RC(__m256i_u*, data), x);}
IL u64x4 fus_mul(const u32x8 &x, const u32x8 &y){return RC(u64x4, _mm256_mul_epu32(RC(I256, x), RC(I256, y)));}
}
namespace field_Z{
using SIMD::x8;
using SIMD::u32x8;
using SIMD::setu32x8;
using Zx8 = u32x8;
constexpr u32x8 modx8 = setu32x8(mod), mod2x8 = setu32x8(mod2), NInvx8 = setu32x8(Space.NInv);
constexpr Zx8 one_Zx8 = setu32x8(one_Z), zerox8 = setu32x8(0u);
IL Zx8 dilate2x8(const Zx8 &x){return x + (RC(Zx8, RC(SIMD::i32x8, x) < RC(SIMD::i32x8, zerox8)) & mod2x8);}
IL Zx8 shrinkx8(const Zx8 &x){return x - ((x >= modx8) & modx8);}
namespace calc{
//Author : killer_joke
using namespace SIMD;
IL Zx8 addZx8(const Zx8 &x, const Zx8 &y){return dilate2x8(x + y - mod2x8);}
IL Zx8 subZx8(const Zx8 &x, const Zx8 &y){return dilate2x8(x - y);}
IL Zx8 mulZx8(const Zx8 &x, const Zx8 &y){
u32x8 z = (NInvx8 * x * y);
return blend<0xaa>(RC(u32x8, (fus_mul(x, y) + fus_mul(z, modx8)) >> 32), RC(u32x8, (fus_mul(u32x8(u64x4(x) >> 32), u32x8(u64x4(y) >> 32)) + fus_mul(shuffle<0xf5>(z), modx8))));
}
IL Zx8 powZx8(const Zx8 &y, u32 b, const Zx8 &_r = one_Zx8){
Zx8 x = y, r = _r;
for(; b; x = mulZx8(x, x), b >>= 1){if(b & 1){r = mulZx8(r, x);}}
return r;
}
IL Zx8 invZx8(const Zx8 &x){return powZx8(x, mod - 2);}
IL Zx8 divZx8(const Zx8 &x, const Zx8 &y){return powZx8(y, mod - 2, x);}
template<bool strict = true>IL Zx8 mulZsx8(const Zx8 &x, const Zx8 &y){
if constexpr (strict){
u32x8 z = (NInvx8 * x * y);
z = blend<0xaa>(RC(u32x8, (fus_mul(x, y) + fus_mul(z, modx8)) >> 32), RC(u32x8, (fus_mul(u32x8(u64x4(x) >> 32), u32x8(u64x4(y) >> 32)) + fus_mul(shuffle<0xf5>(z), modx8)))) - modx8;
return z + (RC(Zx8, RC(i32x8, z) < RC(i32x8, zerox8)) & modx8);
}
return mulZx8(x, y);
}
IL Zx8 negZx8(const Zx8 &x){return mod2x8 - x;}
IL Zx8 _AmulZx8(const Zx8 &a, const Zx8 &b, const Zx8 &c){return mulZx8(a + b, c);}
IL Zx8 _SmulZx8(const Zx8 &a, const Zx8 &b, const Zx8 &c){return mulZx8(a - b + mod2x8, c);}
template<int typ>IL u32x8 Neg(const u32x8 &x){return blend<typ>(x, mod2x8 - x);}
}
}
#undef VEC
#define STATI(ifo, dt, num) if constexpr (Stat_Info::Detail > dt){Stat_Info::ifo += (num);}
namespace poly_base{
using namespace field_Z;
using __qed::bit_up;
constexpr int mp2 = __qed::crz_32(mod - 1);
void fl0(Z *f, int l){STATI(fill0_size, 1, l)std::fill_n(f, l, zero_Z);}
void fl0(Z *bg, Z *ed){STATI(fill0_size, 1, ed - bg)std::fill(bg, ed, zero_Z);}
void Cpy(const Z *f, int l, Z *g){STATI(copy_size, 1, l)std::copy_n(f, l, g);}
void Cpy(const Z *bg, const Z *ed, Z *g){STATI(copy_size, 1, ed - bg)std::copy(bg, ed, g);}
void rev(Z *bg, Z *ed){STATI(rev_size, 1, ed - bg)std::reverse(bg, ed);}
void rev(Z *f, int l){STATI(rev_size, 1, l)std::reverse(f, f + l);}
void Cpy_fl0(const Z *f, int n, Z *g, int lim){Cpy(f, n, g), fl0(g + n, g + lim);}
void Cpy_rev(const Z *f, int l, Z *g){STATI(rev_size, 1, l)STATI(copy_size, 1, l)std::reverse_copy(f, f + l, g);}
#define flx(nam, opt) void nam(const Z *A, int l, Z *B, Z t){int i=0;if(l>16){Zx8 tx8=setu32x8(t);for(;i+7<l;i+=8){x8(B+i)=calc::opt##Zx8(x8(A+i),tx8);}}for(;i<l;++i){B[i]=calc::opt##Z(A[i],t);}}
flx(mul_t_s, mul)
flx(add_t_s, add)
flx(sub_t_s, sub)
#undef flx
#define flx(nam, opt) void nam(Z *A,int l,const Z *B){int i=0;for(;i+7<l;i+=8){x8(A + i)=calc::opt##Zx8(x8(A+i),x8(B+i));}for(;i<l;++i){A[i]=calc::opt##Z(A[i],B[i]);}} void nam(const Z *A,const Z *B,int l,Z *C){int i=0;for(;i+7<l;i+=8){x8(C+i)=calc::opt##Zx8(x8(A+i),x8(B+i));}for(;i<l;++i){C[i]=calc::opt##Z(A[i],B[i]);}}
flx(dot, mul)
flx(addP, add)
flx(subP, sub)
#undef flx
void NegP(const Z *A, int l, int r, Z *B){
int i = l;
if(r - l >= 16){
for(; i & 7; ++i){B[i] = calc::negZ(A[i]);}
for(; i + 7 < r; i += 8){x8(B + i) = calc::negZx8(x8(A + i));}
}
for(; i < r; ++i){B[i] = calc::negZ(A[i]);}
}
}
namespace poly_base{
using mem_helper::pre_aloc;
using mem_helper::Freed;
const std::function<Z*(size_t)> alocS = mem_helper::AlocS<Z, 32>;
const std::function<Z*(size_t)> alocH = mem_helper::AlocH<Z, 32>;
struct pre_base{
const std::function<void(Z*, int, int)> jok;
Z *p;
int sz;
pre_base(std::function<void(Z*, int, int)> f) : jok(f), p(nullptr), sz(0){
}
~pre_base(){
if(p != nullptr){
Freed(p);
}
}
void expand(int len){
if(len > sz){
len = (len + 7) & -8;
Z *nxtp = alocH(len);
if(p != nullptr){
Cpy(p, sz, nxtp), Freed(p);
}
jok(p = nxtp, sz, len), sz = len;
}
}
void clear(){
if(p != nullptr){
Freed(p), sz = 0;
}
}
const Z *operator()(int len){
return expand(len), p;
}
const Z *operator()(){
return p;
}
};
void multi_iv_Zx8(const Zx8 *g, int n, Zx8 *f){
using namespace calc;
if(n == 0){
return ;
}
f[0] = g[0];
for(int i = 1; i < n; ++i){
f[i] = mulZx8(f[i - 1], g[i]);
}
f[n - 1] = invZx8(f[n - 1]);
for(int i = n - 1; i; --i){
Zx8 ivi = f[i];
f[i] = mulZx8(ivi, f[i - 1]);
f[i - 1] = mulZx8(ivi, g[i]);
}
}
void multi_iv(const Z *g, int n, Z *f){
int N = n >> 3, i = N << 3;
multi_iv_Zx8(RC(const Zx8*, g), N, RC(Zx8*, f));
for(; i < n; ++i){
f[i] = calc::invZ(g[i]);
}
}
namespace cal_helper{
using namespace calc;
void _cal_i(Z *dus, int l, int r){
for(; l < 8; ++l){
dus[l] = InZ(l);
}
constexpr Zx8 eight_Zx8 = setu32x8(InZs(8));
for(int i = l; i < r; i += 8){
x8(dus + i) = addZx8(x8(dus + i - 8), eight_Zx8);
}
}
pre_base seqi(_cal_i);
void _cal_iv(Z *kil, int l, int r){
const Z* _i = seqi(r);
if(l < 8){
x8(kil) = invZx8(x8(_i)), l = 8;
}
multi_iv_Zx8(RC(const Zx8*, _i + l), (r - l) >> 3, RC(Zx8*, kil + l));
}
pre_base seqiv(_cal_iv);
void deriv(const Z *g, int n, Z *f){
const Z *_i = seqi(n + 1);
int i = 0;
if(n > 7){
for(; i < 7; ++i){f[i] = mulZ(g[i + 1], _i[i + 1]);}
for(; i + 7 < n; i += 8){storeu(RC(I256, mulZx8(x8(g + i + 1), x8(_i + i + 1))), f + i);}
}
for(; i < n; ++i){f[i] = mulZ(g[i + 1], _i[i + 1]);}
f[n] = zero_Z;
}
void integ(const Z *g, int n, Z *f){
const Z *_iv = seqiv(n + 1);
int i = n - 1;
if(i > 7){
for(; (i & 7) != 6; --i){f[i + 1] = mulZ(g[i], _iv[i + 1]);}
for(i -= 7; ~i; i -= 8){x8(f + i + 1) = mulZx8(RC(Zx8, loadu(g + i)), x8(_iv + i + 1));}
i += 7;
}
for(; ~i; --i){f[i + 1] = mulZ(g[i], _iv[i + 1]);}
f[0] = zero_Z;
}
}
using cal_helper::seqi;
using cal_helper::seqiv;
using cal_helper::deriv;
using cal_helper::integ;
}
namespace poly_base{
namespace Tools{
int compP(const Z *f, int l, const Z *g){
for(int i = 0; i < l; ++i){if(!equalZ(f[i], g[i])){return i;}}
return -1;
}
template<int fixes = 1, typename Tf = std::istream>void scanP(Z *f, int l, Tf &inf = std::cin){
for(int i = 0; i < l; ++i){inf >> f[i], f[i] = trans<fixes, false>(f[i]);}
}
template<int fixes = -1, typename Tf = std::ostream>void printP(const Z *f, int l, Tf &outf = std::cout){
if(l){
outf << trans<fixes>(f[0]);
for(int i = 1; i < l; ++i){outf << ' ' << trans<fixes>(f[i]);}
outf << '\n';
}
}
template<typename T = u32, typename lT = u64>T stoi_m(const std::string &s, const T m){
T r(0);
for(auto ch:s){r = ((lT)r * 10 + (ch - '0')) % m;}
return r;
}
int clzP(const Z *f, int n){
int i = 0;
while(i < n && equalZ(f[i], zero_Z)){++i;}return i;
}
int crzP(const Z *f, int n){
int i = n;
while(i && equalZ(f[i - 1], zero_Z)){--i;}return n - i;
}
}
}
namespace poly_base{
//Author : QedDust413
namespace f_n_t_t_base{
using namespace calc;
using namespace __qed;
constexpr int base4thre = 4;
template<bool strict = false>void mul_t(Z *A, int l, Z t){
for(int j = 0; j < l; ++j){A[j] = mulZs<strict>(A[j], t);}
}
constexpr Z _g(InZ(primitive_root(mod)));
struct P_R_Tab{
Z t[mp2 + 1];
constexpr P_R_Tab(Z G) : t(){
t[mp2] = shrink(powZ(G, (mod - 1) >> mp2));
for(int i = mp2 - 1; ~i; --i){t[i] = mulZs(t[i + 1], t[i + 1]);}
}
constexpr Z operator [] (int i) const {
return t[i];
}
}constexpr rt1(_g), rt1_I(invZ(_g));
template<int lim>struct omega_info_base2{
Z o[lim >> 1];
constexpr omega_info_base2(const P_R_Tab &Tb) : o(){
o[0] = one_Z;
for(int i = 1, l = lim >> 1; i < l; ++i){
o[i] = mulZ(o[i & (i - 1)], Tb[crz_32(i) + 2]);
}
}
constexpr Z operator [] (int i) const {
return o[i];
}
};
const omega_info_base2<base4thre> w(rt1), wI(rt1_I);
constexpr Zx8 powXx8(Z X){
Z X2 = mulZs(X, X), X3 = mulZs(X2, X), X4 = mulZs(X3, X), X5 = mulZs(X4, X), X6 = mulZs(X5, X), X7 = mulZs(X6, X);
return (Zx8){one_Z, X, X2, X3, X4, X5, X6, X7};
}
struct ntt_info_base4x8{
Z rt3[mp2 - 2], rt3_I[mp2 - 2];
Zx8 rt4ix8[mp2 - 3], rt4ix8_I[mp2 - 3];
constexpr ntt_info_base4x8() : rt3(), rt3_I(), rt4ix8(), rt4ix8_I(){
Z pr = one_Z, pr_I = one_Z;
for(int i = 0; i < mp2 - 2; ++i){
rt3[i] = mulZs(pr, rt1[i + 3]), rt3_I[i] = mulZs(pr_I, rt1_I[i + 3]);
pr = mulZs(pr, rt1_I[i + 3]), pr_I = mulZs(pr_I, rt1[i + 3]);
}
pr = one_Z, pr_I = one_Z;
for(int i = 0; i < mp2 - 3; ++i){
rt4ix8[i] = powXx8(mulZs(pr, rt1[i + 4])), rt4ix8_I[i] = powXx8(mulZs(pr_I, rt1_I[i + 4]));
pr = mulZs(pr, rt1_I[i + 4]), pr_I = mulZs(pr_I, rt1[i + 4]);
}
}
};
}
//fast_number_theoretic_transform
//Author : QedDust413
namespace f_n_t_t{
using namespace f_n_t_t_base;
template<bool strict = false, int fixes = 0>void dif_base2(Z *A, int lim){
STATI(ntt_size, 0, lim)
for(int L = lim >> 1, R = lim; L; L >>= 1, R >>= 1){
for(int i = 0, k = 0; i < lim; i += R, ++k){
for(int j = 0; j < L; ++j){
Z x = dilate2(A[i + j] - mod2) , y = mulZ(w[k], A[i + j + L]);
A[i + j] = x + y, A[i + j + L] = x - y + mod2;
}
}
}
if constexpr(fixes){mul_t<strict>(A, lim, trans<fixes>(one_Z));}
else{
for(int j = 0; j < lim; ++j){
A[j] = dilate2(A[j] - mod2);
if constexpr(strict){A[j] = shrink(A[j]);}
}
}
}
template<bool strict = false, int fixes = 0>void dit_base2(Z *A, int lim){
STATI(ntt_size, 0, lim)
for(int L = 1, R = 2; L < lim; L <<= 1, R <<= 1){
for(int i = 0, k = 0; i < lim; i += R, ++k){
for(int j = 0; j < L; ++j){
Z x = A[i + j], y = A[i + j + L];
A[i + j] = addZ(x, y), A[i + j + L] = mulZ(x - y + mod2, wI[k]);
}
}
}
mul_t<strict>(A, lim, trans<fixes + 1>(mod - ((mod - 1) >> crz_32(lim))));
}
constexpr Zx8 w14x8 = setu32x8(rt1[2]), w14Ix8 = setu32x8(rt1_I[2]);
constexpr Z w38 = mod - rt1_I[3], w38I = mod - rt1[3];
constexpr ntt_info_base4x8 iab4;
template<bool strict = false, int fixes = 0>void dif_base4x8(Z *A, int lim){
STATI(ntt_size, 0, lim)
int n = lim >> 3, L = n >> 1;
Zx8 *f = RC(Zx8*, A);
if(crz_32(n) & 1){
for(int j = 0; j < L; ++j){
Zx8 x = f[j], y = f[j + L];
f[j] = x + y, f[j + L] = x - y + mod2x8;
}
L >>= 1;
}
L >>= 1;
for(int R = L << 2; L; L >>= 2, R >>= 2){
Z _r = one_Z, _r2 = one_Z, _r3 = one_Z;
for(int i = 0, k = 0; i < n; i += R, ++k){
Zx8 r = setu32x8(_r), r2 = setu32x8(_r2), r3 = setu32x8(_r3);
for(int j = 0; j < L; ++j){
#define F(c) f[i + j + c * L]
Zx8 f0 = dilate2x8(F(0) - mod2x8), f1 = mulZx8(F(1), r), f2 = mulZx8(F(2), r2), f3 = mulZx8(F(3), r3);
Zx8 f1f3 = _SmulZx8(f1, f3, w14x8), f02 = addZx8(f0, f2), f13 = addZx8(f1, f3), f_02 = subZx8(f0, f2);
F(0) = f02 + f13, F(1) = f02 - f13 + mod2x8, F(2) = f_02 + f1f3, F(3) = f_02 - f1f3 + mod2x8;
#undef F
}
_r = mulZs(_r, iab4.rt3[cro_32(k)]), _r2 = mulZs(_r, _r), _r3 = mulZs(_r2, _r);
}
}
{
constexpr Zx8 _r = setu32x8(trans<fixes>(one_Z)), pr2 = {one_Z, one_Z, one_Z, rt1[2], one_Z, one_Z, one_Z, rt1[2]}, pr4 = {one_Z, one_Z, one_Z, one_Z, one_Z, rt1[3], rt1[2], w38};
Zx8 r = _r;
for(int i = 0; i < n; ++i){
Zx8& fi = f[i];
fi = mulZx8(fi, r);
fi = _AmulZx8(Neg<0xf0>(fi), RC(Zx8, swaplohi128(RC(I256, fi))), pr4);
fi = _AmulZx8(Neg<0xcc>(fi), shuffle<0x4e>(fi), pr2);
fi = addZx8(Neg<0xaa>(fi), shuffle<0xb1>(fi));
if constexpr(strict){fi = shrinkx8(fi);}
r = mulZsx8(r, iab4.rt4ix8[cro_32(i)]);
}
}
}
template<bool strict = false, int fixes = 0>void dit_base4x8(Z *A, int lim){
STATI(ntt_size, 0, lim)
int n = lim >> 3, L = 1;
Zx8 *f = RC(Zx8*, A);
{
constexpr Zx8 pr2 = {one_Z, one_Z, one_Z, rt1_I[2], one_Z, one_Z, one_Z, rt1_I[2]}, pr4 = {one_Z, one_Z, one_Z, one_Z, one_Z, rt1_I[3], rt1_I[2], w38I};
Zx8 r = setu32x8(trans<fixes + 1>(mod - ((mod - 1) >> crz_32(lim))));
for(int i = 0; i < n; ++i){
Zx8& fi = f[i];
fi = _AmulZx8(Neg<0xaa>(fi), shuffle<0xb1>(fi), pr2);
fi = _AmulZx8(Neg<0xcc>(fi), shuffle<0x4e>(fi), pr4);
fi = _AmulZx8(Neg<0xf0>(fi), RC(Zx8, swaplohi128(RC(I256, fi))), r);
r = mulZsx8(r, iab4.rt4ix8_I[cro_32(i)]);
}
}
for (int R = L << 2; L < (n >> 1) ; L <<= 2, R <<= 2){
Z _r = one_Z, _r2 = one_Z, _r3 = one_Z;
for(int i = 0, k = 0; i < n; i += R, ++k){
Zx8 r = setu32x8(_r), r2 = setu32x8(_r2), r3 = setu32x8(_r3);
for(int j = 0; j < L; ++j){
#define F(c) f[i + j + c * L]
Zx8 f0 = F(0), f1 = F(1), f2 = F(2), f3 = F(3);
Zx8 f2f3 = _SmulZx8(f2, f3, w14Ix8), f01 = addZx8(f0, f1), f23 = addZx8(f2, f3), f_01 = subZx8(f0, f1);
F(0) = addZx8(f01, f23), F(1) = _AmulZx8(f_01, f2f3 ,r), F(2) = _SmulZx8(f01, f23, r2), F(3) = _SmulZx8(f_01, f2f3, r3);
#undef F
}
_r = mulZs(_r, iab4.rt3_I[cro_32(k)]), _r2 = mulZs(_r, _r), _r3 = mulZs(_r2, _r);
}
}
if(L != n){
for(int j = 0; j < L; ++j){
Zx8 x = f[j], y = f[j + L];
f[j] = addZx8(x, y), f[j + L] = subZx8(x, y);
}
}
if constexpr (strict){for(int i = 0; i < n; ++i){f[i] = shrinkx8(f[i]);}}
}
template<bool strict = false, int fixes = 0>void dif(Z *A, int lim){
lim > base4thre ? dif_base4x8<strict, fixes>(A, lim) : dif_base2<strict, fixes>(A, lim);
}
template<bool strict = false, int fixes = 0>void dit(Z *A,int lim){
lim > base4thre ? dit_base4x8<strict, fixes>(A, lim) : dit_base2<strict, fixes>(A, lim);
}
}
using f_n_t_t::dif;
using f_n_t_t::dit;
void Conv(Z *f, int lim, Z *g){
dif(f, lim), dif(g, lim), dot(f, lim, g), dit(f, lim);
}
}
#undef STATI
namespace poly{
using namespace poly_base;
//Author : QedDust413
void Mul(const Z *f, int n, const Z *g, int m, Z *h){
int lim = bit_up(n + m);
pre_aloc mem;
Z *F = alocS(lim), *G = alocS(lim);
Cpy_fl0(f, n + 1, F, lim), Cpy_fl0(g, m + 1, G, lim), Conv(F, lim, G), Cpy(F, n + m + 1, h);
}
void Inv(const Z *g, int n, Z *f){
int lim = bit_up(n);
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
f[0] = calc::invZ(g[0]);
for(int t = 2, mid = 1; t <= lim; t <<= 1, mid <<= 1){
int xl = std::min<int>(t, n + 1);
Cpy_fl0(g, xl, o, t), Cpy_fl0(f, mid, h, t), Conv(o, t, h), fl0(o, mid), dif(o, t), dot(o, t, h), dit(o, t), NegP(o, mid, xl, f);
}
}
void _Div(const Z *g, int n, const Z *f, Z *q){
if(n == 0){
return q[0] = calc::divZ(g[0], f[0]), void();
}
int lim = bit_up(n) << 1;
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
Inv(f, n, o), fl0(o + n + 1, o + lim), Cpy_fl0(g, n + 1, h, lim), Conv(h, lim, o), Cpy(h, n + 1, q);
}
void Div(const Z *g, int n, const Z *f, Z *q){
using namespace calc;
if(n <= 64){return _Div(g, n, f, q);}
int bn = bit_up(n >> 4), bt = (n / bn) + 1, bn2 = bn << 1;
pre_aloc mem;
Z *o = alocS(bn2), *jok = alocS(bn2);
Inv(f, bn - 1, o), fl0(o + bn, bn), Cpy_fl0(g, bn, jok, bn2), Conv(jok, bn2, o);
Z *nf = alocS(bn2 * bt), *ng = alocS(bn2 * (bt - 1));
Cpy(jok, bn, q), Cpy_fl0(f, bn, nf, bn2), dif(nf, bn2);
for(int bi = 1; bi < bt; ++bi){
int xl = std::min<int>(bn, n + 1 - bn * bi);
Cpy_fl0(f + bi * bn, xl, nf + bi * bn2, bn2), dif(nf + bi * bn2, bn2);
Cpy_fl0(q + (bi - 1) * bn, bn, ng + (bi - 1) * bn2, bn2), dif<true>(ng + (bi - 1) * bn2, bn2), fl0(jok, bn2);
for(int j = 0; j < bi; ++j){
Z *nF = nf + (bi - j) * bn2, *nF1 = nF - bn2, *ngj = ng + j * bn2;
for(int i = 0; i < bn; i += 8){
x8(jok + i) = subZx8(x8(jok + i), _AmulZx8(x8(nF + i), x8(nF1 + i), x8(ngj + i)));
}
for(int i = bn; i < bn2; i += 8){
x8(jok + i) = subZx8(x8(jok + i), _SmulZx8(x8(nF + i), x8(nF1 + i), x8(ngj + i)));
}
}
dit(jok, bn2), fl0(jok + bn, bn), addP(jok, xl, g + bn * bi), dif(jok, bn2), dot(jok, bn2, o), dit(jok, bn2), Cpy(jok, xl, q + bn * bi);
}
}
void Div(const Z *g, int n, const Z *f, int m, Z *q){
int lim = n - m + 1, R = std::min<int>(m + 1, lim);
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
Cpy_rev(f + std::max<int>(2 * m - n, 0), R, o), fl0(o + R, o + lim);
Cpy_rev(g + m, lim, h), Div(h, n - m, o, q), rev(q, lim);
}
void Div(const Z *g, int n, const Z *f, int m, Z *q, Z *r){
Div(g, n, f, m, q);
int lim = bit_up(std::min<int>(n, m + m - 1));
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
Cpy_fl0(f, m, o, lim), Cpy_fl0(q, std::min<int>(n - m + 1, m), h, lim), Conv(o, lim, h), subP(g, o, m, r);
}
void Ln(const Z *g, int n, Z *f){
dot(g, seqi(n + 1), n + 1, f), Div(f, n, g, f), dot(f, n + 1, seqiv(n + 1));
}
template<bool calc_inv = true>void Expi(const Z *g, int n, Z *f, Z *h){
f[0] = h[0] = one_Z;
if(n < 1){return ;}
int lim = bit_up(n);
seqiv.expand(lim);
pre_aloc mem;
Z *o = alocS(lim), *A = alocS(lim), *B = alocS(lim);
A[0] = A[1] = one_Z;
for(int t = 2, mid = 1; t <= lim; t <<= 1, mid <<= 1){
int xl = std::min(n + 1, t), dl = std::min(t << 1, lim);
deriv(g, mid - 1, o), dif(o, mid), dot(o, mid, A), dit(o, mid), deriv(f, mid - 1, o + mid);
subP(o + mid, mid - 1, o), o[t - 1] = zero_Z, fl0(o, mid - 1), o[mid - 1] = calc::negZ(o[mid - 1]);
Cpy_fl0(h, mid, B, t), dif(B, t), dif(o, t), dot(o, t, B), dit(o, t);
integ(o, t - 1, o), fl0(o, mid), subP(o + mid, xl - mid, g + mid);
dif(o, t), dot(o, t, A), dit(o, t), NegP(o, mid, xl, f);
if(t != lim || calc_inv){
Cpy_fl0(f, xl, A, dl), dif(A, dl), dot(A, B, t, o), dit(o, t), fl0(o, mid);
dif(o, t), dot(o, t, B), dit(o, t), NegP(o, mid, xl, h);
}
}
}
void Exp(const Z *g, int n, Z *f){
using namespace calc;
pre_aloc mem;
if(n <= 64){
Z *h = alocS(n + 1);
return Expi<false>(g, n, f, h);
}
int bn = bit_up(n >> 4), bt = (n / bn) + 1, bn2 = bn << 1;
seqiv.expand(n + 1);
Z *h = alocS(bn2);
Expi(g, bn - 1, f, h);
fl0(h + bn, bn), dif(h, bn2);
Z *jok = alocS(bn2), *nf = alocS(bn2 * bt), *ng = alocS(bn2 * (bt - 1));
dot(g, seqi(), bn, nf), fl0(nf + bn, nf + bn2), dif(nf, bn2);
for(int bi = 1; bi < bt; ++bi){
int xl = std::min<int>(bn, n + 1 - bn * bi);
dot(g + bi * bn, seqi() + bi * bn, xl, nf + bi * bn2), fl0(nf + bi * bn2 + xl, bn2 - xl), dif(nf + bi * bn2, bn2);
Cpy_fl0(f + (bi - 1) * bn, bn, ng + (bi - 1) * bn2, bn2), dif<true>(ng + (bi - 1) * bn2, bn2), fl0(jok, bn2);
for(int j = 0; j < bi; ++j){
Z *nF = nf + (bi - j) * bn2, *nF1 = nF - bn2, *ngj = ng + j * bn2;
for(int i = 0; i < bn; i += 8){
x8(jok + i) = addZx8(x8(jok + i), _AmulZx8(x8(nF + i), x8(nF1 + i), x8(ngj + i)));
}
for(int i = bn; i < bn2; i += 8){
x8(jok + i) = addZx8(x8(jok + i), _SmulZx8(x8(nF + i), x8(nF1 + i), x8(ngj + i)));
}
}
dit(jok, bn2), fl0(jok + bn, bn), dif(jok, bn2), dot(jok, bn2, h), dit(jok, bn2), fl0(jok + bn, bn);
dot(jok, xl, seqiv() + bn * bi), dif(jok, bn2), dot(jok, bn2, ng), dit(jok, bn2), Cpy(jok, xl, f + bn * bi);
}
}
void InvSqrt(const Z *g, int n, Z *f){
int lim = bit_up(n);
pre_aloc mem;
Z *o = alocS(lim << 1), *h = alocS(lim << 1);
f[0] = calc::invZ(calc::extend_::sqrtZ(g[0]));
for(int t = 2, mid = 1; t <= lim; t <<= 1, mid <<= 1){
int xl = std::min<int>(t, n + 1);
Cpy_fl0(f, mid, o, t << 1), Cpy_fl0(g, xl, h, t << 1), dif(o, t << 1), dif(h, t << 1), dot(h, t << 1, o), dot(o, t << 1, o), dot(h, t << 1, o), dit(h, t << 1), mul_t_s(h + mid, xl - mid, f + mid, Const::_neghalf);
}
}
template<bool calc_inv = true>void Sqrti(const Z *g, int n, Z *f, Z *h){
int lim = bit_up(n);
f[0] = calc::extend_::sqrtZ(g[0]), h[0] = calc::invZ(f[0]);
if(n == 0)return ;
pre_aloc mem;
Z *o = alocS(lim), *H = alocS(lim), *ff = alocS(lim);
ff[0] = f[0];
for(int t = 2, mid = 1; t <= lim; t <<= 1, mid <<= 1){
int xl = std::min<int>(t, n + 1);
dot(ff, mid, ff), dit(ff, mid), subP(ff, g, mid, ff + mid), subP(ff + mid, g + mid, mid, ff + mid), fl0(ff, mid), Cpy_fl0(h, mid, H, t), Conv(ff, t, H), mul_t_s(ff + mid, xl - mid, f + mid, Const::_neghalf);
if(t != lim || calc_inv){
Cpy(f, t, o), dif(o, t), Cpy(o, t, ff), dot(o, t, H), dit(o, t), fl0(o, mid), dif(o, t), dot(o, t, H), dit(o, t), NegP(o, mid, xl, h);
}
}
}
void Sqrt(const Z *g, int n, Z *f){
using namespace calc;
pre_aloc mem;
if(n <= 64){
Z *h = alocS(n + 1);
return Sqrti<false>(g, n, f, h);
}
int bn = bit_up(n >> 4), bt = (n / bn) + 1, bn2 = bn << 1;
Z *o = alocS(bn2);
Sqrti(g, bn - 1, f, o), mul_t_s(o, bn, o, Const::_neghalf), fl0(o + bn, bn), dif(o, bn2);
Z *jok = alocS(bn2), *ng = alocS(bn2 * (bt - 1));
for(int bi = 1; bi < bt; ++bi){
int xl = std::min<int>(bn, n + 1 - bn * bi);
Cpy_fl0(f + (bi - 1) * bn, bn, ng + (bi - 1) * bn2, bn2), dif<true>(ng + (bi - 1) * bn2, bn2), fl0(jok, bn2);
{
Z* nG1 = ng + (bi - 1) * bn2;
for(int i = 0; i < bn; i += 8){
x8(jok + i) = addZx8(x8(jok + i), mulZx8(x8(nG1 + i), x8(ng + i)));
}
for(int i = bn; i < bn2; i += 8){
x8(jok + i) = subZx8(x8(jok + i), mulZx8(x8(nG1 + i), x8(ng + i)));
}
}
for(int j = 1; j < bi; ++j){
Z *nG = ng + (bi - j) * bn2, *nG1 = nG - bn2, *ngj = ng + j * bn2;
for(int i = 0; i < bn; i += 8){
x8(jok + i) = addZx8(x8(jok + i), _AmulZx8(x8(nG1 + i), x8(nG + i), x8(ngj + i)));
}
for(int i = bn; i < bn2; i += 8){
x8(jok + i) = subZx8(x8(jok + i), _SmulZx8(x8(nG1 + i), x8(nG + i), x8(ngj + i)));
}
}
dit(jok, bn2), fl0(jok + bn, bn), subP(jok, xl, g + bi * bn), dif(jok, bn2), dot(jok, bn2, o), dit(jok, bn2), Cpy(jok, xl, f + bi * bn);
}
}
void Pow(const Z *g, int n, int k, Z *f){
if(k == 0){f[0] = one_Z, fl0(f + 1, n);}
else if(k == 1){Cpy(g, n + 1, f);}
else{
pre_aloc mem;
Z *h = alocS(n + 1);
Ln(g, n, f), mul_t_s(f, n + 1, h, InZ(k)), Exp(h, n, f);
}
}
void Pow(const Z *g, int n, int k, int k2, Z *f){
if(g[0] != one_Z){
pre_aloc mem;
Z *h = alocS(n + 1);
mul_t_s(g, n + 1, f, calc::invZ(g[0])), Pow(f, n, k, h), mul_t_s(h, n + 1, f, calc::powZ(g[0], k2));
}
else{Pow(g, n, k, f);}
}
}
namespace poly{
constexpr int good_poly = 1, empty_poly = 0, not_exist_poly = -1;
int Safe_Sqrt(const Z *g, int n, Z *f){
int shift = Tools::clzP(g, n + 1);
if(shift > n){
return empty_poly;
}
if((shift & 1) || (!calc::extend_::isQR(g[shift]))){
return not_exist_poly;
}
if(shift){
pre_aloc mem;
Z *h = alocS(n + 1);
Cpy_fl0(g + shift, n - shift + 1, f, n - (shift >> 1) + 1), Sqrt(f, n - (shift >> 1), h);
Cpy(h, n - (shift >> 1) + 1, f + (shift >> 1)), fl0(f, shift >> 1);
}
else{
Sqrt(g, n, f);
}
return good_poly;
}
int Safe_Pow(const Z *g, int n, std::string s, Z *f){
using namespace Tools;
int shift = clzP(g, n + 1);
if(shift){
if((s.size() > size_t(8)) || ((stoll(s) * shift) > i64(n))){
return empty_poly;
}
int rf = (stoi(s) * shift), l = n + 1 - rf;
pre_aloc mem;
Z *h = alocS(l);
Cpy(g + shift, l, f), Pow(f, l - 1, stoi_m(s, mod), stoi_m(s, mod - 1), h), Cpy(h, l, f + rf), fl0(f, rf);
}
else{
Pow(g, n, stoi_m(s, mod), stoi_m(s, mod - 1), f);
}
return good_poly;
}
}
namespace poly{
//Author : killer_joke
namespace Trigonometric_Function{
void Trig(const Z *g, int n, Z *sint, Z *cost){
pre_aloc mem;
Z *A = alocS(n + 1), *o = alocS(n + 1), *h = alocS(n + 1);
mul_t_s(g, n + 1, A, Const::img), Expi<true>(A, n, h, o);
if(sint != nullptr){subP(o, h, n + 1, sint), mul_t_s(sint, n + 1, sint, Const::_imghalf);}
if(cost != nullptr){addP(h, o, n + 1, cost), mul_t_s(cost, n + 1, cost, Const::_half);}
}
void Sin(const Z *g, int n, Z *f){Trig(g, n, f, nullptr);}
void Cos(const Z *g, int n, Z *f){Trig(g, n, nullptr, f);}
void Tan(const Z *g, int n, Z *f){
pre_aloc mem;
Z* h = alocS(n + 1);
Trig(g, n, f, h), Div(f, n, h, f);
}
void aSin(const Z *g, int n, Z *f){
int lim = bit_up(n) << 1;
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
Cpy_fl0(g, n + 1, o, lim), dif(o, lim), dot(o, lim, o), dit(o, lim), NegP(o, 0, n + 1, o), o[0] = calc::addZ(one_Z, o[0]);
InvSqrt(o, n, h), deriv(g, n, o), fl0(o + n, o + lim), fl0(h + n + 1, h + lim), Conv(o, lim, h), integ(o, n, f);
}
void aCos(const Z* g, int n, Z* f){
aSin(g, n, f), NegP(f, 0, n + 1, f);
}
void aTan(const Z *g, int n, Z *f){
int lim = bit_up(n);
pre_aloc mem;
Z *o = alocS(lim << 1);
Cpy_fl0(g, n + 1, o, lim << 1), dif(o, lim << 1), dot(o, lim << 1, o), dit(o, lim << 1);
o[0] = calc::addZ(o[0], one_Z), deriv(g, n, o + lim), Div(o + lim, n, o, f), integ(f, n, f);
}
}
}
namespace poly_base{
namespace fac_helper{
void _fac_i(Z *_fac, int l, int r){
const Z *_i = seqi(r);
if(l == 0){_fac[0] = one_Z, l = 1;}
for(int i = l; i < r; ++i){_fac[i] = calc::mulZ(_fac[i - 1], _i[i]);}
}
void _fac_iv(Z *_ifac, int l, int r){
const Z *_iv = seqiv(r);
if(l == 0){_ifac[0] = one_Z, l = 1;}
for(int i = l; i < r; ++i){_ifac[i] = calc::mulZ(_ifac[i - 1], _iv[i]);}
}
pre_base seqfac(_fac_i), seqifac(_fac_iv);
}
using fac_helper::seqfac;
using fac_helper::seqifac;
}
namespace poly_base{
//Author : killer_joke
namespace makeseq{
using namespace calc;
void Ci(Z *f, int n, Z c){
Z fx = one_Z;
int i = 0;
for(; i < std::min(n, 8); ++i, fx = mulZ(fx, c)){f[i] = fx;}
if(n > 16){
Zx8 C8 = setu32x8(fx);
for(; i + 7 < n; i += 8){x8(f + i) = mulZx8(x8(f + i - 8), C8);}
}
for(; i < n; ++i){f[i] = mulZ(f[i - 1], c);}
}
void E_x(Z *f, int n){
const Z *ifac = seqifac(n);
int i = 0;
for(; i + 7 < n; i += 8){x8(f + i) = Neg<0xaa>(x8(ifac + i));}
for(; i < n; ++i){f[i] = (i & 1) ? negZ(ifac[i]) : ifac[i];}
}
void mRisingfac(Z *f, int n, Z m){
f[0] = one_Z;
for(int i = 1; i < n; ++i){f[i] = mulZ(f[i - 1], m), m = addZ(m, one_Z);}
}
void mFallingfac(Z *f, int n, Z m){
f[0] = one_Z;
for(int i = 1; i < n; ++i){f[i] = mulZ(f[i - 1], m), m = subZ(m, one_Z);}
}
}
}
namespace poly_base{
//Author : killer_joke
namespace f_n_t_t{
void _dfx2_f(Z *ka, int l, int r){
for(int i = std::max(1, l); i < r; i <<= 1){
makeseq::Ci(ka + i, i, rt1[__qed::lg2d(i) + 1]);
}
}
pre_base predifx2(_dfx2_f);
template<bool strict = false>void difx2(Z *f, int lim){
Cpy(f, lim, f + lim), dit(f + lim, lim);
dot(f + lim, lim, predifx2() + lim), dif<strict>(f + lim, lim);
}
template<bool strict = false>void difx2fxC(Z *f, int lim){
constexpr Z two_Z = InZs(2);
Cpy(f, lim, f + lim), dit(f + lim, lim);
dot(f + lim, lim, predifx2() + lim), f[lim] = subZ(f[lim], two_Z), dif<strict>(f + lim, lim);
}
void dif_neg(const Z *g, int lim, Z *f){
if(lim < 8){
if(lim == 1){
*f = *g;
}
else{
for(int i = 0; i < lim; ++i){
f[i] = g[i ^ 1];
}
}
}
else{
for(int i = 0; i < lim; i += 8){
x8(f + i) = shuffle<0xb1>(x8(g + i));
}
}
}
}
using f_n_t_t::predifx2;
using f_n_t_t::difx2;
using f_n_t_t::difx2fxC;
}
namespace field_Z{
namespace calc{
constexpr IL Z _AdMulZ(Z a, Z b, Z c, Z d){
return Space.reduce_s(u64(a) * b + u64(c) * d);
}
IL Zx8 _AdMulZx8(const Zx8 &a, const Zx8 &b, const Zx8 &c, const Zx8 &d){
u32x8 z = (NInvx8 * ((a * b) + (c * d)));
z = blend<0xaa>(RC(u32x8, (fus_mul(a, b) + fus_mul(c, d) + fus_mul(z, modx8)) >> 32), RC(u32x8, (fus_mul(u32x8(u64x4(a) >> 32), u32x8(u64x4(b) >> 32)) + fus_mul(u32x8(u64x4(c) >> 32), u32x8(u64x4(d) >> 32)) + fus_mul(shuffle<0xf5>(z), modx8)))) - modx8;
return z + (RC(Zx8, RC(i32x8, z) < RC(i32x8, zerox8)) & modx8);
}
}
}
namespace poly_base{
void adddot(const Z *A, const Z *B, const Z *C, const Z *D, int n, Z *E){
int i = 0;
for(; i + 7 < n; i += 8){
x8(E + i) = calc::_AdMulZx8(x8(A + i), x8(B + i), x8(C + i), x8(D + i));
}
for(; i < n; ++i){
E[i] = calc::_AdMulZ(A[i], B[i], C[i], D[i]);
}
}
}
namespace poly{
namespace tech{
//Author : QedDust413
void Chirp_Z(const Z *g, int n, Z c, int m, Z *f){
using namespace calc;
pre_aloc mem;
int lim = bit_up((++n) + m), t = std::max<int>(n, m);
Z *F = alocS(lim), *G = alocS(lim), *ipwc = alocS(t + 1), cI = invZ(c);
makeseq::Ci(F, n + m, c);
for(int i = 0; i < n + m; ++i){F[i + 1] = mulZ(F[i], F[i + 1]);}
fl0(F + n + m, F + lim - 1), F[lim - 1] = one_Z;
makeseq::Ci(ipwc, t + 1, cI);
for(int i = 1; i <= t; ++i){ipwc[i] = mulZ(ipwc[i - 1], ipwc[i]);}
for(int i = 1; i < n; ++i){G[i] = mulZ(g[n - i], ipwc[n - i - 1]);}
G[n] = g[0], G[0] = zero_Z, fl0(G + n + 1, G + lim), Conv(F, lim, G), f[0] = F[n - 1];
for(int i = 1; i < m; ++i){f[i] = mulZ(F[i - 1 + n], ipwc[i - 1]);}
}
void translate(const Z *g, int n, Z c, Z *f){
int lim = bit_up(n << 1);
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
dot(g, seqfac(n + 1), n + 1, o), rev(o, n + 1), fl0(o + n + 1, o + lim);
makeseq::Ci(h, n + 1, c), dot(h, n + 1, seqifac(n + 1)), fl0(h + n + 1, h + lim);
Conv(o, lim, h);
Cpy_rev(o, n + 1, f), dot(f, n + 1, seqifac());
}
void multi_eva(const Z *f, int n, const Z *a, int m, Z *o){
using namespace calc;
pre_aloc mem;
int lgn = __qed::lg2u(m), lim = 1 << lgn, lim2 = lim << 1;
predifx2.expand(lim);
Z *G = alocS(lgn * lim2), *g = G, *_d = G;
for(const Z *p = a; p != a + m; ++p, g += 2){
g[0] = subZ(one_Z, *p), g[1] = subZ(Const::neg_one, *p);
}
std::fill(g, _d += lim2, one_Z), g = _d;
for(int t = 2, t2 = 4, M = m & -t; t < lim; t = t2, t2 <<= 1, M &= -t){
Z *lc = g - lim2, *rc = lc + t, *ed = g + (M << 1);
for(; g != ed; g += t2, lc += t2, rc += t2){dot(lc, rc, t, g), difx2fxC(g, t);}
if(M < m){dot(lc, rc, t, g), difx2(g, t), g += t2;}
std::fill(g, _d += lim2, one_Z), g = _d;
}
int Lim = bit_up(n + m);
Z *w = alocS(Lim << 1), *h = w + Lim;
dot(g - lim2, g - lim, lim, w), dit(w, lim);
if(m == lim){w[0] = subZ(w[0], one_Z), w[lim] = one_Z;}
rev(w, m + 1), fl0(w + m + 1, w + Lim), Inv(w, n, h), rev(h, n + 1), fl0(h + n + 1, h + Lim);
Cpy_fl0(f, n + 1, w, Lim), Conv(h, Lim, w), Cpy(h + n, m, w), fl0(w + m, w + lim);
for(int t = lim, t2 = t << 1, M = m & -t, mid = t >> 1; t > 1; t2 = t, t = mid, mid >>= 1, M |= (m & -t)){
Z *lc = (g -= lim2), *rc = lc + t, *p = w, *ed = w + (M << 1);
for(; p != ed; p += t2, lc += t2, rc += t2){
dif(p, t), dot(rc, p, t, p + t), dit(p + t, t), dot(p, t, lc), dit(p, t);
Cpy(p + mid, mid, p + t), Cpy(p + t + mid, mid, p);
}
int rl = m - M - mid;
if(rl > 0){
dif(p, t), dot(rc, p, t, p + t), dit(p + t, t), dot(p, t, lc), dit(p, t);
Cpy(p + mid, rl, p + t), Cpy(p + t + rl, mid, p);
}
}
for(int i = 0; i < m; ++i){
o[i] = w[i << 1];
}
}
void fst_intpol(const Z *x, const Z *y, int n, Z *f){
using namespace calc;
pre_aloc mem;
int lgn = __qed::lg2u(n), lim = 1 << lgn, lim2 = lim << 1;
predifx2.expand(lim);
Z *G = alocS(lgn * lim2), *g = G, *_d = G;
for(const Z *p = x; p != x + n; ++p, g += 2){
g[0] = subZ(one_Z, *p), g[1] = subZ(Const::neg_one, *p);
}
std::fill(g, _d += lim2, one_Z), g = _d;
for(int t = 2, t2 = 4, M = n & -t; t < lim; t = t2, t2 <<= 1, M &= -t){
Z *lc = g - lim2, *rc = lc + t, *ed = g + (M << 1);
for(; g != ed; g += t2, lc += t2, rc += t2){dot(lc, rc, t, g), difx2fxC(g, t);}
if(M < n){dot(lc, rc, t, g), difx2(g, t), g += t2;}
std::fill(g, _d += lim2, one_Z), g = _d;
}
Z *h = alocS(lim2), *w = alocS(lim2), *o = alocS(n + 1);
dot(g - lim2, g - lim, lim, w), dit(w, lim);
if(n == lim){w[0] = subZ(w[0], one_Z), w[lim] = one_Z;}
Cpy_rev(w, n + 1, o), Inv(o, n, h), rev(h, n + 1), fl0(h + n + 1, h + lim2);
deriv(w, n, w), fl0(w + n, w + lim2), Conv(h, lim2, w), Cpy(h + n, n, w);
for(int t = lim, t2 = t << 1, M = n & -t, mid = t >> 1; t > 1; t2 = t, t = mid, mid >>= 1, M |= (n & -t)){
Z *lc = (g -= lim2), *rc = lc + t, *p = w, *ed = w + (M << 1);
for(; p != ed; p += t2, lc += t2, rc += t2){
dif(p, t), dot(rc, p, t, p + t), dit(p + t, t), dot(p, t, lc), dit(p, t);
Cpy(p + mid, mid, p + t), Cpy(p + t + mid, mid, p);
}
int rl = n - M - mid;
if(rl > 0){
dif(p, t), dot(rc, p, t, p + t), dit(p + t, t), dot(p, t, lc), dit(p, t);
Cpy(p + mid, rl, p + t), Cpy(p + t + rl, mid, p);
}
}
for(int i = 0; i < n; ++i){
w[i] = w[i << 1];
}
multi_iv(w, n, o);
for(int i = 0; i < n; ++i){
w[i << 1] = w[i << 1 | 1] = mulZs(y[i], o[i]);
}
fl0(w + (n << 1), w + lim2);
for(int t = 2, t2 = 4; t < lim; t = t2, t2 <<= 1, g += lim2){
int M = (n + (t - 1)) & -t;
Z *lc = g, *rc = lc + t, *p = w, *ed = w + (M << 1);
for(; p != ed; p += t2, lc += t2, rc += t2){
adddot(lc, p + t, rc, p, t, p), difx2<true>(p, t);
}
}
adddot(g, w + lim, g + lim, w, lim, w), dit(w, lim), Cpy(w, n, f);
}
Z linear(const Z *f, const Z *a, int k, i64 n){
auto _cal_linear = [](Z *p, int l, int r){
if(l == 0){
p[0] = Const::_half, l = 1;
}
for(int i = l; i < r; ++i){
p[i] = calc::mulZs(p[i & (i - 1)], f_n_t_t::rt1_I[__qed::crz_32(i) + 2]);
}
};
static pre_base pre_linear(_cal_linear);
int hlim = bit_up(k), lim = hlim << 1;
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim), *q = alocS(lim);
const Z *w = pre_linear(hlim);
predifx2.expand(lim), Cpy_fl0(a, k, h, lim), q[0] = one_Z, Cpy_fl0(f, k, q + 1, lim - 1);
NegP(q, 1, k + 1, q), Conv(h, lim, q), fl0(h + k, h + lim), dif(h, lim);
for(; n; n >>= 1){
f_n_t_t::dif_neg(q, lim, o), dot(q, lim, o), dot(h, lim, o);
if(n & 1){
for(int i = 0; i < hlim; ++i){
h[i] = calc::mulZ(h[i << 1] - h[i << 1 | 1] + mod2, w[i]);
q[i] = q[i << 1];
}
}
else{
for(int i = 0; i < hlim; ++i){
h[i] = calc::mulZ(h[i << 1] + h[i << 1 | 1], Const::_half);
q[i] = q[i << 1];
}
}
n == 1 ? dit(h, hlim) : difx2(h, hlim);
n == 1 ? dit(q, hlim) : difx2(q, hlim);
}
return calc::divZ(h[0], q[0]);
}
}
}
namespace poly_base{
void poi_translate(const Z *p, int n, Z c, int m, Z *o){
int l = n + m, lim = bit_up(l);
pre_aloc mem;
Z *A = alocS(lim), *B = alocS(lim), *fc = alocS(l + 1), *fic = alocS(l + 1);
const Z* ifac = seqifac(n + 1);
for(int i = 0; i <= n; ++i){
Z fx = calc::mulZ(ifac[i], ifac[n - i]);
A[i] = calc::mulZ(p[i], ((n - i) & 1) ? mod2 - fx : fx);
}
fl0(A + n + 1, A + lim);
makeseq::mRisingfac(fc, l + 1, calc::subZ(c, InZs(n)));
multi_iv(fc, l + 1, fic);
for(int i = 0; i < l; ++i){B[i] = calc::mulZ(fc[i], fic[i + 1]);}
fl0(B + l, B + lim), Conv(A, lim, B);
for(int i = 0; i < m; ++i){o[i] = calc::mulZ(A[i + n], fc[i + n + 1]);}
dot(o, m, fic);
}
//Author : killer_joke
namespace EGF{
void dot(Z *f, int lim, const Z *g){
poly_base::dot(f, lim, seqfac(lim));
poly_base::dot(f, lim, g);
}
void dot(const Z *f, const Z *g, int lim, Z *h){
poly_base::dot(f, seqfac(lim), lim, h);
poly_base::dot(h, lim, g);
}
void fromOGF(Z *f, int n){
poly_base::dot(f, n + 1, seqifac(n + 1));
}
void toOGF(Z *f, int n){
poly_base::dot(f, n + 1, seqifac(n + 1));
}
}
}
namespace Falling_Factorial{
using namespace poly_base;
namespace f_d_t{
//Author : killer_joke
Z *EX[mp2 + 1][2];
void prefdt(int lim){
int lg = __qed::lg2d(lim);
if(EX[lg][0]){
return ;
}
EX[lg][0] = alocH(lim << 1), EX[lg][1] = alocH(lim << 1);
Cpy_fl0(seqifac(lim), lim, EX[lg][0], lim << 1), dif(EX[lg][0], lim << 1);
f_n_t_t::dif_neg(EX[lg][0], lim << 1, EX[lg][1]);
}
void fdt(Z *f, int lim){
dif(f, lim << 1), dot(f, lim << 1, EX[__qed::lg2d(lim)][0]), dit(f, lim << 1);
}
void ifdt(Z *f, int lim){
dif(f, lim << 1), dot(f, lim << 1, EX[__qed::lg2d(lim)][1]), dit(f, lim << 1);
}
void _fdtA(Z *p, int l, int r){
for(int i = std::max(1, l); i < r; i <<= 1){
makeseq::E_x(p + i, i), rev(p + i, i);
}
}
void _fdtB(Z *p, int l, int r){
seqiv.expand(r);
for(int i = std::max(2, l); i < r; i <<= 1){
Cpy(seqiv(), i, p + i), dif(p + i, i);
}
}
pre_base dA(_fdtA), dB(_fdtB);
void prefdtx2(int lim){
dA.expand(lim), dB.expand(lim << 1);
}
void fdtx2(const Z *g, int lim, Z *f){
dot(g, dA() + lim, lim, f), fl0(f + lim, lim), dif(f, lim << 1), dot(f, lim << 1, dB() + (lim << 1));
dit(f, lim << 1), dot(f + lim, lim, seqifac()), Cpy(g, lim, f);
}
}
using f_d_t::prefdt;
using f_d_t::fdt;
using f_d_t::ifdt;
void Mul(const Z *f, int n, const Z *g, int m, Z *h){
int lim = bit_up(n + m);
prefdt(lim);
pre_aloc mem;
Z *F = alocS(lim << 1), *G = alocS(lim << 1);
Cpy_fl0(f, n + 1, F, lim << 1), Cpy_fl0(g, m + 1, G, lim << 1);
fdt(F, lim), fdt(G, lim), EGF::dot(F, lim, G);
fl0(F + lim, lim), ifdt(F, lim), Cpy(F, n + m + 1, h);
}
void translate(const Z *g, int n, Z c, Z *f){
int lim = bit_up(n << 1);
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
makeseq::mFallingfac(o, n + 1, c);
EGF::fromOGF(o, n), fl0(o + n + 1, o + lim);
dot(g, seqfac(n + 1), n + 1, h), rev(h, n + 1), fl0(h + n + 1, h + lim);
Conv(o, lim, h), Cpy_rev(o, n + 1, f), dot(f, n + 1, seqifac(n + 1));
}
void fall_down(const Z *g, int n, Z *f){
int lim = bit_up(n);
pre_aloc mem;
Z *o = alocS(lim << 1), *h = alocS(lim), *u = alocS(lim), *r = alocS(lim);
prefdt(lim), f_d_t::prefdtx2(lim);
Cpy_fl0(g, n + 1, o, lim), Cpy(seqi(n + 1), n + 1, u);
for(int i = 0; i <= n; i += 2){
o[i + 1] = calc::addZ(o[i], o[i + 1]);
}
for(int t = 4, mid = 2, t2 = 8; t <= lim; t <<= 1, mid <<= 1, t2 <<= 1){
dot(u, n + 1, u);
for(int j = 0; j <= n; j += t){
f_d_t::fdtx2(o + j, mid, r), f_d_t::fdtx2(o + j + mid, mid, h);
dot(h, u, t, o + j), addP(o + j, t, r);
}
}
fl0(o + lim, lim), ifdt(o, lim), Cpy(o, n + 1, f);
}
void stand_up(const Z *g, int n, Z *f){
int lim = bit_up(n);
pre_aloc mem;
Z *F = alocS(lim), *G = alocS(lim), *u = alocS(lim), *r = alocS(lim);
predifx2.expand(lim);
std::fill(G, G + lim, one_Z), subP(G, lim, seqi(lim));
Cpy_fl0(g, n + 1, F, lim);
for(int t = 2, mid = 1; t <= lim; t <<= 1, mid <<= 1){
for(int j = 0; j <= n; j += t){
Cpy(F + j + mid, mid, u), difx2(u, mid), difx2(F + j, mid);
Cpy(G + j + mid, mid, r), difx2fxC(r, mid), difx2fxC(G + j, mid);
dot(u, t, G + j), addP(F + j, t, u), dot(G + j, t, r);
}
}
dit(F, lim), Cpy(F, n + 1, f);
}
}
namespace Command{
using Stat_Info::report;
void cut_string(){
_Exit(0);
}
}
struct auto_reporter{
~auto_reporter(){
Command::report();
}
}jack;
using poly_base::alocS;
using poly_base::pre_aloc;
using namespace poly_base::Tools;
using namespace field_Z;
void solve(){
int n, k;
std::cin >> n >> k;
pre_aloc mem;
auto F = alocS(n + 1), G = alocS(n + 1), H = alocS(n + 1);
scanP(F, n + 1);
poly::InvSqrt(F, n, G), poly::integ(G, n, H), poly::Exp(H, n, G), poly::subP(F, G, n + 1, H), H[0] = one_Z, poly::Ln(H, n, G), G[0] = one_Z, poly::Pow(G, n, k, H), poly::deriv(H, n, G);
printP(G, n);
}
int main(){
std::cin.tie(nullptr) -> sync_with_stdio(false);
solve();
return 0;
}
/*
不是所有人都能在现实中达成自己的人生目标，岁月会改变许多人。
人的动力不会久留。
虽然一开始总是十分强烈，可这份热情会慢慢退却，最终消磨殆尽。
但是在这模拟的世界里，基于某一时刻的精神状态，可以让他们一时的动机贯穿一生。
*/设做长度为lim的ntt复杂度为E(lim)
Setting 设置模数等需要外部告知的编译期常量
设置 模数(mod), 临时空间限制(sta_l_MB(单位 MB)), 统计信息的详细程度(Detail)
__qed 提供一些辅助函数
Stat_Info 统一管理各类统计信息
ntt_size 进行 ntt 的长度(Detail>0)
fill0_size 填充0的总长度(Detail>1)
copy_size 进行拷贝的总长度(Detail>1)
rev_size 进行翻转的总长度(Detail>1)
max_cost_sta_l 最大临时空间使用量(单位 字节)(Detail>1)
report 向输出流 outf 中报告统计信息
mem_helper 提供关于对齐内存申请的函数
包括动态内存申请 (AlocH) 和临时内存申请 (AlocS)
AlocS 用于申请临时的对齐内存(全局区)
使用 AlocS 前应当先设置 pre_aloc
pre_aloc 是一个预申请器
一个大括号内只能存在一个
当析构的时候会释放所有在这个大括号内申请的临时内存
AlocH 用于申请动态的对齐内存(堆区)
使用 Freed 释放
注 : poly_base内定义了 alocH 和 alocS 用于封装
Montgo 提供常规猛哥马莉约碱的实现
_s 后缀表示返回值 [0, mod) 否则 [0, 2mod]
另外 OutZ 返回值也是 [0, mod)
SIMD 提供了对于指令集命令的一些支持
field_Z 定义了运算所处的域Z 以及常规数字出入域的工具 还有域对指令集的基础支持
InZs 返回值[0, mod)
OutZ 返回值[0, mod)
field_Z::calc 提供了在域下的各类基本运算函数(如 +,-,*,/) 以及域对指令集进行基本运算的支持
大部分的函数都是要求输入 [0, 2mod] 返回值也是
x8一般就是对于向量进行一次8个和普通函数一致的运算
特殊的:
模板函数trans, mulZs(x8)在strict为true时返回值[0, mod)
_AmulZx8要求向量c中的数 [0, mod) 否则返回值不保证
_AdMulZ(x8)要求a,b 中至少有一个 < mod, 且c,d 中至少有一个 < mod 且其返回值[0, mod) 否则返回值不保证
field_Z::calc::extend_ 提供了在域下的各类扩展运算函数(如绝对值)
主要是关于模意义下的平方根
field_Z::Const 定义了域内一些有意义的常数
包括 1 / 2, - 1 / 2, - 1, i, -i, i / 2
poly_base 提供了对于多项式操作的基础
fl0 用于填充0 [bg, ed) / [f, f + l)
Cpy 用于拷贝 [bg, ed) -> [g, g + ed - bg) / [f, f + l) -> [g, g + l)
rev 用于翻转 [bg, ed) / [f, f + l)
Cpy_fl0 / Cpy_rev 上述函数的复合
mul / add / sub _t_s 把一段序列乘以/加上/减去t (支持原地)
dot 点乘 (支持原地)
sub / add P 减 / 加 一个序列 (支持原地)
NegP 把一个序列部分取负并存到另一个序列上 (支持原地)
adddot a b + c d
要求a,b 中至少有一个 < mod, 且c,d 中至少有一个 < mod
则其返回值[0, mod) 否则返回值不保证
multi_iv_Zx8
是一个求取多向量乘法逆元的函数
复杂度 O(l + log mod)
不支持原地 且其中不能含有0
multi_iv
求取多个逆元的函数
复杂度 O(l + log mod)
不支持原地 且其中不能含有0
pre_base
一个封装类, 用于管理各种需要预处理的资源
构造函数传入一个 std::function
用于在 a[l, r) 上构造序列
expand 保证可用长度至少为 len
clear 清除所持有的资源
operator() 提供所有的资源 可传入需要的至少可用长度
poly_base::cal_helper 提供了对于积分和求导的支持
seqi 关于I的连续序列 (0, 1, 2 ...)
seqiv 关于I_iv的连续序列 (0, 1 / 1, 1 / 2 ...)
deriv 求导 (支持原地)
ps:如果对于下降阶乘幂就是差分
integ 积分 (支持原地)
ps:如果对于下降阶乘幂就是和式
poly_base::Tools 提供了对于多项式操作的一些方便的工具
compP 比较两个序列是否一致 返回第一个不同之处或者-1(表示一致)
scanP 使用输入流inf读入一段序列并加上对于域的修正
printP 使用输出流outf输出一段序列并加上对于域的修正
上述函数复杂度均为线性
clzP/crzP 前导0 / 后导0
复杂度为0的个数
Conv 进行长度为lim的卷积, 返回值存f 并且g为其 dif 3E(lim)
poi_translate
p = f(0), f(1), f(2) ... f(n)
o -> f(c), f(c + 1), f(c + 2) ... f(c + n)
也就是连续点值平移
要求 c > n , 可以原地求
6E(n)
poly_base::f_n_t_t_base 提供了关于快速数论变换的工具和信息
base4thre 基4阈值
长度大于阈值时使用基4的 dif / dit
mul_t 把一段序列乘以t 并且包括一个返回值是否严格的模板 不会使用指令集进行优化
mp2 表示 mod - 1 中因子2的数量
_g 表示模数的最小原根
P_R_Tab 2 ^ n 次本原单位根表
omega_info_base2 基2 dif / dit 的单位根信息
powXx8 返回一个由X的[0, 8)次幂组成的向量
ntt_info_base4x8 基4 指令集加速 dif / dit 的信息
Neg 用于把一个向量中的一部分取负
poly_base::f_n_t_t 实现了快速数论变换
dif_base2 基2 dif
dit_base2 基2 dit
dif_base4x8 基4 指令集加速 dif 变换长度至少为8且序列需对32对齐
dit_base4x8 基4 指令集加速 dit 变换长度至少为8且序列需对32对齐
dif 封装的dif 自动选择 dif 方式
函数模板表示返回值是否严格, 还有是否进行域的修正
dit 封装的dit 自动选择 dit 方式
函数模板表示返回值是否严格, 还有是否进行域的修正
predifx2
关于对于difx2预处理的资源
扩展时必须传入2的幂次
预处理后可进行的最大difx2后长度为lim
some E(n) 与之前的长度有关
difx2
f = dif(?, l) -> dif(?, 2l)
2E(lim)
difx2fxC
f = dif(? + x ^ l, l) -> dif(? + x ^ l, 2l)
2E(lim)
dif_neg
dif(f(x)) -> dif(f(-x))
lim > 1
复杂度线性
poly_base::fac_helper 提供了对于阶乘的支持
seqfac 关于 i! 的连续序列 (0!, 1!, 2! ...)
seqifac 关于 1 / i! 的连续序列 (1 / 0!, 1 / 1!, 1 / 2! ...)
poly_base::makeseq 提供了用于生成特定序列的工具
Ci C ^ i -> 1 / (1 - cx)
E_x (-1) ^ i / i! -> e ^ -x
by the way e ^ x -> ifac
mRisingfac m的i次上升幂
mFallingfac m的i次下降幂
以上函数复杂度均为线性
poly 实现了各种多项式(形式幂级数)操作
Mul 进行多项式乘法 3E(n + m)
n, m传入次数, 返回一个 n + m 次的多项式 支持原地
Inv
求 f = 1 / g (mod x ^ {n + 1}) 10E(n)
要求常数项可逆 不支持原地
使用经典的 10E(n) 牛迭求倒数。这个相信大家都会。。
Div(const Z g, int n, const Z f, Z *q)
求 q = f / g (mod x ^ {n + 1}) 10E(n) + eps
要求f可逆 支持原地
使用分块计算除法，可见
https://rogeryoungh.blog.uoj.ac/blog/7530
Div(const Z g, int n, const Z f, int m, Z *q)
求 q 满足 g = f * q + r, q 是一个 n - m 次的多项式 要求f最高次不为0, 10E(n - m) + eps
不支持原地
常规方法 使用 Div
Div(const Z g, int n, const Z f, int m, Z q, Z r)
求 q, r 满足 g = f * q + r, q 是一个 n - m 次的多项式 r是一个低于m次的多项式 要求f最高次不为0
10E(n - m) + 3E(min(2 * m, n)) + eps + O(n)
不支持原地
常规方法 使用 Div
Ln
求 f(x) = Ln(g(x)) mod x ^ {n + 1} 要求g的常数项为1
10E(n) + eps
不支持原地
使用 Div
并基于一个小优化
考虑到求导后再积分的过程中的左右平移可以抵消 因此不实际进行求导和积分中的平移减少内存搬迁导致的低效。
Expi
求 f(x) = Exp(g(x)) mod x ^ {n + 1} (h(x) = 1 / f(x) mod x ^ {n + 1} ) 要求g的常数项为0
如果计算逆 21E(n), 否则 17E(n)
不支持原地
牛顿迭代
方法基于 https://negiizhao.blog.uoj.ac/blog/4671
Exp
求 f(x) = Exp(g(x)) mod x ^ {n + 1} 要求g的常数项为0
14E(n) + eps
不支持原地
使用分块计算Exp，可见
https://rogeryoungh.blog.uoj.ac/blog/7530
并使用Expi加速计算第一块。
InvSqrt
求 f(x), f(x) * f(x) = 1 / g(x) (mod x ^ (n + 1)) 要求g的常数项是非0二次剩余
12E(n)
不支持原地
很简单且经典的方法。
Sqrti
求 f(x), f(x) * f(x) = g(x) (h(x) = 1 / f(x)) (mod x ^ (n + 1)) 要求g的常数项是非0二次剩余
如果计算逆 15E(n), 否则 11E(n)
不支持原地
牛顿迭代
方法基于 https://negiizhao.blog.uoj.ac/blog/4671
Sqrt
求 f(x), f(x) * f(x) = g(x) (mod x ^ {n + 1}) 要求g的常数项是非0二次剩余
8E(n) + eps
使用分块计算Sqrt，可见
https://rogeryoungh.blog.uoj.ac/blog/7530
并使用Sqrti加速计算第一块。
Pow(const Z g, int n, int k, Z f)
求 f(x) = g(x) ^ {k} (mod x ^ {n + 1}) 要求g(x) 常数项为1 k % mod
24E(n) + eps
不支持原地
Pow(const Z g, int n, int k, int k2, Z f)
求 f(x) = g(x) ^ {k} (mod x ^ {n + 1}) 要求g(x) 常数项非0
k % mod, k2 % phi(mod)
24E(n) + eps
不支持原地
Safe_Sqrt
求 f(x), f(x) * f(x) = g(x) (mod x ^ {n + 1})
如果正常返回good_poly 如果多项式是空的返回empty_poly 如果不存在返回not_exist_poly
不支持原地(这叫safe?)
Safe_Pow
求 f(x) = g(x) ^ {k} (mod x ^ {n + 1})
如果正常返回good_poly 如果多项式是空的返回empty_poly
不支持原地(这叫safe?)
poly::Trigonometric_Function 实现了多项式三角函数
Trig 同时求一个多项式的Sin和Cos, 21E(n)
Sin / Cos Trig的封装, 21E(n)
Tan 求一个多项式的Tan, 31E(n) + eps
aSin 求一个多项式的arcsin 22E(n)
aCos 求一个多项式的arccos 22E(n)
aTan 求一个多项式的arctan 14E(n) + eps
均支持原地
照着公式做的
poly::tech 实现了一些多项式科技(多点求值 快速插值 Chirp-Z变换 多项式平移)
Chirp_Z
获取 g(1), g(c^1), g(c^2), g(c^3) .... g(c^(m-1))
3E(n + m)
支持原地
translate
f(x) = g(x + c)
c > n, 6E(n)
支持原地
multi_eva
多点求值
o -> g(a_0), g(a_1) ... g(a_{m - 1})
大概 5 lg n E(n) ?
使用转置算法而不是分治取模。
支持原地
使用了点值维护铺在连续序列上的 product tree 以降低常数
其构建时类似bfs 而非 dfs 保证连续访问性减少了内存访问的常数
并额外使用了一点小技巧
维护的product tree 为正常转置算法实现的转置 把初始向量的每一个翻转一下即可实现
这样可以无需求点值时乘转置 而是直接相乘 大大减少常数。
ps : 点值时对一个多项式进行转置(翻转)亦可行，不过常数略大。
此外由于这样设计之后可以保证最高次为1
因此维护的点值相乘循环卷积溢出后重新复原了一下
额外处理多出来的一项以减少常数。
fst_intpol
快速插值
f -> f(x_0) = y_0, f(x_1) = y_1 ... f(x_{n - 1}) = y_{n - 1}
大概 7 lg n E(n) ?
最后搓原多项式时使用了一份连续序列维护其 减少了一倍的空间常数。
这个过程想要dfs实现也是可以 不需要双倍空间的 怎么实现类似于线段树合并如何线性空间。
但是 bfs 访问连续性上有一定优势。
支持原地
linear
常系数齐次线性递推
传入一个k阶齐次的的递推式f并提供前k项a，计算其第 n 项 (n > k)
大概 2 lg n E(k) ?
方法 :
https://arxiv.org/pdf/2008.08822.pdf
poly_base::EGF 提供了一些关于EGF的工具
dot EGF点乘
fromOGF 从OGF转换为EGF
toOGF 从OGF转换为EGF
以上函数复杂度均为线性
Falling_Factorial 实现了基本下降阶乘幂操作
Falling_Factorial::f_d_t 实现了点值的EGF和FFP的转换
prefdt 预处理长度为 lim 的 fdt / ifdt(注意不是不超过)
具体来说就是把 $e^x$ 和 $e^-x$ 的点值求出来了
并使用其为逆的关系一次 $dif$ 求出两份点值。
2E(lim)
fdt 把一个下降阶乘幂转换为其点值的EGF
[f + lim, f + 2lim) must all zero
大力乘法
4E(lim)
ifdt 把点值的EGF转换为下降阶乘幂
[f + lim, f + 2lim) must all zero
大力乘法
4E(lim)
prefdtx2
预处理x2后长度不超过lim的fdtx2
some E(n) 取决于之前的预处理长度
fdtx2
倍增点值 EGF
fdt(?, lim) -> fdt(?, 2lim)
4E(lim)
预处理的多项式连续点值平移 也可以用 生成函数 推出来一样的式子。
Mul
下降阶乘幂乘法 12E(n + m)
使用了2次fdt和1次ifdt
translate
平移下降阶乘幂
6E(n)
随便用下降幂二项式定理推一推能得到和普通多项式平移差不多的式子。
fall_down
多项式 -> 下降阶乘幂
大概 4 lgn E(n)
stand_up
维护其点值的 EGF
下降阶乘幂 -> 多项式
大概 4 lgn E(n)
Command 用来控制程序的运行状态
report() 用于报告目前的统计信息
再提供一些不太重要的实现细节
真的会有人需要吗?
建议看之前先读一遍，这个不是独立介绍
为了支持C++14写的非常魔怔
其实很多C++20的特性可以优化代码
关于取模的实现
为了性能考虑这个代码使用了大量的惰性取模 因此自己参考时一定要小心
如果为了实现的简洁以及易搬运, 还是建议使用封装取模类
我们假设读者有基础的对猛哥马莉约碱的了解,因为作者显然没有
简而言之就是通过 猛哥马莉约碱 实现快速计算模意义下的 (a * b / R)
因此提前给每个数带一个R 这样 $(aR * bR / R) = abR$
因此需要进入和离开数域 进入数域即为 *= R, 离开即为 /= R
reduce 的返回值 ∈ $[0, 2mod)$ 因此所有的其他计算函数也是以这个为基础的
注意 如果要进行 $(a + b)c$ 不对 $a + b$ 进行取模是有前提的, 必须要求c ∈ $[0, mod)$
否则对于接近 1 << 30 (如998244353) 的模数不能保证约减结果 ∈ $[0, 2mod)$
而所有运算函数都是默认输入 ∈ $[0, 2mod)$，因此可能导致错误
我实现的时候每次 _Amul / _Smul 时保证所乘严格 从而保证了正确性
这两个函数加下划线也是这个考虑 实际上不建议使用这两个函数
dif 时为了减少取模允许中间结果达到 $[0, 4mod)$
然后通过每次取的时候预设结果位于 $[0, 4mod)$ 然后做乘法或者提前处理保证正确性
$dif$ 的返回值一定是正常的。
做严格乘法的代价也不高，所以我实现上是保证了正确性的。
注意: trans函数 的返回值默认 ∈ $[0, mod)$
其他的一些细节
这个实现并不是很吝惜清空以及复制 这两个无论如何比fntt快多了
而且进行一次清空和复制有助于目标区域在缓存里从而提高fntt的效率
实测表明9E(n)的Inv常数不如10E(n)的Inv
所有需要倒数迭代的也同样使用10E(n)的Inv做法
实际上很多函数的复杂度可以更优秀
这里放一下9E(n) Inv 的实现
void Inv_9En(const Z* g, int n, Z *f){
using namespace calc;
int lim = bit_up(n);
pre_aloc mem;
Z *o = alocS(lim), *h = alocS(lim);
f[0]=invZ(g[0]);
constexpr Z img = f_n_t_t::rt1[2], imgI = f_n_t_t::rt1_I[2], img2 = mulZs(f_n_t_t::rt1[2], Const::_half);
for(int t = 2, mid = 1; t <= lim; t <<= 1, mid <<= 1){
int xl = std::min<int>(n + 1, t);
Cpy_fl0(g, xl, o, t), Cpy_fl0(f, mid, h, t);
dif(o, t), dif(h, t), dot(o,t,h), dot(o,t,h), subP(o,t,h), dit(o, t);
Z tr = f_n_t_t::rt1[std::__lg(t)+1];
auto gR = g + mid;
auto hR = h + mid;
for(int i = 0, fx = one_Z; i < mid; ++i, fx = mulZ(fx, tr)){
h[i]=mulZ(f[i],fx), hR[i]=mulZ(fx,i<(xl-mid)?subZ(g[i],mulZ(gR[i],imgI)):g[i]);
}
dif(h, mid), dif(hR, mid), dot(hR,mid,h), dot(hR,mid,h), subP(hR,mid,h), dit(hR, mid);
tr = f_n_t_t::rt1_I[std::__lg(t)+1];
for(int i = mid, fx = one_Z; i < xl; ++i, fx = mulZ(fx, tr)){
f[i] = mulZ(addZ(o[i - mid], mulZ(h[i], fx)) + mulZ(img, o[i]), img2);
}
}
}但是实测这个东西打不过 10E(n) 的 Inv
关于效率的一些测试 可见 https://judge.yosupo.jp/
测试账号名: QedDust413
一些设计中的大坑
1.究竟传入次数还是传入项数
- (a + b) * c 的溢出问题
3.比较两个Z相等不要用 == 请使用equalZ
鄂公网安备 42010202000505 号