文章目录
多项式的运算
多项式的加减法,数乘
这个大家应该都会吧
不推了。
直接上代码:
friend inline poly operator+(const poly&a,const poly&b){
poly ret(max(a.deg(),b.deg()));
for(ri i=0;i<=a.deg();++i)ret[i]=add(ret[i],a[i]);
for(ri i=0;i<=b.deg();++i)ret[i]=add(ret[i],b[i]);
return ret;
}
friend inline poly operator-(const poly&a,const poly&b){
poly ret(max(a.deg(),b.deg()));
for(ri i=0;i<=a.deg();++i)ret[i]=add(ret[i],a[i]);
for(ri i=0;i<=b.deg();++i)ret[i]=dec(ret[i],b[i]);
return ret;
}
friend inline poly operator*(const int&a,const poly&b){
poly ret(b.deg());
for(ri i=0;i<=b.deg();++i)ret[i]=mul(a,b[i]);
return ret;
}
多项式乘法
这个大家应该都会吧
还是推一推吧。
已知的:
A ( x ) = ∑ i = 0 n a i x i A(x)=\sum_{i=0}^na_ix^i A(x)=∑i=0naixi
B ( x ) = ∑ i = 0 m b i x i B(x)=\sum_{i=0}^mb_ix^i B(x)=∑i=0mbixi
要求的:
C ( x ) = A ( x ) B ( x ) = ∑ i = 0 n + m ( ∑ j = 0 m i n { i , n } a j ∗ b i − j ) x i C(x)=A(x)B(x)=\sum_{i=0}^{n+m}(\sum_{j=0}^{min\{i,n\}}a_j*b_{i-j})x^i C(x)=A(x)B(x)=∑i=0n+m(∑j=0min{
i,n}aj∗bi−j)xi
显然直接暴力做是 O ( n 2 ) O(n^2) O(n2)的,考虑如何优化。
那么我们使用 f f t fft fft或者 n t t ntt ntt来实现点值表示法和系数表示法之间的快速转化。
为了方便起见,我们将 A , B A,B A,B的最高次数统一成一个 2 2 2的幂(对于超过 n / m n/m n/m的项的系数看成0即可)
所谓的系数表示法就是我们平常用的那种。
而点值表示法,就是把这个多项式理解成一个函数,用这个函数上的若干个点的坐标来描述这个多项式:
f ( x ) = ( x 0 , y 0 ) , ( x 1 , y 1 ) , . . . , ( x n , y n ) = p 0 , p 1 , p 2 , . . . , p n f(x)=(x_0,y_0),(x_1,y_1),...,(x_n,y_n)=p_0,p_1,p_2,...,p_n f(x)=(x0,y0),(x1,y1),...,(xn,yn)=p0,p1,p2,...,pn
假设我们已经将 A , B A,B A,B两个函数转化成了点值表示,于是就可以马上求出 C C C的点值表示:
A ( x ) = ( x a , 0 , y a , 0 ) , ( x a , 1 , y a , 1 ) , . . . , ( x a , n − 1 , y a , n − 1 ) = p a , 0 , p a , 1 , . . . , p a , n − 1 A(x)=(x_{a,0},y_{a,0}),(x_{a,1},y_{a,1}),...,(x_{a,n-1},y_{a,n-1})=p_{a,0},p_{a,1},...,p_{a,n-1} A(x)=(xa,0,ya,0),(xa,1,ya,1),...,(xa,n−1,ya,n−1)=pa,0,pa,1,...,pa,n−1
B ( x ) = ( x b , 0 , y b , 0 ) , ( x b , 1 , y b , 1 ) , . . . , ( x b , n − 1 , y b , n − 1 ) = p b , 0 , p b , 1 , . . . , p b , n − 1 B(x)=(x_{b,0},y_{b,0}),(x_{b,1},y_{b,1}),...,(x_{b,n-1},y_{b,n-1})=p_{b,0},p_{b,1},...,p_{b,n-1} B(x)=(xb,0,yb,0),(xb,1,yb,1),...,(xb,n−1,yb,n−1)=pb,0,pb,1,...,pb,n−1
那么 C ( x ) = p a , 0 ∗ p b , 0 , p a , 1 ∗ p b , 1 , . . . , p a , n − 1 ∗ p b , n − 1 = p c , 0 , p c , 1 , . . . , p c , n − 1 C(x)=p_{a,0}*p_{b,0}, p_{a,1}*p_{b,1},...,p_{a,n-1}*p_{b,n-1}=p_{c,0},p_{c,1},...,p_{c,n-1} C(x)=pa,0∗pb,0,pa,1∗pb,1,...,pa,n−1∗pb,n−1=pc,0,pc,1,...,pc,n−1
然后再把 C ( x ) C(x) C(x)还原成系数表达式即可。
注意:我们需要保证 x a i , x b i x_{ai},x_{bi} xai,xbi互不相同)
现在就只用考虑如何实现点值表示和系数表示的互换了。
也就是如何用更少的计算次数来求出 n n n个不同的 x x x值对应的 y y y值。
考虑有两个具有特殊性质的东西:
单位根保证了 w n 0 , w n 1 , . . . , w n n − 1 w_n^0,w_n^1,...,w_n^{n-1} wn0,wn1,...,wnn−1是互不相同的并且有 w n i j = ( w n i ) j w_n^{ij}=(w_n^i)^j wnij=(wni)j
而原根在模数为质数 p p p的时候也有 g 0 , g 1 , g 2 , . . . , g n − 1 g^0,g^1,g^2,...,g^{n-1} g0,g1,g2,...,gn−1是互不相同的并且 g i j ≡ ( g i ) j m o d    p g^{ij}\equiv (g^i)^j \mod p gij≡(gi)jmodp
这满足了我们上面的性质,因此我们考虑将 w n 0 , w n 1 , . . . w n n − 1 w_n^0,w_n^1,...w_n^{n-1} wn0,wn1,...wnn−1作为 x 0 , x 1 , . . . x n − 1 x_0,x_1,...x_{n-1} x0,x1,...xn−1带入求点值。
然后要用到两个引理:
- 折半引理: w n k ∗ 2 = w n 2 k w_n^{k*2}=w_{\frac n2}^k wnk∗2=w2nk(n为偶数 )
- 消去引理: w n k = − w n k + n 2 w_n^{k}=-w_n^{k+\frac n2} wnk=−wnk+2n
这两个引理可以画个单位圆简单证明(建议各位神犇自己简要证明一下)
然后利用按照下标的奇偶性来进行分治处理:
f ( x ) = ∑ i = 0 n a i x i f(x)=\sum_{i=0}^na_ix^i f(x)=∑i=0naixi
=> f ( w n k ) = ∑ i = 0 n a i w n i k f(w_n^k)=\sum_{i=0}^na_i w_n^{ik} f(wnk)=∑i=0naiwnik
=> f ( w n k ) = ∑ i = 0 n 2 − 1 a 2 i w n 2 i k + w n k ∑ i = 0 n 2 − 1 a 2 i + 1 w n 2 i k f(w_n^k)=\sum_{i=0}^{\frac n2-1}a_{2i}w_n^{2ik}+w_n^k\sum_{i=0}^{\frac n2-1}a_{2i+1}w_n^{2ik} f(wnk)=∑i=02n−1a2iwn2ik+wnk∑i=02n−1a2i+1wn2ik
同时又有:
f ( w n k + n 2 ) = ∑ i = 0 n 2 − 1 a 2 i w n 2 i k − w n k ∑ i = 0 n 2 − 1 a 2 i + 1 w n 2 i k f(w_n^{k+\frac n2})=\sum_{i=0}^{\frac n2-1}a_{2i}w_n^{2ik}-w_n^k\sum_{i=0}^{\frac n2-1}a_{2i+1}w_n^{2ik} f(wnk+2n)=∑i=02n−1a2iwn2ik−wnk∑i=02n−1a2i+1wn2ik
这一步需要用到引理
所以我们只要算出两个重新分配了系数的多项式的值就可以了。
显然一直分下去只有 l o g log log层。
于是总时间复杂度 O ( n l o g n ) O(nlogn) O(nlogn)
注意到递归的效率很低,我们可以预处理最后一层的系数的位置然后用迭代的方式还原回去来优化常数
代码:
vector<int>pos;
int lim,tim;
inline void init(int up){
lim=1,tim=0;
while(lim<=up)lim<<=1,++tim;
pos.resize(lim),pos[0]=0;
for(ri i=0;i<lim;++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<(tim-1));
}
原理:
我们发现分治完之后相当于将 i i i和 i i i对应的二进制数在 l i m lim lim下反转之后对应的新二进制数 i ′ i' i′这两个位置的系数 a i , a i ′ a_i,a_i' ai,ai′交换了位置。
这就成功的实现了系数转点值。
下面来看点值转系数:
注:下面的图片有几点没写清楚:
- 最后的时候 e i , i = 0 , e i , j ∣ i ̸ = j e_{i,i}=0,e_{i,j|i\not=j} ei,i=0,ei,j∣i̸=j
- 那个 I n I_n In指的就是对角矩阵。
以上是 f f t fft fft的证明, n t t ntt ntt同理。
代码:
fft:
//复数定义
struct Complex{
double x,y;
friend