|
|
作者：中二攻子
链接：https://ac.nowcoder.com/discuss/175409
来源：牛客网
本文含NTT、MTT、拆系数FFT、共轭优化FFT、多项式求逆与ln
约定：
1、 F(x) 表示一个普通的项数为 2 的幂次多项式， F_D(x) 是他的点值表示。
2、 w 代表单位根， w_m 表示 m 次单位根。
3、 A 代表一个数列。
4、 g 表示原根。
多项式系列之零 底层知识：
多项式的表示：
多项式可以通过系数数列 A 表示， a_i 是 x_i 的系数。
多项式可以通过点值表示，对于一个 n 次多项式，取 n 种不同的
取值带入 F(x) ，得到 n 个值，在取相同这 n 个数的意义下，可以唯一的表示这个多项式。
多项式乘法：
定义 F(x)\oplus G(x)=\sum_{i=0}^n\sum_{j=0}^i f_ix^i\times g_{j-i}x^{j-i} ，在系数表示之下相乘复杂度 \Theta(n^2) ，在点值表示之下 F(x)\oplus G(x)=A_f\times A_g=\sum_{i=1}^n a_{fi}\times a_{gi} ，复杂度 \Theta(n) 。
复数：
复数一般情况下可以表示成 a+bi 的形式， a,b 是实数， i=\sqrt{-1} 。
复数的幅角：平面直角坐标系上点 (a,b) 所在的任意角。
复数的模长： \sqrt{a^2+b^2}
两个复数相乘： (a+bi)\times(c+di)=ac+adi+bci-bd=(ac-bd)+(ad+bc)i ，复数相乘之后，模长等于原来两个复数的模长的乘积，幅角的角度等于原来两个幅角的和。
复数可以加减乘除，可以和实数一样的带入 F(x) 。
w ：
在单位圆上从 w_m^0=(1,0) 开始平均取 m 个点，从 0 开始编号，分别是 w_m^0,w_m^1,w_m^2,w_m^3\cdots w_m^{m-1} 。
画图观察可得：
1.w_m^k=(cos(\frac k m 2\pi),sin(\frac k m 2\pi)) 所代表的复数
2.w_m^{-k}=(cos(\frac {-k} n 2\pi),sin(\frac {-k} n 2\pi)) 所代表的复数
3.w_m^m=w_m^0=(1,0)
4.w_m^m=-w_m^{\frac m 2}
5.w_{2m}^{2k}=w_m^k
6.w_m^{k+\frac m 2}=-w_m^k
DFT&IDFT：
在科学的数学函数意义上DFT是讲一个函数转化成三角函数的加减乘除的形式，三角函数的系数是原函数系数与点值之间的变换规律。IDFT是DFT的逆变换。
g ：
1、什么是 g ：在 \mod~p 意义下 g^i(i\in[0,p-1]) 互不相同，即 g 可以张成整个 \mod~p 下的域。
2、 g 存在的条件： p=2,4,q^a,2q^a ， q 是奇素数。
3、如何求 g ：把 \phi(p) 进行质因数分解 \phi(p)=\prod p_i^{a_i} ，如果对于任意的 p_i ，总有 g^{\frac {\phi(p)} {p_i}}\neq 1(mod~p) ，暴力枚举即可。
CRT合并：
求解 {\begin{cases}x\equiv a_1 (\mod~p_1)\\x\equiv a_2 (\mod~p_2)\end{cases}},\gcd (p_1,p_2)=1
由 x\equiv a_1 (mod~p_1) ，得
x-p_1y=a_1
x=a_1+p_1y
带入二式，得
a_1+p_1y\equiv a_2(mod~p_2)
p_1y\equiv a_2-a_1(mod~p_2)
若 \gcd (p1,p2)==1 ，用逆元直接除便可；否则通过 exgcd 可求得 y ，若无解则方程组无解。
最后 x=p_1y+a_1(mod~p_1p_2) 。
多项式全集之一 FFT：
什么是FFT：
FFT是利用DFT的特殊性质，把 w 带入 x 从而 \Theta(nlogn) 求一个系数多项式的点值表示，所以叫FDFT。
w 的具体应用：
1、可以方便的IDFT：
设 F(x) 的系数是 A ，在 w_m 的DFT下点值是 B ， G(x) 的系数是 B ，在 w_m^{-1} 的DFT下点值是 C 。
c_k=\sum_{i=0}^{m-1}b_iw_m^{-ki}
c_k=\sum_{i=0}^{m-1}(\sum_{j=0}^{m-1}a_jw_m^{ji})w_m^{-ki}
c_k=\sum_{j=0}^{m-1}a_j\sum_{i=0}^{m-1}w_m^{(j-k)i}
当 j-k=0 时 \sum_{i=0}^{m-1}w_m^{(j-k)i}=m ，否则根据等比数列求和公式得
\frac{w_m^{(j-k)m}-1}{w_m^{j-k}-1}=\frac{w_m^{m(j-k)}-1}{w_m^{j-k}-1}=\frac{1-1}{w_m^{j-k}-1}=0
由此可得： c_k=m\times a_k ， a_k=\frac{c_k}{m}
综上所述，对于点值取的 w 相反数做DFT再除以 m 可得到系数。
2、可以快速的DFT：
直接将 w 带入多项式做DFT需要复杂度 \Theta(n^2) ，我们利用 m 的性质优化：
把 F(x) 按照奇偶分裂， F(x)=(a_0x^0+a_2x^2+a_4x^4+a_6x^6\cdots)+(a_1x^1+a_3x^3+a_5x^5+a_7x^7\cdots)
我们令 F_0(x)=a_0x^0+a_2x^1+a_4x^2+a_6x^3\cdots+a_{m-2}x^{\frac m 2}
F_1(x)=a_1x^0+a_3x^1+a_5x^2+a_6x^3\cdots+a_{m-1}x^{\frac m 2-1}
我们可以发现 F(x)=F_0(x^2)+xF_1(x^2)
现在我们把 w_m 带入，令 k<\frac m 2 ：
F(w_m^k)=F_0(w_m^{2k})+w_m^kF_1(w_m^{2k})
F(w_m^k)=F_0(w_{\frac m 2}^{k})+w_m^kF_1(w_{\frac m 2}^{k})
F(w_m^{k+\frac m 2})=F_0(w_m^{2k+m})+w_m^{k+\frac m 2}F_1(w_m^{2k+m})
F(w_m^{k+\frac m 2})=F_0(w_{\frac m 2}^{k} )-w_m^{k}F_1(w_{\frac m 2}^{k})
我们知道取 w_{\frac m 2} 时， A_0,A_1 的取值，就可以算出 A(w_m) ，而 A_0,A_1 的长度都为 A 的一半，所以可以递归计算。
非递归优化FFT：
1、优化原理：
画图可知，递归版FFT最底层结束状态第 i 个位置的项是 i 二进制翻转后的结果。我们可以 \Theta(n) 的得到最底层的结果，然后向上模拟回溯合并即可。
2、蝴蝶变换：
由上述式子： F(w_m^k)=F_0(w_{\frac m 2}^{k})+w_m^kF_1(w_{\frac m 2}^{k})
可得在迭代时 w_m^k 与 w_m^{k+\frac m 2} 都只与 w_{\frac m 2}^k,w_{\frac m 2}^{k+\frac m 2} 有关，所以我们可以用临时变量记录下一层的两个信息向上迭代。
共轭复数优化FFT：
1、优化原理：
（在DFT时）
令 D(x)=A(x)+iB(x)
E(x)=A(x)-iB(x)
那么
E_D(x)=conj(D_D(m-x))
A_D(x)=\frac{D_D(x)+E_D(x)} 2=\frac{D_D(x)+conj(D_D(m-x))} 2
B_D(x)=\frac{D_D(x)-E_D(x)} {2i}=-i\frac{D_D(x)-conj(D_D(m-x))} 2
2、证明：
step~1 ：
D_D(w_m^k)=A_D(w_m^k)+iB_D(w_m^k)
D_D(w_m^k)=\sum_{j=0}^{m-1}a_jw_m^{jk}+ib_jw_m^{jk}
D_D(w_m^k)=\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{jk}
step~2 ：为方便起见，我们用 x 代替 \frac{2\pi j k} {m}
\begin{equation} E_D(w_m^k)=A_D(w_m^k)-iB_D(w_m^k)\\ E_D(w_m^k)=\sum_{j=0}^{m-1}a_jw_m^{jk}-ib_jw_m^{jk}\\ E_D(w_m^k)=\sum_{j=0}^{m-1}(a_j-ib_j)w_m^{jk}\\ E_D(w_m^k)=\sum_{j=0}^{m-1}(a_j-ib_j)(cosX+isinX)\\ E_D(w_m^k)=\sum_{j=0}^{m-1}(a_jcosX+b_jsinX)+i(a_jsinX-b_jcosX)\\ E_D(w_m^k)=conj\big(\sum_{j=0}^{m-1}(a_jcosX+b_jsinX)-i(a_jsinX-b_jcosX)\big)\\ E_D(w_m^k)=conj\big(\sum_{j=0}^{m-1}(a_jcosX+b_jsinX)-i(a_jsinX-b_jcosX)\big)\\ E_D(w_m^k)=conj\big(\sum_{j=0}^{m-1}(a_jcos(-X)-b_jsin(-X))+i(a_jsin(-X)+b_jcos(-X))\big)\\ E_D(w_m^k)=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)(cos(-X)+isin(-X)))\big)\\ E_D(w_m^k)=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{-jk}\big)\\ E_D(w_m^k)=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{jm-jk}\big)=conj(D_D(w_m^{m-k}))\\ \end{equation}
而在IDFT时，我们需要
D_D(w_m^{-k})=\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{-jk}
E_D(w_m^{-k})=conj\big(\sum_{j=0}^{m-1}(a_j+ib_j)w_m^{jk}\big)=conj(D_D(w_m^{k}))
数论优化FFT(NTT)：
g^{\frac {p-1} m} 与 w_m 的共性：
1、 (g^{\frac {p-1} m})^k 和 w_m^k 都互不相同 (k\in[0,m-1]) 。
2、 (g^{\frac {p-1} m})^m=g^{p-1}=1 ， w_m^m=1 。
3、 (g^{\frac {p-1} m})^{\frac m 2}=g^{\frac{p-1} 2}=\sqrt{g^{p-1}} ，由于原根的互不相同， (g^{\frac {p-1} m})^{\frac m 2}=-1=-g^{p-1}=-(g^{\frac {p-1} m})^m ， w_m^m=-w_m^{\frac m 2} 。
4、 (g^{\frac {p-1} {2m}})^{2k}=(g^{\frac {p-1} {m}})^{k} ， w_{2m}^{2k}=w_m^k
5、 (g^{\frac {p-1} {m}})^{k+\frac m 2}=(g^{\frac {p-1} {m}})^{k}\times (g^{\frac {p-1} {m}})^{\frac m 2}=-(g^{\frac {p-1} {m}})^{k} ， w_m^{k+\frac m 2}=-w_m^k
因为有这些共性，所以 g^{\frac {p-1} m} 可以代替 w_m 。
喜闻乐见的模板：
FFT模板（共轭优化）
namespace FFT{
const double pi = acos(-1);
struct cp{
double x, y;
cp() {x = y = 0;}
cp(double X,double Y) {x = X; y = Y; }
cp conj() {return (cp) {x, -y};}
}a[3000005], b[3000005], c[3000005], I(0, 1), d[3000005];
cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
cp operator* (const cp &a, double b) {return (cp){a.x * b, a.y * b};}
cp operator/ (const cp &a, double b) {return (cp){a.x / b, a.y / b};}
struct p_l_e{
int wz[3000005];
void FFT(cp *a, int N, int op) {
for(int i = 0; i < N; i++)
if (i<wz) swap(a,a[wz]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
for(int i = 0; i < N; i += le) {
cp x, y, w = (cp) {1, 0};
cp wn = (cp){cos(op * pi / mid), sin(op * pi / mid)};
for(int j = 0 ; j < mid; j++) {
x = a[i + j]; y = a[i + j + mid] * w;
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void D_FFT(cp *a, cp *b, int N, int op){
for(int i = 0; i < N; i++) d = a + I * b;
FFT(d, N, op);
d[N] = d[0];
for(int i = 0; i < N; i++){
a = (d + d[N - i].conj()) / 2;
b = I * (-1) * (d - d[N - i].conj()) / 2;
}
d[N] = cp(0, 0);
}
void mult(cp *a, cp *b, cp *c, int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
D_FFT(a, b, N, 1);
for(int i = 0; i < N; i++) c = a * b;
FFT(c, N, -1);
for(int i = 0; i < N; i++) c.x = c.x / N;
}
}PLE;
int n, m;
void main() {
scanf(&#34;%d%d&#34;, &n, &m); n++; m++;
for(int i = 0; i < n; i++) scanf(&#34;%lf&#34;, &a.x);
for(int i = 0; i < m; i++) scanf(&#34;%lf&#34;, &b.x);
PLE.mult(a, b, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf(&#34;%d &#34;, (int)round(c.x));
return;
}
}
NTT模板：
namespace NTT{
typedef long long LL;
const int mod = 998244353;
const int g = 3;
LL a[3000005], b[3000005], c[3000005];
int n, m;
LL qpow(LL a, LL b){
LL ans = 1;
while(b){
if(b & 1) ans = ans * a % mod;
a = a * a % mod;
b >>= 1;
}
return ans;
}
struct p_l_e{
int wz[3000005];
void NTT(LL *a, int N, int op) {
for(int i = 0; i < N; i++)
if(i < wz) swap(a, a[wz]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
LL wn = qpow(g, (mod - 1) / le);
if(op == -1) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += le) {
int w = 1, x, y;
for(int j = 0; j < mid; j++) {
x = a[i + j];
y = a[i + j + mid] * w % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
void mult(LL *a, LL *b, LL *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
NTT(a, N, 1); NTT(b, N, 1);
for(int i = 0; i < N; i++) c = a * b % mod;
NTT(c, N, -1);
LL t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) c = c * t % mod;
}
}PLE;
void main() {
scanf(&#34;%d%d&#34;, &n, &m); n++; m++;
for(int i = 0; i < n; i++) scanf(&#34;%lld&#34;, &a);
for(int i = 0; i < m; i++) scanf(&#34;%lld&#34;, &b);
PLE.mult(a, b, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf(&#34;%lld &#34;, c);
}
}
多项式全集之二 任长任模的FFT:
三模NTT实现任模FFT：
1、为什么要用MTT：当 p 不是NTT模数或者多项式长度大于模数限制时，就要使用MTT。
2、MTT的使用原理：我们对初始多项式取模，那么如果在不取模卷积情况下，答案 x 不会超过 N\times p^2 。我们取三个NTT模数 p_1,p_2,p_3 ，分别做多项式乘法，得到 x 分别 \mod~p_1,p_2,p_3 的答案，通过CRT合并可以得到 x~\mod~p_1p_2p_3 的答案，如果 x<p_1p_2p_3 那么就可以得到准确的答案，再对 p 取模即可。
3、CRT合并的小优化：
step~0: 初始式子
{\begin{cases}x\equiv c_1(mod~p_1)\\x\equiv c_2(mod~p_2)\\x\equiv c_3(mod~p_3)\end{cases}}
step~1: 把一式二式合并(LL范围内)。
{\begin{cases}x\equiv a(mod~p_1p_2)\\x\equiv c_3(mod~p_3)\end{cases}}
step~2: 再次合并(不需要 long~double 快速乘)。
4、常用NTT模数：
以下模数的共同 g=3189
拆系数FFT(CFFT)实现任模FFT：
1、实现原理：运用实数FFT不取模做乘法，然后取模回归到整数。但是由于误差较大（值域是
），我们令 t=\sqrt{m} 把系数 a_i=k_it+b_i ，对 k_i,t_i 交叉做四遍卷积，求出答案按系数贡献取模加入。
2、可按合并DFT的方法优化DFT次数。
算法实现任长FFT：
当 m 不是 2 的幂次的时候，我们从式子入手：
A_k=\sum_{j=0}^{m-1}a_jw_m^{jk}
A_k=\sum_{j=0}^{m-1}a_jw_m^{\frac{j^2+k^2-{(k-j)}^2}{2}}
A_k=w_m^{\frac {k^2} 2}\sum_{j=0}^{m-1}a_jw_m^{\frac{j^2} 2}w_m^{\frac{-{(k-j)}^2}{2}}
令 X_i=a_iw_m^{\frac {i^2} 2},Y_i=w_m^{\frac{-i^2}2}
A_k=w_m^{\frac {k^2} 2}\sum_{j=0}^{m-1}X_jY_{k-j}
喜闻乐见的模板：
三模NTT模板（注意：不可以MTT回来，因为系数会取模）
namespace MTT{
typedef long long LL;
int n, m;
LL p, mod;
const LL p1 = 998244353;
const LL p2 = 1004535809;
const LL p3 = 104857601;
const int g = 3189;
LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005];
LL c3[300005], c1[300005], c2[300005];
LL qpow(LL a, LL b, LL mod) {
LL ans = 1;
while(b) {
if(b & 1) ans = ans * a % mod;
a = a * a % mod;
b >>= 1;
}
return ans;
}
const LL inv12 = qpow(p1, p2 - 2, p2);
const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3);
struct p_l_e{
int wz[300005];
void MTT(LL *a, int N, int op) {
for(int i = 0; i < N; i++)
if(i < wz) swap(a, a[wz]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
LL wn = qpow(g, (mod - 1) / le, mod);
if(op == -1) wn = qpow(wn, mod - 2, mod);
for(int i = 0; i < N ;i += le) {
LL w = 1, x, y;
for(int j = 0; j < mid; j++) {
x = a[i + j];
y = a[i + j + mid] * w % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
void mult(LL *a, LL *b, LL *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
MTT(a, N, 1); MTT(b, N, 1);
for(int i = 0; i < N; i++) c = a * b % mod;
MTT(c, N, -1);
LL t = qpow(N, mod - 2, mod);
for(int i = 0; i < N; i++) c = c * t % mod;
}
}PLE;
LL CRT(LL c1, LL c2, LL c3) {
LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2));
LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p;
return y;
}
void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) {
for(int i = 0; i < N; i++)
c = CRT(c1, c2, c3);
return;
}
void main() {
scanf(&#34;%d%d%lld&#34;, &n, &m, &p); n++; m++;
for(int i = 0; i < n; i++) scanf(&#34;%lld&#34;, &a);
for(int i = 0; i < m; i++) scanf(&#34;%lld&#34;, &b);
mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1);
mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1);
mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1);
merge(c1, c2, c3, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf(&#34;%lld &#34;, (c % p + p) % p);
return;
}
}
拆系数FFT模板（注意：相同系数的两项可以合并一起IDFT。采用共轭优化法，只进行四次DFT）
namespace CFFT{
typedef long long LL;
int n, m, p ,sqrp;
int a[300005], b[300005];
const long double pi = acos(-1);
struct cp{
long double x, y;
cp() {x = y = 0;}
cp(long double X,long double Y) {x = X; y = Y; }
cp conj() {return (cp) {x, -y};}
}ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005];
cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};}
cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};}
struct p_l_e{
int wz[300005];
void FFT(cp *a, int N, int op){
for(int i = 0; i < N; i++)
if(i < wz) swap(a, a[wz]);
for(int le = 2; le <= N; le <<= 1){
int mid = le >> 1;
cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)};
for(int i = 0; i < N; i += le){
w = (cp){1, 0};
for(int j = 0; j < mid; j++){
x = a[i + j];
y = a[i + j + mid] * w;
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void D_FFT(cp *a, cp *b, int N, int op){
for(int i = 0; i < N; i++) d = a + I * b;
FFT(d, N, op);
d[N] = d[0];
if(op == 1){
for(int i = 0; i < N; i++){
a = (d + d[N - i].conj()) / 2;
b = I * (-1) * (d - d[N - i].conj()) / 2;
}
} else {
for(int i = 0; i < N; i++){
a = cp(d.x, 0);
b = cp(d.y, 0);
}
}
d[N] = cp(0, 0);
}
void mult(int *a, int *b, int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
for(int i = 0; i < N; i++){
ka.x = a >> 15;
kb.x = b >> 15;
ta.x = a & 32767;
tb.x = b & 32767;
}
D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1);
for(int i = 0; i < N; i++){
kk = ka * kb;
kt = ka * tb + ta * kb;
tt = ta * tb;
}
D_FFT(tt, kk, N, -1); FFT(kt, N, -1);
for(int i = 0; i < N; i++){
tt = tt / N;
kt = kt / N;
kk = kk / N;
}
}
}PLE;
void main() {
scanf(&#34;%d%d%d&#34;, &n, &m, &p); n++; m++;
for(int i = 0; i < n; i++) scanf(&#34;%d&#34;, &a),a = a % p;
for(int i = 0; i < m; i++) scanf(&#34;%d&#34;, &b),b = b % p;
PLE.mult(a, b, n + m - 1);
for(int i = 0; i < n + m - 1; i++)
printf(&#34;%lld &#34;,(((((LL)round(kk.x)) % p) << 30) + ((((LL)round(kt.x)) % p) << 15) + ((LL)round(tt.x)) % p) % p);
}
}
blue\_stein 模板：
struct polynie {
CP getw(int m, int k, int op) {
return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m));
}
int wz[MAXN];
CP A[MAXN], B[MAXN], C[MAXN];
void FFT(CP *a, int N, int op) {
rop(i, 0, N) if(i < wz) swap(a, a[wz]);
for(int l = 2; l <= N; l <<= 1) {
int mid = l >> 1;
CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid));
for(int i = 0; i < N; i += l) {
w = CP(1, 0);
rop(j, 0, mid) {
x = a[i + j];
y = w * a[i + j + mid];
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void mult(CP *a, CP *b, CP *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
rop(i, 0, N) wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
FFT(a, N, 1); FFT(b, N, 1);
rop(i, 0, N) c = a * b;
FFT(c, N, -1);
rop(i, 0, N) c.x = c.x / N, c.y = c.y / N;
}
void blue_stein(CP *a, int M, int op) {
int M2 = M << 1;
memset(A, 0, sizeof(A));
memset(B, 0, sizeof(B));
memset(C, 0, sizeof(C));
rop(i, 0, M) A = a * getw(M2, 1ll * i * i % M2, op);
rop(i, 0, M2) B = getw(M2, 1ll * (i - M) * (i - M) % M2, -op);
mult(A, B, C, M2 + M - 1);
rop(i, 0, M) a = C[i + M] * getw(M2, 1ll * i * i % M2, op);
if(op == -1) rop(i, 0, M) a.x = a.x / M, a.y = a.y / M;
}
}PLE;
多项式全集之三 多项式求逆:
多项式求逆：
1、问题描述：
已知 F(x) ，且 F(x)G(x)\equiv 1 (mod~x^n) ，求 G(x)
2、推导过程：
设 B(x)\equiv F(x)^{-1}(mod~x^{\lceil \frac n 2 \rceil})
由于 G(x)\equiv F(x)^{-1} (mod~x^n)
所以 G(x)\equiv F(x)^{-1} (mod~x^{\lceil \frac n 2 \rceil})
那么
G(x)\equiv B(x)(mod~x^{\lceil \frac n 2 \rceil})
G(x)-B(x)\equiv 0(mod~x^{\lceil \frac n 2 \rceil})
两边平方，得：
由于 [G(x)-B(x)]^2 的第 k<n 项为
\sum_{i=0}^k[g_ix^i-b_ix^i][g_{k-i}x^{k-i}-b_{k-i}x^{k-i}]
i,k-i 一定有一项 <\frac n 2 ，所以
[G(x)-B(x)]^2\equiv 0 (mod~x^n)
G^2(x)+B^2(x)-2G(x)B(x)\equiv 0(mod~x^n)
两边同乘 A(x) ，得：
A(x)G^2(x)+A(x)B^2(x)-2A(x)G(x)B(x)\equiv 0 (mod~x^n)
G(x)+A(x)B^2(x)-2B(x)\equiv 0 (mod~x^n)
G(x)\equiv 2B(x)-A(x)B^2(x)(mod~x^n)
G(x)\equiv B(x)[2-A(x)B(x)](mod~x^n)
喜闻乐见的模板：
namespace INV{
typedef long long LL;
int n, a[300005], b[300005];
const int mod = 998244353;
const int g = 3189;
int qpow(int a, int b){
int ans = 1;
while(b){
if(b & 1) ans = 1ll * ans * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return ans;
}
struct p_l_e{
int wz[300005], i_c[300005];
void NTT(int *a, int N, int op){
for(int i = 0; i < N; i++)
if(i < wz) swap(a, a[wz]);
for(int le = 2; le <= N; le <<= 1){
int mid = le >> 1, wn = qpow(g, (mod - 1) / le);
if(op == -1) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += le){
LL w = 1; int x, y;
for(int j = 0; j < mid; j++){
x = a[i + j];
y = w * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
int init(int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
return N;
}
void INV(int *a, int *b, int deg){
if(deg == 1){b[0] = qpow(a[0], mod - 2); return;}
INV(a, b, (deg + 1) >> 1);
int N = init(deg + deg - 1);
for(int i = 0; i < deg; i++) i_c = a;
for(int i = deg; i < N; i++) i_c = 0;
NTT(b, N, 1);NTT(i_c, N, 1);
for(int i = 0; i < N; i++) b = 1ll * b * (2 - 1ll * b * i_c % mod + mod) % mod;
NTT(b, N, -1);
int t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) b = 1ll * b * t % mod;
for(int i = deg; i < N; i++) b = 0;
}
}PLE;
void main(){
scanf(&#34;%d&#34;, &n);
for(int i = 0; i < n; i++) scanf(&#34;%d&#34;, &a);
PLE.INV(a, b, n);
for(int i = 0; i < n; i++) printf(&#34;%d &#34;,b);
}
}
多项式全集之四 多项式ln:
1、做法：
设 G(x) =\ln F(x)
两边求导得
G&#39;(x)=\frac {F&#39;(x)} {F(x)}
积分回去即可。
2、应用：
e^{F(x)}=\sum_{k \ge 0}\frac{F^k(x)}{k!}=G(x)
F(x) = \ln G(x)
这个的组合意义是：无序组合。
设 F(x) ， f_i 表示一些东西，那么这些东西有序组合的方案数为
F^0(x) + F^1(x)+F^2(x)+\cdots=\frac 1 {1 - F(x)}
而无序组成的方案数为：
\frac{F^0(x)}{0!} + \frac {F^1(x)}{1!}+\frac{F^2(x)}{2!}+\cdots=e^{F(x)}
如果无序组合方案数好求，那么求 \ln 就能得到 F(x) 。
例题
喜闻乐见的代码：
多项式 \ln :
namespace PLE_ln{
struct polyme {
int li[SZ], wz[SZ];
void NTT(int *a, int N, int op) {
rop(i, 0, N) if(i < wz) swap(a, a[wz]);
for(int l = 2; l <= N; l <<= 1) {
int mid = l >> 1;
int x, y, w, wn = qpow(g, (mod - 1) / l);
if(op) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += l) {
w = 1;
for(int j = 0; j < mid; ++j) {
x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = 1ll * w * wn % mod;
}
}
}
}
void qd(int *a, int *b, int n) {
rop(i, 0, n) b = 1ll * a[i + 1] * (i + 1) % mod;
}
void jf(int *a, int *b, int n) {
rop(i, 1, n) b = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
}
void mult(int *a, int *b, int *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len ++;
rop(i, 0, N) wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
NTT(a, N, 0); NTT(b, N, 0);
rop(i, 0, N) c = 1ll * a * b % mod;
NTT(c, N, 1);
int t = qpow(N, mod - 2);
rop(i, 0, N) c = 1ll * c * t % mod;
}
void inv(int *a, int *b, int deg) {
if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;}
inv(a, b, (deg + 1) >> 1);
rop(i, 0, deg) li = a;
int N = 1, len = 0;
while(N < deg + deg - 1) N <<= 1, len ++;
rop(i, 0, N) wz = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
rop(i, deg, N) li = 0;
NTT(li, N, 0); NTT(b, N, 0);
rop(i, 0, N) b = 1ll * b * (2 - 1ll * li * b % mod + mod) % mod;
NTT(b, N, 1);
int t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) b = 1ll * b * t % mod;
rop(i, deg, N) b = 0;
}
}PLE;
int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n;
void main() {
scanf(&#34;%d&#34;, &n);
rop(i, 0, n) scanf(&#34;%d&#34;, &a);
PLE.qd(a, da, n);
PLE.inv(a, ia, n);
PLE.mult(ia, da, dla, n + n - 1);
PLE.jf(dla, la, n);
rop(i, 0, n) printf(&#34;%d &#34;, la);
}
}
与作者交流：https://ac.nowcoder.com/discuss/175409
更多算法和题解：https://ac.nowcoder.com/acm/contest/discuss |
本帖子中包含更多资源
您需要 登录 才可以下载或查看，没有账号？立即注册
×
|
窥视卡
雷达卡
发表于 2023-2-12 10:49
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
千斤顶
显身卡