參考資料
牛頓迭代法在多項式運算的應用-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;
}