多项式全家桶
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) \times B(x) \equiv1\pmod {x^{n}}\)
则有\(A(x) \times B(x) \equiv1\pmod {x^{ \left\lceil\dfrac{n}{2}\right\rceil}}\)
又设\(A(x) \times C(x) \equiv1\pmod {x^{ \left\lceil\dfrac{n}{2}\right\rceil}}\)
那么有\(B(x)-C(x)\equiv0\pmod { x^{\left\lceil\dfrac{n}{2}\right\rceil}}\)
平方可得\(B(x)^2-2\times B(x)\times C(x)-C(x)^2\equiv 0\pmod{x^n}\)
那么有\(A(x)\times B(x)^2-2\times A(x) \times B(x)\times C(x) +A(x)\times C(x)^2\equiv 0\pmod{x^n}\)
\(B(x)-2\times C(x)+A(x)\times C(x)^2\equiv 0\pmod{x^n}\)
\(B(x)=2C(x)-A(x)\times C(x)^2\pmod{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)^2\equiv A(x)\pmod{x^n}\)
有\(C(x)^2-A(x)\equiv 0\pmod{x^n}\)
\(C(x)^4-2\times C(x)^2\times A(x)+A(x)^2\equiv 0\pmod{x^{2n}}\)
两边同时加上\(4\times C(x)^2\times A(x)\)可得
\(C(x)^4+2\times C(x)^2\times A(x)+A(x)^2\equiv 4\times C(x)^2\times A(x)\pmod{x^{2n}}\)
\(\frac{C(x)^4+2\times C(x)^2\times A(x)+A(x)^2}{4C(x)^2}\equiv A(x)\pmod{x^{2n}}\)
那么最后答案为\(B(x)=\frac{C(x)^2+ A(x)}{2\times 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))^{'}\times (g(x)-g(x_0))\equiv 0\pmod{x^{2n}}\)
\(g(x)\equiv g_0(x)\times(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+i\sin x\)
我们同理有\(e^{-ix}=\cos x-i\sin x\)
于是我们相加或者相减即可得到\(\begin{aligned} \sin x&=\frac{e^{ix}-e^{-ix}}{2i}\end{aligned}\) \(\begin{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)\equiv\arcsin F(x)\pmod{x^n}\)
两边求导可得\(G'(x)\equiv\frac{F'(x)}{\sqrt{1-F(x)^2}}\pmod{x^n}\)
然后再积回去\(G(x)\equiv\int\frac{F'(x)}{\sqrt{1-F(x)^2}}\mathrm{d}x\pmod{x^n}\)
对于\(arccos\)同理计算即可最后为\(G(x)\equiv-\int\frac{F'(x)}{\sqrt{1-F(x)^2}}\mathrm{d}x\pmod{x^n}\)
对于\(arctan\)也一样最后为\(G(x)\equiv\int\frac{F'(x)}{1+F(x)^2}\mathrm{d}x\pmod{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\)的性质
\[\sum_{i=0}^{\infty}\frac{i^{\underline{n}}}{i!}x^i=\sum_{i=0}^{\infty}\frac{1}{(i-n)!}x^i=\sum_{i=n}^{\infty}\frac{1}{i!}x^{i-n}=x^n\sum_{i=0}^{\infty}\frac{1}{i!}x^i=x^ne^x
\]
考虑\(A\)的下降幂为\(F(n)=\sum_{i=0}^{\infty}F[i]n^{\underline{i}}\)
那么考虑构造其\(EGF\)
\[\sum_{i=0}^{\infty}\frac{x^i}{i!}\sum_{j=0}^{\infty}F[j]i^{\underline{j}}
\]
\[=\sum_{j=0}^{\infty}F[j]\sum_{i=0}^{\infty}\frac{i^{\underline{j}}}{i!}x^i
\]
\[=\sum_{j=0}^{\infty}F[j]x^je^x
\]
卷积然后卷上\(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\)，那么：
\[\begin{aligned}x\equiv x_1\pmod{A} \\ x\equiv x_2\pmod{B} \\ x\equiv x_3\pmod{C}\end{aligned}
\]
先对前两个用\(CRT\)合并我们有
\[\begin{aligned}x_1+k_1A=x_2+k_2B\\x_1+k_1A\equiv x_2\pmod{B}\\k_1\equiv \frac{x_2-x_1}A\pmod{B}\end{aligned}
\]
于是求出了 \(k_1\)，也就求出了 \(x\equiv x_1+k_1A\pmod{AB}\)，记 \(x_4=x_1+k_1A\)
\[\begin{aligned}x_4+k_4AB=x_3+k_3C\\x_4+k_4AB\equiv x_3\pmod{C}\\k_4\equiv \dfrac{x_3-x_4}{AB}\pmod{C}\end{aligned}
\]
求出了 \(k_4\)， \(x\equiv x_4+k_4AB\pmod{ABC}\)，因为 \(x\lt 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;
}
呐，这份感情，什么时候可以传达呢
