多项式牛顿迭代及其运用
参考资料
牛顿迭代法在多项式运算的应用-by Miskcoo
如何通俗易懂地讲解牛顿迭代法？
前言-牛顿迭代
在食用本文之前，建议先学习这篇博客：多项式常用操作归纳
同样的，本文还是建议从前往后进行学习~~~
实数意义下的
首先是看了马老师的博客，然后就了解了求不规则函数根的方法。
下面是博主自己的概括和理解：
大概就是随便在x轴上找一个点，然后向上作x轴的垂线，交函数于一点y，然后再作(x,y)处的切线，交x轴于(x',0)。又从(x',0)这个点开始不断地重复。
最终我们找到的交x轴的那个点，有极大的概率是方程的根（函数的零点）。
现在我们来看一下，在已知\((x_0,f(x_0))\)的情况下，如何求出x'的值：
设原函数为\(f(x)\)，然后在\((x_0,f(x_0))\)的切线方程为：\(y=f'(x_0)(x-x_0)+f(x_0)\)，然后就可以简单的令y=0，就可以得到：\(0=f'(x_0)(x'-x_0)+f(x_0)\)->\(x'=x_0-\frac{f(x_0)}{f'(x_0)}\)，然后就可以通过这个式子进行不断地迭代了。
多项式意义下的
这个...其实还没有发现和上面的那个有什么关系...可能是我研究的还不够深入吧...先道个歉...
重点是记住结论就好了
首先看一个式子:
\[G(F(x))\equiv 0(mod\ x^n)
\]
已知的是G(x)，要求的是F(x)。
咋搞呢？
首先是多项式问题的常见套路：设\(F_0(x)\)是在\(mod\ x^{\lfloor\frac{n}{2}\rfloor}\)意义下的解，而且已经求出来了。
也就是有：
\[G(F_0(x))\equiv 0(mod\ x^{\lfloor\frac{n}{ 2}\rfloor})
\]
由泰勒展开可得：
\[G(F(x))=\frac{G(F_0(x))}{0!}+\frac{G'(F_0(x))}{1!}(F(x)-F_0(x))+\frac{G''(F_0(x))}{2!}(F(x)-F_0(x))^2+...
\]
然后我们进一步可以发现，其实从式子中的第三项开始就可以省略了。
为什么呢？？
因为在高次的情况下满足的式子，低次也是满足的：
\[G(F(x))\equiv 0(mod\ x^{\lfloor\frac{n}{2}\rfloor})
\]
进一步也就是说\(F(x)\)和\(F_0(x)\)的后n/2位是相等的。
然后你就会发现\((F(x)-F_0(x))^2\)的最低次项的次数都是(n/2)*(n/2)=n的，在模\(x^n\)的意义下，就变成0了！！！
所以说只有前面两项是有意义的。
这样子之后就变得炒鸡简单啦：
\[G(F(x))=\frac{G(F_0(x))}{0!}+\frac{G'(F_0(x))}{1!}(F(x)-F_0(x))
\]
然后我们进行进一步的变形，就能够得出F(x)的表达式：
\[F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}
\]
特别需要注意的是，这里的F(x)最好是看成是一个变元（通俗地讲，就看成是'x'），在这个意义下在进行求导的运算。
如何具体的使用这个蕴藏着丰富力量的公式呢？
每次都将题目中的同余式化成右边等于0 的形式，然后再令左边等于G(F(x))，再带入上述的结论式中，就能够得到你想要的结论了~~~
下面给出的常见操作大都是使用牛顿迭代得到的，其实也有常规的推导方式，只不过使用牛顿迭代更加简单易懂罢了。
多项式对数函数
这个是用不到牛顿迭代的...这里就先说了。
通常题目都是在模x^n的意义下进行求解的，后面就不再重复了。
题目就是：
\[B(x)\equiv ln(A(x))(mod\ x^n)
\]
其中A(x)是已知的，B(x)是我们要求的。
对两边同时求导可得：
\[B'(x)\equiv \frac{A'(x)}{A(x)}(mod\ x^n)
\]
然后我们发现后面的那个式子是能够使用之前的方法解决的（求逆+求导）。求导的话，不懂的可以出门百度一下，比较简单。
最后得到B'(x)之后，再积分回来，就能够得到最终的答案了。（积分和求导都是有比较简单的公式的，且是O(n)的）
板题：【模板】多项式对数函数
#include<cstdio>
#include<cstring>
#include<algorithm>
#define MAXN 2000000
#define MO 998244353
#define G 3
using namespace std;
int a[MAXN+5],a2[MAXN+5];
int n;
int PowMod(int x,int y)
{
int ret=1;
while(y)
{
if(y&1)
ret=1LL*ret*x%MO;
x=1LL*x*x%MO;
y>>=1;
}
return ret;
}
void NTT(int P[],int len,int oper)
{
for(int i=1,j=0;i<len-1;i++)
{
for(int s=len;j^=s>>=1,~j&s;);
if(i<j) swap(P[i],P[j]);
}
int unit,unit_p0;
for(int d=0;(1<<d)<len;d++)
{
int m=(1<<d),m2=m*2;
unit_p0=PowMod(G,(MO-1)/m2);
if(oper==-1)
unit_p0=PowMod(unit_p0,MO-2);
for(int i=0;i<len;i+=m2)
{
unit=1;
for(int j=0;j<m;j++)
{
int &P1=P[i+j+m],&P2=P[i+j];
int t=1LL*unit*P1%MO;
P1=((1LL*P2-1LL*t)%MO+MO)%MO;
P2=(1LL*P2+1LL*t)%MO;
unit=1LL*unit*unit_p0%MO;
}
}
}
if(oper==-1)
{
int inv=PowMod(len,MO-2);
for(int i=0;i<len;i++)
P[i]=1LL*P[i]*inv%MO;
}
}
void Mul(int ret[],int _x[],int l1,int _y[],int l2)
{
static int X[MAXN+5],Y[MAXN+5],RET[MAXN+5];
int len=1;
while(len<l1+l2) len<<=1;
copy(_x,_x+l1,X);
copy(_y,_y+l2,Y);
fill(X+l1,X+len,0);
fill(Y+l2,Y+len,0);
NTT(X,len,1);NTT(Y,len,1);
for(int i=0;i<len;i++)
RET[i]=1LL*X[i]*Y[i]%MO;
NTT(RET,len,-1);
for(int i=0;i<l1+l2-1;i++)
ret[i]=RET[i];
}
void Polynomial_Inverse(int deg,int A[],int B[])
{
static int tmpA[MAXN+5],tmpB[MAXN+5];
if(deg==1)
{
B[0]=PowMod(A[0],MO-2);
return;
}
Polynomial_Inverse((deg+1)/2,A,B);
int p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);
copy(B,B+deg,tmpB);
fill(tmpA+deg,tmpA+p,0);
fill(tmpB+deg,tmpB+p,0);
NTT(tmpA,p,1);NTT(tmpB,p,1);
for(int i=0;i<p;i++)
tmpB[i]=(1LL*tmpB[i]*(2LL-1LL*tmpA[i]*tmpB[i]%MO)%MO+MO)%MO;
NTT(tmpB,p,-1);
copy(tmpB,tmpB+deg,B);
}
int main()
{
scanf("%d",&n);
for(int i=0;i<n;i++)
scanf("%d",&a[i]);
Polynomial_Inverse(n,a,a2);
for(int i=0;i<n-1;i++)
a[i]=1LL*a[i+1]*(1LL*i+1LL)%MO;
a[n-1]=0;
Mul(a,a,n-1,a2,n);
for(int i=n-1;i>=1;i--)
{
int inv=PowMod(i,MO-2);
a[i]=1LL*a[i-1]*inv%MO;
}
a[0]=0;
for(int i=0;i<n;i++)
printf("%d ",a[i]);
return 0;
}
/*
2
1 1
*/
多项式指数函数
从这里开始就需要用到牛顿迭代了。
我们的问题是：
\[e^{A(x)}\equiv B(x)(mod\ x^n)
\]
其中A(x)是已知的，B(x)是我们要求的。
将问题转化一下，就变成了：
\[ln(B(x))\equiv A(x)(mod\ x^n)
\]
然后将A(x)移到左边来，就有了：
\[ln(B(x))-A(x)\equiv 0(mod\ x^n)
\]
再令\(G(F(x))=ln(B(x))-A(x)\)
最后来一波牛顿迭代，就有了：
\[B(x)=B_0(x)-\frac{ln(B_0(x))-A(x)}{\frac{1}{B_0(x)}}
\]
再变一下，就有了：
\[B(x)=B_0(x)(1-ln(B_0(x))+A(x))$$其中$B_0(x)$是我们在模$x^{\lfloor\frac{n}{2}\rfloor}$意义下已经求得的答案。求ln(B0(x))的话，可以用之前的**多项式对数函数**搞就可以了。
这里需要特别注意的是，多项式求指数函数是有前提的：A(x)的常数项必须为0，否则在底层的时候B[0]无法赋初值
这样子就有B[0]=1(因为ln(1)=0)。
板题：[【模板】多项式指数函数](https://www.luogu.org/problemnew/show/P4726)
```cpp
#include<cstdio>
#include<cstring>
#include<algorithm>
#define MAXN 5000000
#define MO 998244353
#define G 3
using namespace std;
typedef long long LL;
LL a[MAXN+5],b[MAXN+5];
LL inv[MAXN+5],n;
void Init()
{
inv[1]=1;
for(LL i=2;i<=MAXN;i++)
inv[i]=(1LL*MO-1LL*MO/i)*inv[MO%i]%MO;
}
void Diff(LL A[],LL len)
{
for(LL i=0;i<len-1;i++)
A[i]=1LL*A[i+1]*(1LL*i+1LL)%MO;
A[len-1]=0;
}
void Calc(LL A[],LL len)
{
for(LL i=len-1;i>=1;i--)
A[i]=1LL*A[i-1]*inv[i]%MO;
A[0]=0;
}
LL PowMod(LL x,LL y)
{
LL ret=1;
while(y)
{
if(y&1)
ret=1LL*ret*x%MO;
x=1LL*x*x%MO;
y>>=1;
}
return ret;
}
void NTT(LL P[],LL len,LL oper)
{
for(LL i=1,j=0;i<len-1;i++)
{
for(LL s=len;j^=s>>=1,~j&s;);
if(i<j) swap(P[i],P[j]);
}
LL unit,unit_p0;
for(LL d=0;(1<<d)<len;d++)
{
LL m=(1<<d),m2=m*2;
unit_p0=PowMod(G,(MO-1)/m2);
if(oper==-1)
unit_p0=PowMod(unit_p0,MO-2);
for(LL i=0;i<len;i+=m2)
{
unit=1;
for(LL j=0;j<m;j++)
{
LL &P1=P[i+j+m],&P2=P[i+j];
LL t=1LL*unit*P1%MO;
P1=((1LL*P2-1LL*t)%MO+1LL*MO)%MO;
P2=(1LL*P2+1LL*t)%MO;
unit=1LL*unit*unit_p0%MO;
}
}
}
if(oper==-1)
{
LL INV=PowMod(len,MO-2);
for(LL i=0;i<len;i++)
P[i]=1LL*P[i]*INV%MO;
}
}
void Polynomial_Inverse(LL deg,LL A[],LL B[])
{
static LL tmpA[MAXN+5],tmpB[MAXN+5];
if(deg==1)
{
B[0]=PowMod(A[0],MO-2);
return;
}
Polynomial_Inverse((deg+1)/2,A,B);
LL p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);
copy(B,B+deg,tmpB);
fill(tmpA+deg,tmpA+p,0);
fill(tmpB+deg,tmpB+p,0);
NTT(tmpA,p,1);NTT(tmpB,p,1);
for(LL i=0;i<p;i++)
tmpB[i]=(1LL*tmpB[i]*(2LL-1LL*tmpB[i]*tmpA[i]%MO)%MO+1LL*MO)%MO;
NTT(tmpB,p,-1);
copy(tmpB,tmpB+deg,B);
}
void Mul(LL ret[],LL _x[],LL l1,LL _y[],LL l2)
{
static LL X[MAXN+5],Y[MAXN+5],RET[MAXN+5];
LL len=1;
while(len<l1+l2) len<<=1;
copy(_x,_x+l1,X);copy(_y,_y+l2,Y);
fill(X+l1,X+len,0);fill(Y+l2,Y+len,0);
NTT(X,len,1);NTT(Y,len,1);
for(LL i=0;i<len;i++)
RET[i]=1LL*X[i]*Y[i]%MO;
NTT(RET,len,-1);
copy(RET,RET+l1+l2,ret);
}
void Polynomial_Logarithm(LL deg,LL A[],LL B[])
{
static LL tmpA[MAXN+5],tmpA2[MAXN+5];
LL p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);
fill(tmpA+deg,tmpA+p,0);
fill(tmpA2,tmpA2+p,0);
Polynomial_Inverse(deg,tmpA,tmpA2);
Diff(tmpA,deg);
Mul(tmpA,tmpA,deg-1,tmpA2,deg);
Calc(tmpA,deg);
copy(tmpA,tmpA+deg,B);
}
void Polynomial_Exponential(LL deg,LL A[],LL B[])
{
static LL tmpA[MAXN+5],tmpB[MAXN+5],tmpB2[MAXN+5];
if(deg==1)
{
B[0]=1;
return;
}
Polynomial_Exponential((deg+1)/2,A,B);
LL p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);
copy(B,B+deg,tmpB);
fill(tmpA+deg,tmpA+p,0);
fill(tmpB+deg,tmpB+p,0);
Polynomial_Logarithm(deg,tmpB,tmpB2);
tmpA[0]=(1LL*tmpA[0]+1LL)%MO;
for(LL i=0;i<deg;i++)
tmpA[i]=((1LL*tmpA[i]-1LL*tmpB2[i])%MO+1LL*MO)%MO;
NTT(tmpA,p,1);NTT(tmpB,p,1);
for(LL i=0;i<p;i++)
tmpB[i]=1LL*tmpA[i]*tmpB[i]%MO;
NTT(tmpB,p,-1);
copy(tmpB,tmpB+deg,B);
fill(B+deg,B+p,0);
}
int main()
{
Init();
scanf("%lld",&n);
for(LL i=0;i<n;i++)
scanf("%lld",&a[i]);
Polynomial_Exponential(n,a,b);
for(LL i=0;i<n;i++)
printf("%lld ",b[i]);
printf("\n");
return 0;
}
```
#多项式开根
直接上问题：
$$\sqrt{A(x)}\equiv B(x)(mod\ x^{\frac{n}{2}})\]
其中A(x)已知，B(x)是要求的。
把根号去掉就有：
\[A(x)\equiv B^2(x)(mod\ x^n)
\]
然后又用常见的套路：
\[G(F(x))=A(x)-(B(x))^2
\]
就有牛顿迭代：
\[B(x)=B_0(x)-\frac{A(x)-(B_0(x))^2}{-2*B_0(x)}
\]
化简一下就有：
\[B(x)=\frac{B_0(x)}{2}+\frac{A(x)}{2*B_0(x)}
\]
然后就可以做了(求逆就可以了)。
多项式求幂
这个就比较简单（但是也想不到这个思路）。
之前我们就有一个比较简单的一个思路，就是多项式快速幂，是\(O(n\log^2 n)\)的算法。这里利用之前的方法能够优化到\(O(n\log n)\)（但是常数巨大...）。
思路其实就是将\((F(x))^k\)转化为\(e^{k*ln(F(x))}\)。
如果觉得不显然的话，展开之后，你就会发现他们是一样的!!!
那就是nlogn了（是不是感觉有点震惊...）。
组合板题：帕秋莉的超级多项式
（多项式运算五合一）
代码：
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define MAXN 500000
#define MO 998244353
#define G 3
using namespace std;
int a[MAXN+5],b[MAXN+5],inv[MAXN+5];
void Prepare()
{
inv[1]=1;
for(int i=2;i<=MAXN;i++)
inv[i]=1LL*inv[MO%i]*(1LL*MO-1LL*MO/i)%MO;
}
int PowMod(int x,int y)
{
int ret=1;
while(y)
{
if(y&1)
ret=1LL*ret*x%MO;
y>>=1;
x=1LL*x*x%MO;
}
return ret;
}
void Diff(int A[],int len)
{
for(int i=0;i<len-1;i++)
A[i]=1LL*A[i+1]*(1LL*i+1LL)%MO;
A[len-1]=0;
}
void Calc(int A[],int len)
{
for(int i=len-1;i>=1;i--)
A[i]=1LL*A[i-1]*inv[i]%MO;
A[0]=0;
}
void NTT(int P[],int len,int oper)
{
for(int i=1,j=0;i<len-1;i++)
{
for(int s=len;j^=s>>=1,~j&s;);
if(i<j) swap(P[i],P[j]);
}
int unit,unit_p0;
for(int d=0;(1<<d)<len;d++)
{
int m=(1<<d),m2=m*2;
unit_p0=PowMod(G,(MO-1)/m2);
if(oper==-1)
unit_p0=PowMod(unit_p0,MO-2);
for(int i=0;i<len;i+=m2)
{
unit=1;
for(int j=0;j<m;j++)
{
int &P1=P[i+j+m],&P2=P[i+j];
int t=1LL*unit*P1%MO;
P1=((1LL*P2-1LL*t)%MO+MO)%MO;
P2=(1LL*P2+1LL*t)%MO;
unit=1LL*unit*unit_p0%MO;
}
}
}
if(oper==-1)
{
int INV=PowMod(len,MO-2);
for(int i=0;i<len;i++)
P[i]=1LL*P[i]*INV%MO;
}
}
void Mul(int ret[],int _x[],int l1,int _y[],int l2)
{
static int RET[MAXN+5],X[MAXN+5],Y[MAXN+5];
int len=1;
while(len<l1+l2) len<<=1;
copy(_x,_x+l1,X);copy(_y,_y+l2,Y);
fill(X+l1,X+len,0);fill(Y+l2,Y+len,0);
NTT(X,len,1);NTT(Y,len,1);
for(int i=0;i<len;i++)
RET[i]=1LL*X[i]*Y[i]%MO;
NTT(RET,len,-1);
copy(RET,RET+l1+l2-1,ret);
}
void Polynomial_Inverse(int deg,int A[],int B[])
{
static int tmpA[MAXN+5],tmpB[MAXN+5];
if(deg==1)
{
B[0]=PowMod(A[0],MO-2);
return;
}
Polynomial_Inverse((deg+1)/2,A,B);
int p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);copy(B,B+deg,tmpB);
fill(tmpA+deg,tmpA+p,0);fill(tmpB+deg,tmpB+p,0);
NTT(tmpA,p,1);NTT(tmpB,p,1);
for(int i=0;i<p;i++)
tmpB[i]=(1LL*tmpB[i]*(2LL-1LL*tmpA[i]*tmpB[i]%MO)%MO+MO)%MO;
NTT(tmpB,p,-1);
copy(tmpB,tmpB+deg,B);
fill(B+deg,B+p,0);
}
void Polynomial_Logarithm(int deg,int A[],int B[])
{
static int tmpA[MAXN+5],tmpA2[MAXN+5];
int p=1;
while(p<2*deg) p<<=1;
copy(A,A+deg,tmpA);
fill(tmpA+deg,tmpA+p,0);
fill(tmpA2,tmpA2+p,0);
Polynomial_Inverse(deg,tmpA,tmpA2);
Diff(tmpA,deg);
Mul(tmpA,tmpA,deg-1,tmpA2,deg);
Calc(tmpA,deg);
copy(tmpA,tmpA+deg,B);
}
void Polynomial_Exponential(int deg,int A[],int B[])
{
static int tmpA[MAXN+5],tmpB[MAXN+5],tmpB2[MAXN+5];
if(deg==1)
{
B[0]=1;
return;
}
Polynomial_Exponential((deg+1)/2,A,B);
int p=1;
while(p<2*deg) p<<=1;
copy(A,A+deg,tmpA);copy(B,B+deg,tmpB);
fill(tmpA+deg,tmpA+p,0);fill(tmpB+deg,tmpB+p,0);
fill(tmpB2,tmpB2+p,0);
Polynomial_Logarithm(deg,tmpB,tmpB2);
tmpA[0]=(1LL*tmpA[0]+1LL)%MO;
for(int i=0;i<deg;i++)
tmpA[i]=((1LL*tmpA[i]-1LL*tmpB2[i])%MO+1LL*MO)%MO;
NTT(tmpA,p,1);NTT(tmpB,p,1);
for(int i=0;i<p;i++)
tmpB[i]=1LL*tmpA[i]*tmpB[i]%MO;
NTT(tmpB,p,-1);
copy(tmpB,tmpB+deg,B);
fill(B+deg,B+p,0);
}
void Polynomial_Sqrt(int deg,int A[],int B[])
{
static int tmpA[MAXN+5],tmpB[MAXN+5],tmpB2[MAXN+5];
if(deg==1)
{
B[0]=sqrt(A[0]);
return;
}
Polynomial_Sqrt((deg+1)/2,A,B);
int p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);
copy(B,B+deg,tmpB);
fill(tmpA+deg,tmpA+p,0);
fill(tmpB+deg,tmpB+p,0);
fill(tmpB2,tmpB2+p,0);
Polynomial_Inverse(deg,tmpB,tmpB2);
NTT(tmpB,p,1);NTT(tmpB2,p,1);NTT(tmpA,p,1);
for(int i=0;i<p;i++)
tmpB[i]=(1LL*tmpB[i]+1LL*tmpA[i]*tmpB2[i]%MO)%MO*inv[2]%MO;
NTT(tmpB,p,-1);
copy(tmpB,tmpB+deg,B);
fill(B+deg,B+p,0);
}
void Polynomial_Pow(int deg,int K,int A[],int B[])
{
static int tmpA[MAXN+5],tmpA2[MAXN+5];
int p=1;
while(p<deg*2) p<<=1;
copy(A,A+deg,tmpA);
fill(tmpA+deg,tmpA+p,0);
fill(tmpA2,tmpA2+p,0);
Polynomial_Logarithm(deg,tmpA,tmpA2);
for(int i=0;i<deg;i++)
tmpA2[i]=1LL*K*tmpA2[i]%MO;
fill(tmpA,tmpA+p,0);
Polynomial_Exponential(deg,tmpA2,tmpA);
copy(tmpA,tmpA+deg,B);
fill(B+deg,B+p,0);
}
int main()
{
Prepare();
int n,k;
scanf("%d %d",&n,&k);
for(int i=0;i<n;i++)
scanf("%d",&a[i]);
Polynomial_Sqrt(n,a,b);
memset(a,0,sizeof(a));
Polynomial_Inverse(n,b,a);
Calc(a,n);
memset(b,0,sizeof(b));
Polynomial_Exponential(n,a,b);
memset(a,0,sizeof(a));
Polynomial_Inverse(n,b,a);
a[0]=(1LL*a[0]+1LL)%MO;
memset(b,0,sizeof(b));
Polynomial_Logarithm(n,a,b);
b[0]=(1LL*b[0]+1LL)%MO;
memset(a,0,sizeof(a));
Polynomial_Pow(n,k,b,a);
Diff(a,n);
for(int i=0;i<n;i++)
if(i==0)
printf("%d",a[i]);
else
printf(" %d",a[i]);
return 0;
}
