「学习笔记」多项式的蛇皮操作
写的时候注意各种数组的清空
前置知识
趋近
数学公式中,有类似于 ← \leftarrow ← 或者 → \rightarrow → 的东西,叫做趋近。
其中,前者叫做右趋近,后者叫做左趋近。
x → y x\rightarrow y x→y 表示 x x x 无限接近于 y y y 同时 x < y x<y x<y 。
y ← x y\leftarrow x y←x 表示 x x x 无限接近于 y y y 同时 x > y x>y x>y 。
注意,都是严格小于/大于 而没有取等。
那么, x → ∞ x\rightarrow \infty x→∞ 叫做 x x x 无限接近于 无限大 但同时 x < ∞ x<\infty x<∞ 。
同理, x → 0 x\rightarrow 0 x→0 叫做 x x x 无限接近于 0 0 0 但同时保证 x < 0 x<0 x<0 。
反之亦然。
自然常数
自然常数是一个叫做 e e e 的东西,那么他是什么呢?
它的其中一个定义为 e = lim x → ∞ ( x + 1 x ) x e=\lim_{x\rightarrow \infty}(x+\frac{1}{x})^x e=limx→∞(x+x1)x 。
这个不用怎么理解,记住就好了…
对数
当满足 a n = x a^n=x an=x 时,我们定义 log a x = n \log_a^x=n logax=n
在数学中,对数是对求幂的逆运算,正如除法是乘法的倒数,反之亦然。 这意味着一个数字的对数是必须产生另一个固定数字(基数)的指数。 在简单的情况下,乘数中的对数计数因子。更一般来说,乘幂允许将任何正实数提高到任何实际功率,总是产生正的结果,因此可以对于 b b b 不等于 1 1 1 的任何两个正实数 b b b 和 x x x 计算对数。
如果 a a a 的 x x x 次方等于 N N N ( a > 0 a>0 a>0 ,且 a a a 不等于 1 1 1 ) 即 a x = N a^x=N ax=N ,那么数 x x x 叫做以 a a a 为底 N N N 的对数 ( l o g a r i t h m ) (logarithm) (logarithm) ,记作 x = log a N x=\log_a^N x=logaN。其中, a a a 叫做对数的底数, N N N 叫做真数。
特别地,以无理数 e e e 为底记为 ln \ln ln ,称为自然对数。
逆元
这里指取模意义下的逆元。
被各种数论知识折磨过的大神们应该对这个很熟悉吧,这个就不多说了。
导函数
如果函数 f ( x ) f(x) f(x) 在 ( a , b ) (a,b) (a,b) 中每一点处都可导,则称 f ( x ) f(x) f(x) 在 ( a , b ) (a,b) (a,b) 上可导,则可建立 f ( x ) f(x) f(x) 的导函数,简称导数,记为 f ′ ( x ) f'(x) f′(x) 。
看似高深莫测,实际上就是求一个函数 f ( x ) f(x) f(x) 在 x = x 0 x=x_0 x=x0 时的 Δ \Delta Δ 。
那么,我们有 f ′ ( x ) = lim Δ x → 0 f ( x + Δ x ) − f ( x ) Δ x f'(x)=\lim_{\Delta x\rightarrow 0}\frac{f(x+\Delta x)-f(x)}{\Delta x} f′(x)=limΔx→0Δxf(x+Δx)−f(x) 。
牛顿迭代与泰勒公式
这是什么?怎么全都不知道啊…
首先,牛顿迭代是基于泰勒公式的,那么我们先从泰勒公式说起:
泰勒公式是将一个在 x = x 0 x=x_0 x=x0 处具有 n n n 阶导数的函数 f ( x ) f(x) f(x) 利用关于 ( x − x 0 ) (x-x_0) (x−x0) 的 n n n 次多项式来逼近函数的方法。
请注意,这是是逼近,而非等于。
具体的公式呈现:
f
(
x
)
=
f
(
x
0
)
0
!
+
f
′
(
x
0
)
1
!
(
x
−
x
0
)
+
f
′
′
(
x
0
)
2
!
(
x
−
x
0
)
2
+
…
+
f
(
n
)
(
x
0
)
n
!
(
x
−
x
0
)
n
+
R
n
(
x
)
f(x)=\frac{f(x_0)}{0!}+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\ldots +\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x)
f(x)=0!f(x0)+1!f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+…+n!f(n)(x0)(x−x0)n+Rn(x)
对这个式子进行一些说明:
f ( x ) f(x) f(x) :我们要求的目标函数。
f ( n ) ( x 0 ) f^{(n)}(x_0) f(n)(x0) :函数 f ( x ) f(x) f(x) 的 n n n 阶导数。
等号后的多项式称为函数 f ( x ) f(x) f(x) 在 x 0 x_0 x0 处的泰勒展开式,剩余的 R n ( x ) R_n(x) Rn(x) 是泰勒公式的余项,是 ( x − x 0 ) (x-x_0) (x−x0) 的 n n n 的高阶无穷小。
高阶无穷小:若 lim ( β / α ) = 0 \lim(β/α)=0 lim(β/α)=0 ,则称“ β β β 是比 α α α 较高阶的无穷小”。意思是在某一过程 ( x → x 0 x→x_0 x→x0 或 x → ∞ x→∞ x→∞ 这类过程) 中, β → 0 β→0 β→0 比 α → 0 α→0 α→0 快一些。
解释完泰勒公式,现在来说牛顿迭代…
多数方程不存在求根公式,因此求精确根非常困难,甚至不可能,从而寻找方程的近似根就显得特别重要。方法使用函数 f ( x ) f(x) f(x) 的泰勒级数的前面几项来寻找方程 f ( x ) = 0 f(x)=0 f(x)=0 的根。
用牛顿迭代法解非线性方程,是把非线性方程
f
(
x
)
=
0
f(x)=0
f(x)=0 线性化的一种近似方法。把
f
(
x
)
f(x)
f(x) 在点
x
0
x_0
x0 的某邻域内展开成泰勒级数:
f
(
x
)
=
f
(
x
0
)
0
!
+
f
′
(
x
0
)
1
!
(
x
−
x
0
)
+
f
′
′
(
x
0
)
2
!
(
x
−
x
0
)
2
+
…
+
f
(
n
)
(
x
0
)
n
!
(
x
−
x
0
)
n
+
R
n
(
x
)
=
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
+
f
′
′
(
x
0
)
2
!
(
x
−
x
0
)
2
+
…
+
f
(
n
)
(
x
0
)
n
!
(
x
−
x
0
)
n
+
R
n
(
x
)
f(x)=\frac{f(x_0)}{0!}+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\ldots +\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x) \\ =f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\ldots +\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x)
f(x)=0!f(x0)+1!f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+…+n!f(n)(x0)(x−x0)n+Rn(x)=f(x0)+f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+…+n!f(n)(x0)(x−x0)n+Rn(x)
取其线性部分(即泰勒展开的前两项),并令其等于0,即
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
=
0
f(x_0)+f'(x_0)(x-x_0)=0
f(x0)+f′(x0)(x−x0)=0
以此作为非线性方程
f
(
x
)
=
0
f(x)=0
f(x)=0
的近似方程,若
f
′
(
x
0
)
≠
0
f'(x_0)\neq 0
f′(x0)=0
则其解为
x
1
=
x
0
−
f
(
x
0
)
f
′
(
x
0
)
x_1=x_0-\frac{f(x_0)}{f'(x_0)}
x1=x0−f′(x0)f(x0)
这样,得到牛顿迭代法的一个迭代关系式:
x
n
+
1
=
x
n
−
f
(
x
n
)
f
′
(
x
n
)
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
xn+1=xn−f′(xn)f(xn)
不定积分与定积分
积分似乎是一个很高级的东西,但是不幸的是我们现在必须了解(或者说掌握?)其基本定义与基本运算。
不定积分:在微积分中,一个函数 f f f 的不定积分,或原函数,或反导数,是一个导数等于 f f f 的函数 F F F ,即 F ′ = f F ′ = f F′=f 。
什么意思?其实就是我们根据一个导数 f ′ f' f′ 去求一个函数 F F F 满足 F ′ = f ‘ F'=f‘ F′=f‘ 。
但为什么是 不定 ?
显然,导数与函数的常数项是没有关系的,也就是说,对于一个导数 f ′ f' f′ ,如果我们要去求其原函数 f f f ,是求不出其原函数的 常数项 的,即常数项不定,所以我们只能退而求次,求到 F F F (上文有解释),并称 F F F 是 f ′ f' f′ 的不定积分。
那么,我们就可以把这个写作
∫
f
′
=
F
(
x
)
+
C
\int f'=F(x)+C
∫f′=F(x)+C
或者
∫
f
′
(
x
)
d
x
=
F
(
x
)
+
C
\int f'(x)dx=F(x)+C
∫f′(x)dx=F(x)+C
其中, C C C 是任意常数。
那么,定积分是什么?
打一个比方,假如说 小明 在 2019.1.1 2019.1.1 2019.1.1 有 C C C 元钱,在 2020.1.1 2020.1.1 2020.1.1 他有了 C + Δ C+\Delta C+Δ 元钱,那么,我们可以根据作差,求到 Δ \Delta Δ ,但是,如果我们要求到 C C C ,似乎是不行的。
这个例子有什么用?
C C C 就是我们上文提及的 f f f 的常数项,是不可求的,而 Δ \Delta Δ 是可求得,这就是定积分,是能够确定的。
并且, Δ \Delta Δ 似乎是一个常数,也就是说:
定积分定义:设函数f(x) 在区间[a,b]上连续,将区间[a,b]分成n个子区间[x0,x1], (x1,x2], (x2,x3], …, (xn-1,xn],其中x0=a,xn=b。可知各区间的长度依次是: △ x 1 = x 1 − x 0 △x1=x1-x0 △x1=x1−x0 ,在每个子区间 ( x i − 1 , x i ] (xi-1,xi] (xi−1,xi] 中任取一点 ξ i ( 1 , 2 , . . . , n ) ξi(1,2,...,n) ξi(1,2,...,n) ,作和式 ∑ i = 1 n f ( ξ i ) Δ x i \sum_{i=1}^nf(\xi_i)\Delta x_i ∑i=1nf(ξi)Δxi 。该和式叫做积分和,设 λ = m a x { △ x 1 , △ x 2 , … , △ x n } λ=max\{△x1, △x2, …, △xn\} λ=max{△x1,△x2,…,△xn} (即λ是最大的区间长度),如果当 λ → 0 λ→0 λ→0 时,积分和的极限存在,则这个极限叫做函数 f ( x ) f(x) f(x) 在区间 [ a , b ] [a,b] [a,b] 的定积分,记为 ∫ a b f ( x ) d x \int_a^b f(x)dx ∫abf(x)dx ,并称函数 f ( x ) f(x) f(x) 在区间 [ a , b ] [a,b] [a,b] 上可积。
接下来从几何意义上说明定积分。
终于上图了
设图中的曲线函数为
f
(
x
)
f(x)
f(x) ,现在有一个定积分:
∫
P
0
P
1
f
(
x
)
d
x
\int_{P_0}^{P_1}f(x)dx
∫P0P1f(x)dx
其表示的含义就是函数
f
(
x
)
f(x)
f(x) 的图像在
P
0
<
x
⩽
P
1
P_0< x \leqslant P_1
P0<x⩽P1 所围出来的图形面积。
怎么求?把他们拆成很多很多个矩形,将这些矩形面积和累加,即可求得面积。
但是这样是会有误差的,但是我们发现,当 P 0 、 P 1 P_0、P_1 P0、P1 越来越大,或者是越来越接近时,我们的计算就越来越精准。
多项式乘法
这里,我们有十分优秀的 F F T FFT FFT 以及 N T T NTT NTT 算法可以解决这样的问题,具体的见下面博客,这里只贴板子了…
时间复杂度均为 O ( n log n ) \mathcal O(n\log n) O(nlogn) 。
板子: F F T FFT FFT 的模板, N T T NTT NTT 的模板可改
class fft_task{
private:
struct cplx{
double vr,vi;//实部和虚部
cplx(const double R=0,const double I=0):vr(R),vi(I){}//构造函数
//------------------overload----------------//
cplx operator + (const cplx a)const{return cplx(vr+a.vr,vi+a.vi);}//重载加法
cplx operator - (const cplx a)const{return cplx(vr-a.vr,vi-a.vi);}
cplx operator * (const cplx a)const{return cplx(vr*a.vr-vi*a.vi,vr*a.vi+a.vr*vi);}
cplx operator / (const double var)const{return cplx(vr/var,vi/var);}
};
int n,m;
cplx a[MAXN+5],b[MAXN+5];
void fft(cplx* f,const int len,const short opt=1){
//opt==-1 : FFT 的逆变换
if(!len)return;
cplx f0[len+5],f1[len+5];
for(int i=0;i<len;++i)
f0[i]=f[i<<1],f1[i]=f[i<<1|1];
fft(f0,len>>1,opt);
fft(f1,len>>1,opt);
cplx w=cplx(cos(Pi/len),opt*sin(Pi/len)),buf=cplx(1,0);
for(int i=0;i<len;++i,buf=buf*w){
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
}
}
public:
inline void launch(){
qread(n,m);
rep(i,0,n)scanf("%lf",&a[i].vr);
rep(i,0,m)scanf("%lf",&b[i].vr);
for(m+=n,n=1;n<=m;n<<=1);
fft(a,n>>1);
fft(b,n>>1);
rep(i,0,n-1)a[i]=a[i]*b[i];
fft(a,n>>1,-1);
rep(i,0,m)writc((int)((a[i].vr)/n+0.5),' ');
Endl;
}
}This;
多项式求逆元
给定一个
A
(
x
)
A(x)
A(x) ,求
A
−
1
(
x
)
A^{-1}(x)
A−1(x) 满足
A
(
x
)
A
−
1
(
x
)
≡
1
(
mod
x
n
)
A(x)A^{-1}(x)\equiv 1(\text{mod}\space x^n)
A(x)A−1(x)≡1(mod xn)
其中
(
mod
x
N
)
(\text{mod}\space x^N)
(mod xN) 即为舍去次数
⩾
n
\geqslant\space n
⩾ n 的项,只保留
0
0
0 到
n
−
1
n-1
n−1 次项。
而 A − 1 ( x ) A^{-1}(x) A−1(x) 就是 A ( x ) A(x) A(x) 在 mod x n \text{mod}\space x^n mod xn 下的逆元。
首先,我们设
B
(
x
)
=
A
−
1
(
x
)
B(x)=A^{-1}(x)
B(x)=A−1(x) ,那么就有
A
(
x
)
B
(
x
)
≡
1
(
mod
x
n
)
A(x)B(x)\equiv 1(\text{mod}\space x^n)
A(x)B(x)≡1(mod xn)
假设我们已经知道
A
(
x
)
A(x)
A(x) 在
mod
x
⌈
n
2
⌉
\text{mod}\space x^{\left \lceil \frac{n}{2} \right \rceil}
mod x⌈2n⌉ 情况下的逆元
B
0
(
x
)
B_0(x)
B0(x) ,那么就有
A
(
x
)
B
0
(
x
)
≡
1
(
mod
x
⌈
n
2
⌉
)
A(x)B_0(x)\equiv 1 (\text{mod}\space x^{\left \lceil \frac{n}{2} \right \rceil})
A(x)B0(x)≡1(mod x⌈2n⌉)
而第一个式子我们也可以写为
A
(
x
)
B
(
x
)
≡
1
(
mod
x
⌈
n
2
⌉
)
A(x)B(x)\equiv 1(\text{mod}\space x^{\left \lceil \frac{n}{2} \right \rceil})
A(x)B(x)≡1(mod x⌈2n⌉)
将这两个式子相减可得
A
(
x
)
(
B
(
x
)
−
B
0
(
x
)
)
≡
0
(
mod
x
⌈
n
2
⌉
)
A(x)(B(x)-B_0(x))\equiv 0(\text{mod}\space x^{\left \lceil \frac{n}{2} \right \rceil})
A(x)(B(x)−B0(x))≡0(mod x⌈2n⌉)
即
B
(
x
)
−
B
0
(
x
)
≡
0
(
mod
x
⌈
n
2
⌉
)
B(x)-B_0(x)\equiv 0(\text{mod}\space x^{\left \lceil \frac{n}{2} \right \rceil})
B(x)−B0(x)≡0(mod x⌈2n⌉)
平方之后
B
2
(
x
)
+
B
0
2
(
x
)
−
2
B
(
x
)
B
0
(
x
)
≡
0
(
mod
x
n
)
B^2(x)+B_0^2(x)-2B(x)B_0(x)\equiv 0(\text{mod}\space x^n)
B2(x)+B02(x)−2B(x)B0(x)≡0(mod xn)
由于一个多项式平方之后,次数
<
n
<n
<n 的项至少是由原先一个次数
<
⌈
n
2
⌉
<\lceil \frac{n}{2} \rceil
<⌈2n⌉ 的项乘上其他项得到的,所以这个结果的
0
0
0 到
n
−
1
n-1
n−1 次系数仍然是
0
0
0 ,可以变成
(
m
o
d
x
n
)
\pmod{x^n}
(modxn)
同乘
A
(
x
)
A(x)
A(x) 得
B
(
x
)
−
2
B
0
(
x
)
+
A
(
x
)
B
0
2
(
x
)
≡
0
(
m
o
d
x
n
)
B(x)-2B_0(x)+A(x)B_0^2(x)\equiv0\pmod{x^n}
B(x)−2B0(x)+A(x)B02(x)≡0(modxn)
即
B
(
x
)
≡
B
0
(
x
)
(
2
−
A
(
x
)
B
0
(
x
)
)
(
m
o
d
x
n
)
B(x)\equiv B_0(x)(2-A(x)B_0(x))\pmod{x^n}
B(x)≡B0(x)(2−A(x)B0(x))(modxn)
根据定义,我们可以直接取等,得
B
(
x
)
=
B
0
(
x
)
(
2
−
A
(
x
)
B
0
(
x
)
)
mod
x
n
B(x)= B_0(x)(2-A(x)B_0(x))\space \text{mod}\space {x^n}
B(x)=B0(x)(2−A(x)B0(x)) mod xn
时间复杂度为
T
(
n
)
=
T
(
n
2
)
+
O
(
n
log
n
)
=
O
(
n
log
n
)
T(n)=T(\frac{n}{2})+\mathcal O(n\log n)=\mathcal O(n\log n)
T(n)=T(2n)+O(nlogn)=O(nlogn)
#include<cstdio>
#include<algorithm>
using namespace std;
#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
inline int qkpow(int a,int n,const int mod){
int ret=1;
for(;n>0;n>>=1){
if(n&1)ret=1ll*ret*a%mod;
a=1ll*a*a%mod;
}
return ret;
}
inline int qkpow(int a,int n){
int ret=1;
for(;n>0;n>>=1){
if(n&1)ret*=a;
a*=a;
}
return ret;
}
const int MOD=998244353,g=3,gi=332748118;
const int MAXN=3e5;
int revi[MAXN+5];
inline void ntt(int* A,const int n,const short opt=1){
rep(i,0,n-1)revi[i]=(revi[i>>1]>>1)|((i&1)?n>>1:0);
rep(i,0,n-1)if(i<revi[i])swap(A[i],A[revi[i]]);
for(int s=2,len,gn,Pow,tmp;s<=n;s<<=1){
len=s>>1,gn=qkpow(opt==1?g:gi,(MOD-1)/s,MOD);
for(int i=0;i<n;i+=s){Pow=1;
for(int j=i;j<len+i;++j,Pow=1ll*Pow*gn%MOD){
tmp=1ll*Pow*A[len+j]%MOD;
A[len+j]=(0ll+A[j]-tmp+MOD)%MOD;
A[j]=(0ll+A[j]+tmp)%MOD;
}
}
}
if(opt==-1){
int Inv=qkpow(n,MOD-2,MOD);
rep(i,0,n-1)A[i]=1ll*A[i]*Inv%MOD;
}
}
int tmp[MAXN+5],p;
void polyinv(int* A,int* B,const int n){
if(n==1)return (void)(B[0]=qkpow(A[0],MOD-2,MOD));//到了最下面, 直接取其逆元
polyinv(A,B,(n+1)/2);//n 的上取整
for(p=1;p<=(n<<1);p<<=1);//NTT 标准步骤
rep(i,0,n-1)tmp[i]=A[i];//为了避免 A 数组被更改,
rep(i,n,p)tmp[i]=0;
ntt(tmp,p),ntt(B,p);
rep(i,0,p-1)B[i]=((2-1ll*B[i]*tmp[i]%MOD)%MOD+MOD)*B[i]%MOD;
ntt(B,p,-1);
rep(i,n,p)B[i]=0;
}
int n,A[MAXN+5],B[MAXN+5];
signed main(){
#ifdef FILEOI
freopen("1.in","r",stdin);
freopen("file.out","w",stdout);
#endif
qread(n);
rep(i,0,n-1)qread(A[i]);
polyinv(A,B,n);
rep(i,0,n-1)writc(B[i],' ');
return 0;
}
多项式除法/取模
给定
n
−
1
n-1
n−1 次多项式
A
(
x
)
A(x)
A(x) 与
m
−
1
m-1
m−1 次多项式
B
(
x
)
B(x)
B(x) ,求
D
(
x
)
、
R
(
x
)
D(x)、R(x)
D(x)、R(x) 满足
A
(
x
)
=
D
(
x
)
B
(
x
)
+
R
(
x
)
A(x)=D(x)B(x)+R(x)
A(x)=D(x)B(x)+R(x)
其中
D
(
x
)
D(x)
D(x) 最高次项最多为
n
−
m
n-m
n−m ,
R
(
x
)
R(x)
R(x) 的次数
<
m
−
1
<m-1
<m−1 。
或者满足
A
(
x
)
≡
R
(
x
)
(
m
o
d
B
(
x
)
)
A(x)\equiv R(x)\pmod {B(x)}
A(x)≡R(x)(modB(x))
其实,从某种角度上来说,这就是多项式除法或者是取模,取决于你想求什么。
因为有 R ( x ) R(x) R(x) ,我们考虑将其去掉,减小问题的复杂度。
定义反转数组
A
R
(
x
)
=
x
n
−
1
A
(
1
x
)
=
∑
i
=
0
n
−
1
a
n
−
a
−
i
x
i
A^R(x)=x^{n-1}A(\frac{1}{x})=\sum_{i=0}^{n-1}a_{n-a-i}x^i
AR(x)=xn−1A(x1)=i=0∑n−1an−a−ixi
将
1
x
\frac{1}{x}
x1 代入第一个式子,并同时乘以
x
n
−
1
x^{n-1}
xn−1 ,可得
x
n
−
1
A
(
1
x
)
=
x
n
−
m
D
(
1
x
)
x
m
−
1
B
(
1
x
)
+
x
n
−
m
+
1
∗
x
m
−
2
R
(
1
x
)
x^{n-1}A(\frac{1}{x})=x^{n-m}D(\frac{1}{x})x^{m-1}B(\frac{1}{x})+x^{n-m+1}*x^{m-2}R(\frac{1}{x})
xn−1A(x1)=xn−mD(x1)xm−1B(x1)+xn−m+1∗xm−2R(x1)
根据反转数组的定义,可得
A
R
(
x
)
=
D
R
(
x
)
B
R
(
x
)
+
x
n
−
m
+
1
R
R
(
x
)
A^R(x)=D^R(x)B^R(x)+x^{n-m+1}R^R(x)
AR(x)=DR(x)BR(x)+xn−m+1RR(x)
由于
D
R
(
x
)
D^R(x)
DR(x) 是
n
−
m
n-m
n−m 次的,所以其在
(
m
o
d
x
n
−
m
+
1
)
\pmod{x^{n-m+1}}
(modxn−m+1) 下是没有影响的,那么,我们可以得到
A
R
(
x
)
≡
D
R
(
x
)
B
R
(
x
)
(
m
o
d
x
n
−
m
+
1
)
A^R(x)\equiv D^R(x)B^R(x)\pmod{x^{n-m+1}}
AR(x)≡DR(x)BR(x)(modxn−m+1)
观察上式,知道
A
R
(
x
)
、
B
R
(
x
)
A^R(x)、B^R(x)
AR(x)、BR(x) ,那么我们可以求得
D
R
(
x
)
D^R(x)
DR(x) ,再通过反转得到
D
(
x
)
D(x)
D(x) ,即
A
R
(
x
)
B
R
−
1
(
x
)
≡
D
R
(
x
)
(
m
o
d
x
n
−
m
+
1
)
A^R(x){B^R}^{-1}(x)\equiv D^R(x)\pmod {x^{n-m+1}}
AR(x)BR−1(x)≡DR(x)(modxn−m+1)
有一次逆元,两次乘法,时间复杂度是同样的
O
(
n
log
n
)
\mathcal O(n\log n)
O(nlogn) 。
而 R ( x ) R(x) R(x) 应该很好求了吧,这里不再赘述。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
inline int qkpow(int a,int n,const int mod){
int ret=1;
for(;n>0;n>>=1){
if(n&1)ret=1ll*ret*a%mod;
a=1ll*a*a%mod;
}
return ret;
}
inline int qkpow(int a,int n){
int ret=1;
for(;n>0;n>>=1){
if(n&1)ret*=a;
a*=a;
}
return ret;
}
const int MOD=998244353,g=3,gi=332748118;
const int MAXN=3e5;
int revi[MAXN+5];
inline void ntt(int* A,const int n,const short opt=1){
rep(i,0,n-1)revi[i]=(revi[i>>1]>>1)|((i&1)?n>>1:0);
rep(i,0,n-1)if(i<revi[i])swap(A[i],A[revi[i]]);
for(int s=2,len,gn,Pow,tmp;s<=n;s<<=1){
len=s>>1,gn=qkpow(opt==1?g:gi,(MOD-1)/s,MOD);
for(int i=0;i<n;i+=s){Pow=1;
for(int j=i;j<len+i;++j,Pow=1ll*Pow*gn%MOD){
tmp=1ll*Pow*A[len+j]%MOD;
A[len+j]=(0ll+A[j]-tmp+MOD)%MOD;
A[j]=(0ll+A[j]+tmp)%MOD;
}
}
}
if(opt==-1){
int Inv=qkpow(n,MOD-2,MOD);
rep(i,0,n-1)A[i]=1ll*A[i]*Inv%MOD;
}
}
int tmp[MAXN+5];
void polyinv(int* A,int* B,const int n){
if(n==1)return (void)(B[0]=qkpow(A[0],MOD-2,MOD));//到了最下面, 直接取其逆元
polyinv(A,B,(n+1)/2);//n 的上取整
int p;
for(p=1;p<=(n<<1);p<<=1);//NTT 标准步骤
rep(i,0,n-1)tmp[i]=A[i];//为了避免 A 数组被更改
rep(i,n,p)tmp[i]=0;
ntt(tmp,p),ntt(B,p);
rep(i,0,p-1)B[i]=((2-1ll*B[i]*tmp[i]%MOD)%MOD+MOD)*B[i]%MOD;
ntt(B,p,-1);
rep(i,n,p)B[i]=0;
}
int Ar[MAXN+5],Br[MAXN+5],Br_[MAXN+5];
inline void polydiv(int* A,int* B,int* D,int* R,int n,int m){
memset(Ar,0,sizeof Ar);memset(Br,0,sizeof Br);memset(Br_,0,sizeof Br_);
rep(i,0,n-1)Ar[i]=A[n-i-1];
rep(i,0,m-1)Br[i]=B[m-i-1];
polyinv(Br,Br_,n-m+1);
int p;
for(p=1;p<=n*2-m+1;p<<=1);
/*
因为 Br(x) 的逆元 Br_(x) 是在 mod x^{n-m+1} 的情况下进行的
因而其长度为 n-m+1
再加上 A(x) 的长度 n
总长即为 n+n-m+1=n*2-m+1
*/
ntt(Ar,p),ntt(Br_,p);
rep(i,0,p-1)D[i]=1ll*Ar[i]*Br_[i]%MOD;
ntt(D,p,-1);
rep(i,n-m+1,p)D[i]=0;//因为 D(x) 的最高次为 n-m , 因而后面的都可以全部清零
for(int i=0,j=n-m;i<j;++i,--j)swap(D[i],D[j]);//求得 Dr(x) , 还需反转
for(p=1;p<=n;p<<=1);
/*
n+1=n-m+m
n-m : 多项式 D
m : 多项式 B
*/
rep(i,0,m-1)Br[i]=B[i];rep(i,m,p)Br[i]=0;
rep(i,0,n-m)Br_[i]=D[i];rep(i,n-m+1,p)Br_[i]=0;
ntt(Br,p),ntt(Br_,p);
rep(i,0,p-1)R[i]=1ll*Br[i]*Br_[i]%MOD;
ntt(R,p,-1);
rep(i,0,n-1)R[i]=(0ll+A[i]-R[i]+MOD)%MOD;
}
int n,m,F[MAXN+5],G[MAXN+5],Q[MAXN+5],R[MAXN+5];
signed main(){
#ifdef FILEOI
freopen("1.in","r",stdin);
freopen("file.out","w",stdout);
#endif
qread(n,m);
rep(i,0,n)qread(F[i]);
rep(i,0,m)qread(G[i]);
polydiv(F,G,Q,R,n+1,m+1);
rep(i,0,n-m)printf("%d ",Q[i]);
Endl;
rep(i,0,m-1)printf("%d ",R[i]);
Endl;
return 0;
}
多项式牛顿迭代法
这个部分没有代码,只是一种思想。
有一个关于多项式 f ( x ) f(x) f(x) 的方程 d ( f ( x ) ) = 0 d(f(x))=0 d(f(x))=0
对 g ( f ( x ) ) g(f(x)) g(f(x)) 的说明(如果理解其含义可直接跳过):
这里的 f ( x ) f(x) f(x) 并非是一个值,而是多项式,比如有
g ( x ) = 4 x 2 + 1 g(x)=4x^2+1 g(x)=4x2+1
并且我们令
f ( x ) = x + 4 f(x)=x+4 f(x)=x+4
那么,就有
g ( f ( x ) ) = 4 ( x + 4 ) 2 + 1 g(f(x))=4(x+4)^2+1 g(f(x))=4(x+4)2+1
现在开始说明牛顿迭代法。
假设我们已经知道了
f
(
x
)
f(x)
f(x) 的前
n
n
n 项的多项式
f
0
(
x
)
f_0(x)
f0(x) ,即
f
(
x
)
≡
f
0
(
x
)
(
m
o
d
x
n
)
g
(
f
0
(
x
)
)
≡
0
(
m
o
d
x
n
)
f(x)\equiv f_0(x)\pmod{x^n} \\ g(f_0(x))\equiv 0\pmod{x^n}
f(x)≡f0(x)(modxn)g(f0(x))≡0(modxn)
然后,我们对
g
(
f
(
x
)
)
g(f(x))
g(f(x)) 在
f
0
(
x
)
f_0(x)
f0(x) 上进行泰勒展开
g
(
f
(
x
)
)
=
g
(
f
0
(
x
)
)
+
g
′
(
f
0
(
x
)
)
1
!
(
f
(
x
)
−
f
0
(
x
)
)
1
+
g
′
′
f
0
(
x
)
2
!
(
f
(
x
)
−
f
0
(
x
)
)
2
+
…
…
g(f(x))=g(f_0(x))+\frac{g'(f_0(x))}{1!}(f(x)-f_0(x))^1+\frac{g''f_0(x)}{2!}(f(x)-f_0(x))^2+\ldots \ldots
g(f(x))=g(f0(x))+1!g′(f0(x))(f(x)−f0(x))1+2!g′′f0(x)(f(x)−f0(x))2+……
根据定义,有
f
(
x
)
−
f
0
(
x
)
f(x)-f_0(x)
f(x)−f0(x) 的前
n
n
n 项系数为
0
0
0 ,那么就有
g
(
f
(
x
)
)
≡
g
(
f
0
(
x
)
)
+
g
′
(
f
0
(
x
)
)
(
f
(
x
)
−
f
0
(
x
)
)
≡
0
(
m
o
d
x
2
n
)
g(f(x))\equiv g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\equiv 0\pmod{x^{2n}}
g(f(x))≡g(f0(x))+g′(f0(x))(f(x)−f0(x))≡0(modx2n)
即
f
(
x
)
≡
f
0
(
x
)
−
g
(
f
0
(
x
)
)
g
′
(
f
0
(
x
)
)
(
m
o
d
x
2
n
)
f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^{2n}}
f(x)≡f0(x)−g′(f0(x))g(f0(x))(modx2n)
用这种方法也可以推多项式求逆。
多项式开根
给定多项式
A
(
x
)
A(x)
A(x) ,求
B
(
x
)
B(x)
B(x) 使得
B
2
(
x
)
−
A
(
x
)
≡
0
(
m
o
d
x
n
)
B^2(x)-A(x)\equiv0\pmod{x^n}
B2(x)−A(x)≡0(modxn)
设
B
(
x
)
≡
B
0
(
x
)
(
m
o
d
x
n
)
B(x)\equiv B_0(x)\pmod{x^n}
B(x)≡B0(x)(modxn) ,换句话说,设
B
(
x
)
B(x)
B(x) 与
B
0
(
x
)
B_0(x)
B0(x) 的前
n
n
n 项相同。
再令 g ( x ) = B 0 2 ( x ) − A ( x ) g(x)=B_0^2(x)-A(x) g(x)=B02(x)−A(x) ,那么其导函数 g ′ ( x ) = 2 B 0 ( x ) g'(x)=2B_0(x) g′(x)=2B0(x) 。
这与牛顿迭代的定义相似,那么直接带入牛顿迭代,可得
B
(
x
)
≡
B
0
(
x
)
−
B
0
2
(
x
)
−
A
(
x
)
2
B
0
(
x
)
≡
1
2
(
B
0
(
x
)
+
A
(
x
)
B
0
(
x
)
)
(
m
o
d
x
2
n
)
\begin{aligned} B(x)&\equiv B_0(x)-\frac{B_0^2(x)-A(x)}{2B_0(x)} \\ &\equiv\frac{1}{2}\left( B_0(x)+\frac{A(x)}{B_0(x)}\right)\pmod{x^{2n}} \end{aligned}
B(x)≡B0(x)−2B0(x)B02(x)−A(x)≡21(B0(x)+B0(x)A(x))(modx2n)
这样,我们也可以考虑用类似于倍增的思想。
时间复杂度 O ( n log n ) \mathcal O(n\log n) O(nlogn) 。
说明:牛顿迭代不是 逼近 吗?怎么可以拿来求根,万一不精准呢?
我们想一个问题,在我们普通的一些求根之中,不可能也算到精准吧?
比如 2 = 1.414... \sqrt2=1.414... 2=1.414... 而我们一般都取 1.414 1.414 1.414 等等。
多项式也是如此