[学习笔记] 牛顿迭代法
泰勒展开
假设我们有函数\(f(x)\)，\(x\)在\(x_0\)处有\(\infty\)阶导数。
我们知道\(f(x_0)\) 的值
我们希望构造一个多项式\(g\)，使它尽可能逼近函数\(f\)。
-
那么就使它在\(x_0\)处的1~\(\infty\)阶导数都与\(f\)相同，即：\(f^{(n)}=g^{(n)}\)，那么一直到最后就有\(f(x)=g(x)\)
-
我们还原这个多项式：
\(\begin{aligned}g(x)&=\xi_0+f(x_0)+\frac{f^{(1)}(x_0)}{1!}(x-x_0)^1+\frac{f^{(2)}(x_0)}{2!}(x-x_0)^2+...+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n\\&=\xi_0+f(x_0)+\sum_{i=1}^{\infty}\frac{f^{(i)}(x_0)}{i!}(x-x_0)^i\end{aligned}\)
这即是泰勒展开。
牛顿迭代法
假设我们要求：\(F(B(x))\equiv0(\mod x^{2^t})\)，\(B(x)\)为多项式
设 \(B_t(x)\equiv B(x)(\mod x^{2^t})\)
首先，我们可以很容易的算出\(B_0\)。
假设我们已经求出了\(B_{t-1}\)。
- 我们对\(F(B_{t+1}(x))\)在\(B_t\)处进行泰勒展开。
\(\begin{aligned}F(B_{t+1}(x))&=F(B_t(x))+\sum_{i=1}^n\frac{F^{(i)}(B_t(x))\times(B_{t+1}(x)-B_t(x))^i}{i!}\\&=F(B_t(x))+\frac{F'(B_t(x))\times(B_{t+1}(x)-B_t(x))}{1!}\end{aligned}\)
因为 \(F(B_t(x))\equiv0(\mod x^{2^t})\)，
所以 \((B_{t+1}(x)-B_t(x))^n=0(n\ge2)\)
- 把目标项单独移到左边：
\(B_{t+1}(x)=B_t(x)+\frac{F(B_{t+1}(x))-F(B_t(x))}{F'(B_t(x))}\)
- 又因为\(F(B_{t+1}(x))\)在模意义下等于0（定义），又有：
\(B_{t+1}(x)=B_t(x)-\frac{F(B_{t}(x))}{F'(B_t(x))}\)
这即使牛顿迭代法的式子。
多项式求逆
即求：\(A(x)\times B(x)\equiv1\) （\(A\)为原多项式系数）
有：
\(F(B_t(x))=A(x)\times B_t(x)-1\)
- 代入牛顿迭代法的式子：
\(\begin{aligned}B_{t+1}(x)&=B_t(x)-\frac{A(x)\times B_t(x)-1}{A(x)}\\&=B_t(x)-B_t(x)\times(A(x)\times B_{t}(x)-1)\\&=-A(x)\times B_t^2(x)+2B_t(x)\end{aligned}\)
（因为\(A\times B\equiv 1\)，所以除以\(A\)即是乘上\(B\)）
void Inv(int *f,int *g,int n){
int c[N];
f[0]=pow(g[0],mod-2);
for (int l=2;l<n<<1;l<<=1){
Get(l<<1,0);
for (int i=0;i<l;i++) c[i]=g[i];
for (int i=l;i<L;i++) c[i]=0;
fft(f,1),fft(c,1);
for (int i=0;i<L;i++) f[i]=(2-1ll*c[i]*f[i])%mod*f[i]%mod;
fft(f,-1);
for (int i=l;i<L;++i) f[i]=0;
}
}
多项式求 \(\ln\)
即求：\(\ln A(x)=B(x)\)
- 这次不需要用到牛顿迭代法，直接求导：
\(\begin{aligned}\ln A(x)&=B(x)\\\frac{A'(x)}{A(x)}&=B'(x)\end{aligned}\)
然后只要积分即可。
void Qiud(int *f,int *g,int n){
f[n-1]=0;
for (int i=0;i<n-1;i++) f[i]=1ll*g[i+1]*(i+1)%mod;
}
void Jif(int *f,int *g,int n){
f[0]=0;
for (int i=1;i<n;i++) f[i]=1ll*g[i-1]*inv[i]%mod;
}
void Ln(int *f,int *g,int n){
int c[N],d[N];
for (int i=0;i<n;i++) c[i]=d[i]=0;
Inv(c,g,n),Qiud(d,g,n);
Get(n<<1);
for (int i=n;i<L;i++) c[i]=d[i]=0;
fft(c,1),fft(d,1);
for (int i=0;i<L;i++) d[i]=1ll*c[i]*d[i]%mod;
fft(d,-1);
Jif(f,d,n);
}
\(exp\)
即求：\(e^{A(x)}=B(x)\)
- 我们先用\(\ln\) 把exp给去掉：
\(A(x)=\ln B(x)\)
有：
\(F(B_t(x))=\ln B(x)-A(x)\)
- 套牛顿迭代法：
\(\begin{aligned}B_{t+1}(x)&=B_t(x)-\frac{\ln B_t(x)-A(x)}{\frac{1}{B_t(x)}}\\&=B_t(x)\times(1-\ln B_t(x)+A(x))\end{aligned}\)
- 时间复杂度是 \(\varTheta(n\log^2n)\) 的
void Exp(int *f,int *g,int n){
int c[N],d[N];
f[0]=1;
for (int l=2;l<n<<1;l<<=1){
Ln(d,f,l);
Get(l<<1,0);
for (int i=0;i<l;i++) c[i]=g[i];
for (int i=l;i<L;i++) c[i]=d[i]=0;
fft(f,1),fft(c,1),fft(d,1);
for (int i=0;i<L;i++) f[i]=(1ll-d[i]+c[i])*f[i]%mod;
fft(f,-1);
}
}
以上就是全部内容了，刚发现没写过\(FFT\)的博客，待会儿补一下。
