NTT
太多了不写了,直接贴个板子算了
inline void NTT(int *y , int len , int opt){
int *rev = new int[len];
rev[0] = 0;
for(int i = 1;i < len;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
for(int i = 0;i < len;i ++) if(rev[i] > i) swap(y[rev[i]] , y[i]);
for(int i = 1;i < len;i <<= 1){
int G1 = quick_pow(3 , (mod - 1) / (i << 1));
for(int j = 0;j < len;j += (i << 1)){
for(int k = 0 , g = 1;k < i;k ++ , g = 1LL * g * G1 % mod){
int res = 1LL * y[i + j + k] * g % mod;
y[i + j + k] = ((y[j + k] - res) % mod + mod) % mod;
y[j + k] = (y[j + k] + res) % mod;
}
}
}
if(opt == -1){
reverse(y + 1 , y + len);
for(int i = 0 , inv = quick_pow(len , mod - 2);i < len;i ++) y[i] = 1LL * y[i] * inv % mod;
}
delete []rev; rev = 0x0;
}
多项式求逆
(A(x) imes B(x) equiv1pmod {x^{n}})
则有(A(x) imes B(x) equiv1pmod {x^{ leftlceildfrac{n}{2} ight ceil}})
又设(A(x) imes C(x) equiv1pmod {x^{ leftlceildfrac{n}{2} ight ceil}})
那么有(B(x)-C(x)equiv0pmod { x^{leftlceildfrac{n}{2} ight ceil}})
平方可得(B(x)^2-2 imes B(x) imes C(x)-C(x)^2equiv 0pmod{x^n})
那么有(A(x) imes B(x)^2-2 imes A(x) imes B(x) imes C(x) +A(x) imes C(x)^2equiv 0pmod{x^n})
(B(x)-2 imes C(x)+A(x) imes C(x)^2equiv 0pmod{x^n})
(B(x)=2C(x)-A(x) imes C(x)^2pmod{x^n})
inline void poly_inv(int *a , int len){
if(len == 1){ return void(a[0] = quick_pow(a[0] , mod - 2)); }
int len1 = (len + 1) / 2;
int *g = new int[len * 2];
for(int i = 0;i < len1;i ++) g[i] = a[i];
for(int i = len1;i < 2 * len;i ++) g[i] = 0;
poly_inv(g , len1);
for(int i = len1;i < 2 * len;i ++) g[i] = 0;
NTT(g , 2 * len , 1); NTT(a , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) a[i] = ((((1LL * 2 * g[i]) % mod) - 1LL * a[i] * g[i] % mod * g[i] % mod) % mod + mod) % mod;
NTT(a , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) a[i] = 0;
delete []g; g = 0x0;
}
多项式开根
求(C(x)^2equiv A(x)pmod{x^n})
有(C(x)^2-A(x)equiv 0pmod{x^n})
(C(x)^4-2 imes C(x)^2 imes A(x)+A(x)^2equiv 0pmod{x^{2n}})
两边同时加上(4 imes C(x)^2 imes A(x))可得
(C(x)^4+2 imes C(x)^2 imes A(x)+A(x)^2equiv 4 imes C(x)^2 imes A(x)pmod{x^{2n}})
(frac{C(x)^4+2 imes C(x)^2 imes A(x)+A(x)^2}{4C(x)^2}equiv A(x)pmod{x^{2n}})
那么最后答案为(B(x)=frac{C(x)^2+ A(x)}{2 imes C(x)})
inline void poly_sqrt(int *a , int len){
if(len == 1) return;
int len1 = len / 2;
int *f1 = new int[len * 2];
for(int i = 0;i < len1;i ++) f1[i] = a[i];
for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
poly_sqrt(f1 , len1);
for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
int *f2 = new int[len * 2];
for(int i = 0;i < 2 * len;i ++) f2[i] = 1LL * f1[i] * 2 % mod;
poly_inv(f2 , len);
for(int i = len;i < 2 * len;i ++) f2[i] = 0;
NTT(f1 , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) f1[i] = 1LL * f1[i] * f1[i] % mod;
NTT(f1 , 2 * len , -1);
for(int i = 0;i < len;i ++) a[i] = (a[i] + f1[i]) % mod;
NTT(a , 2 * len , 1); NTT(f2 , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f2[i] % mod;
NTT(a , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) a[i] = 0;
delete []f1; f1 = 0x0;
delete []f2; f2 = 0x0;
}
多项式求(mathcal{ln})
考虑一个多项式的导数(g(x)^{'}=frac{f(x)^{'}}{f(x)})
将(g(x))积分回去可得(int{g(x)}=mathcal{ln}f(x))
inline void poly_dev(int *a , int len){
for(int i = 1;i < len;i ++) a[i - 1] = 1LL * a[i] * i % mod;
a[len - 1] = 0;
}
inline void poly_inv_dev(int *a , int len){
for(int i = len + 1 ; i ; i --) a[i] = 1LL * a[i - 1] * quick_pow(i , mod - 2) % mod;
a[0] = 0;
}
inline void poly_ln(int *a , int len){
int *f = new int[2 * len];
for(int i = 0;i < len;i ++) f[i] = a[i];
for(int i = len;i < len * 2;i ++) f[i] = 0;
poly_dev(f , len); poly_inv(a , len);
NTT(f , 2 * len , 1); NTT(a , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * f[i] * a[i] % mod;
NTT(a , 2 * len , -1);
poly_inv_dev(a , len);
for(int i = len;i < 2 * len;i ++) a[i] = 0;
delete []f; f = 0x0;
}
多项式求(mathcal{exp})
(g(x)equiv e^{f(x)}pmod{x^n})
同时取(mathcal{ln})可得(lng(x)-f(x)=0)
则有(h(g(x))=ln(g(x))-f(x))
(h(g(x))^{'}=frac{1}{g(x)})
(h(g(x))=frac{h(g_0(x))}{0!}+frac{h(g_0(x)^{'})}{1!}+......)
(h(g_0(x))+h(g_0(x))^{'} imes (g(x)-g(x_0))equiv 0pmod{x^{2n}})
(g(x)equiv g_0(x) imes(f(x)-lng_0(x))+g_0(x)pmod{x^{2n}})
inline void poly_exp(int *a , int len){
if(len == 1) return void(a[0] ++);
int *f = new int[len * 2]; int len1 = len / 2;
for(int i = 0;i < len1;i ++) f[i] = a[i];
for(int i = len1;i < len;i ++) f[i] = 0;
poly_exp(f , len1);
for(int i = len1;i < len * 2;i ++) f[i] = 0;
int *g = new int[len * 2];
for(int i = 0;i < len * 2;i ++) g[i] = f[i];
poly_ln(g , len); a[0] ++;
for(int i = 0;i < len;i ++) a[i] = ((a[i] - g[i]) % mod + mod) % mod;
NTT(a , len * 2 , 1); NTT(f , len * 2 , 1);
for(int i = 0;i < len * 2;i ++) a[i] = 1LL * a[i] * f[i] % mod;
NTT(a , len * 2 , -1);
for(int i = len;i < len * 2;i ++) a[i] = 0;
delete []g; g = 0x0;
delete []f; f = 0x0;
}
多项式求三角函数
根据欧拉公式(e^{ix}=cos x+isin x)
我们同理有(e^{-ix}=cos x-isin x)
于是我们相加或者相减即可得到(egin{aligned} sin x&=frac{e^{ix}-e^{-ix}}{2i}end{aligned}) (egin{aligned}cos x&=frac{e^{ix}+e^{-ix}}{2} end{aligned})
直接上(exp)即可,(i)为在(w_4),在膜(998244353)下为(g^{frac{p-1}{4}})
inline void poly_sin(int *a , int len){
static int f[kato] , g[kato];
for(int i = 0;i < len;i ++) f[i] = 1LL * a[i] * img % mod;
poly_exp(f , len);
for(int i = 0;i < len;i ++) g[i] = f[i];
poly_inv(g , len);
for(int i = 0 , inv = quick_pow(2 * img , mod - 2);i < len;i ++) a[i] = 1LL * inv * (((f[i] - g[i]) % mod + mod) % mod) % mod;
}
inline void poly_cos(int *a , int len){
static int f[kato] , g[kato];
for(int i = 0;i < len;i ++) f[i] = 1LL * a[i] * img % mod;
poly_exp(f , len);
for(int i = 0;i < len;i ++) g[i] = f[i];
poly_inv(g , len);
for(int i = 0;i < len;i ++) a[i] = 1LL * inv2 * ((f[i] + g[i]) % mod) % mod;
}
多项式求反三角函数
对于(arcsin)我们有(G(x)equivarcsin F(x)pmod{x^n})
两边求导可得(G'(x)equivfrac{F'(x)}{sqrt{1-F(x)^2}}pmod{x^n})
然后再积回去(G(x)equivintfrac{F'(x)}{sqrt{1-F(x)^2}}mathrm{d}xpmod{x^n})
对于(arccos)同理计算即可最后为(G(x)equiv-intfrac{F'(x)}{sqrt{1-F(x)^2}}mathrm{d}xpmod{x^n})
对于(arctan)也一样最后为(G(x)equivintfrac{F'(x)}{1+F(x)^2}mathrm{d}xpmod{x^n})
求导+开根+求逆+积分一套理论
inline void poly_asin(int *a , int len){
int *f = new int[kato];
for(int i = 0;i < len;i ++) f[i] = a[i];
for(int i = len;i < 2 * len;i ++) f[i] = 0;
NTT(f , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) f[i] = 1LL * f[i] * f[i] % mod;
NTT(f , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) f[i] = 0;
for(int i = 0;i < len;i ++) f[i] = (mod - f[i]) % mod;
f[0] = 1;
poly_sqrt(f , len); poly_inv(f , len); poly_dev(a , len);
NTT(a , 2 * len , 1); NTT(f , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f[i] % mod;
NTT(a , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) a[i] = 0;
poly_dev_inv(a , len);
delete []f; f = 0x0;
}
inline void poly_atan(int *a , int len){
int *f = new int[kato];
for(int i = 0;i < len;i ++) f[i] = a[i];
for(int i = len;i < 2 * len;i ++) f[i] = 0;
NTT(f , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) f[i] = 1LL * f[i] * f[i] % mod;
NTT(f , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) f[i] = 0;
f[0] = 1;
poly_inv(f , len); poly_dev(a , len);
NTT(a , 2 * len , 1); NTT(f , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f[i] % mod;
NTT(a , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) a[i] = 0;
poly_dev_inv(a , len);
delete []f; f = 0x0;
}
下降幂多项式乘法
考虑构造关于(A)和(B)的(EGF)
先来点(EGF)的性质
考虑(A)的下降幂为(F(n)=sum_{i=0}^{infty}F[i]n^{underline{i}})
那么考虑构造其(EGF)
卷积然后卷上(e^{-x})即可,最后记得求出来的是(EGF),乘上阶乘
for(int i = 0;i <= yuni;i ++){
A[i] = fac_inv[i];
B[i] = (i & 1) ? mod - fac_inv[i] : fac_inv[i];
}
for(len = 1;len <= 2 * yuni;len <<= 1);
NTT(a , len , 1); NTT(A , len , 1);
for(int i = 0;i < len;i ++) a[i] = 1LL * a[i] * A[i] % mod;
NTT(a , len , -1);
NTT(b , len , 1); NTT(B , len , 1);
for(int i = 0;i < len;i ++) b[i] = 1LL * b[i] * A[i] % mod;
NTT(b , len , -1);
for(int i = 0;i <= yuni;i ++) ans[i] = 1LL * a[i] * b[i] % mod * fac[i] % mod;
NTT(ans , len , 1);
for(int i = 0;i < len;i ++) ans[i] = 1LL * ans[i] * B[i] % mod;
NTT(ans , len , -1);
for(int i = 0;i <= n + m;i ++) printf("%d " , ans[i]);
任意模数多项式乘法
假设这一位的答案是 (x),三个模数分别为 (A,B,C),那么:
先对前两个用(CRT)合并我们有
于是求出了 (k_1),也就求出了 (xequiv x_1+k_1Apmod{AB}),记 (x_4=x_1+k_1A)
求出了 (k_4), (xequiv x_4+k_4ABpmod{ABC}),因为 (xlt ABC),所以 (x=x_4+k_4AB)
const int mod1 = 998244353 , mod2 = 1004535809 , mod3 = 469762049;
const LL mod_1_2 = 1LL * mod1 * mod2;
const int inv1 = quick_pow(mod1 , mod2 - 2 , mod2) , inv2 = quick_pow(mod_1_2 % mod3 , mod3 - 2 , mod3);
struct Int{
int A , B , C;
Int(int A = 0 , int B = 0 , int C = 0): A(A) , B(B) , C(C){}
inline void init(int _n){ A = _n; B = _n; C = _n; }
inline friend Int operator+(const Int &a , const Int &b){
return Int(add(a.A , b.A , mod1) , add(a.B , b.B , mod2) , add(a.C , b.C , mod3));
}
inline friend Int operator-(const Int &a , const Int &b){
return Int(del(a.A , b.A , mod1) , del(a.B , b.B , mod2) , del(a.C , b.C , mod3));
}
inline friend Int operator*(const Int &a , const Int &b){
return Int(1LL * a.A * b.A % mod1 , 1LL * a.B * b.B % mod2 , 1LL * a.C * b.C % mod3);
}
inline int get_ans(){
LL x = static_cast<LL>(B - A + mod2) % mod2 * inv1 % mod2 * mod1 + A;
return (static_cast<LL>(C - x % mod3 + mod3) % mod3 * inv2 % mod3 * (mod_1_2 % mod) % mod + x) % mod;
}
};
二次剩余多项式开根
struct Complex{
int x , y;
Complex(int x = 0 , int y = 0): x(x) , y(y) { }
friend bool operator ==(const Complex &a , const Complex &b){
return a.x == b.x && a.y == b.y;
}
friend Complex operator *(const Complex &a , const Complex &b){
return Complex((static_cast<LL>(a.x) * b.x % mod + static_cast<LL>(a.y) * b.y % mod * qaq % mod) % mod , (static_cast<LL>(a.x) * b.y % mod + static_cast<LL>(a.y) * b.x % mod) % mod);
}
};
inline Complex quick_pow(Complex a , int b){
Complex res = Complex(1 , 0);
for(; b ; b >>= 1 , a = a * a){
if(b & 1){
res = res * a;
}
}
return res;
}
inline bool judge(int x){
return quick_pow(x , (mod - 1) >> 1) == 1;
}
inline bool Cipolla(int n , int p , int &ans0 , int &ans1){
if(quick_pow(n , (mod - 1) / 2) == mod - 1) return 0;
int x = rand() % mod;
while(!x || judge((static_cast<LL>(x) * x % mod + mod - n) % mod)) x = rand() % mod;
qaq = (static_cast<LL>(x) * x % mod + mod - n) % mod;
ans0 = quick_pow(Complex(x , 1) , (mod + 1) >> 1).x;
ans1 = mod - ans0;
if(ans0 > ans1) swap(ans0 , ans1);
return 1;
}
inline void poly_sqrt(int *a , int len){
if(len == 1){ Cipolla(a[0] , mod , ans0 , ans1) , a[0] = min(ans0 , ans1); return; }
int len1 = len / 2;
int *f1 = new int[len * 2];
for(int i = 0;i < len1;i ++) f1[i] = a[i];
for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
poly_sqrt(f1 , len1);
for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
int *f2 = new int[len * 2];
for(int i = 0;i < 2 * len;i ++) f2[i] = 1LL * f1[i] * 2 % mod;
poly_inv(f2 , len);
for(int i = len;i < 2 * len;i ++) f2[i] = 0;
NTT(f1 , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) f1[i] = 1LL * f1[i] * f1[i] % mod;
NTT(f1 , 2 * len , -1);
for(int i = 0;i < len;i ++) a[i] = (a[i] + f1[i]) % mod;
NTT(a , 2 * len , 1); NTT(f2 , 2 * len , 1);
for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f2[i] % mod;
NTT(a , 2 * len , -1);
for(int i = len;i < 2 * len;i ++) a[i] = 0;
delete []f1; f1 = 0x0;
delete []f2; f2 = 0x0;
}