多项式的ln、exp、快速幂和开根学习小记
-
不妨又學習了一下多項式的求ln、exp、快速冪和開根操作。
-
這些操作比之前的求逆更上了一層臺階,應用同樣很廣。
-
多項式求逆等知識在我的博客里有講:多項式的求逆、取模和多點求值學習小記
多項式ln
-
給出一個次數(shù)界為 nnn 的多項式 F(x)F(x)F(x) ,要求多項式 G(x)G(x)G(x) 滿足:G(x)≡lnF(x)(modxn)G(x)\equiv ln\ F(x)\ (mod\ x^n)G(x)≡ln?F(x)?(mod?xn)
-
我們對其兩邊求導(右邊相當于是復合函數(shù)求導):G′(x)=F′(x)F(x)G'(x)=\frac{F'(x)}{F(x)}G′(x)=F(x)F′(x)?
-
再兩邊積分,即:G(x)=∫F′(x)F(x)dxG(x)=\int\frac{F'(x)}{F(x)}\ dxG(x)=∫F(x)F′(x)??dx
-
于是應用一下多項式求逆即可在 O(nlogn)O(n\ log\ n)O(n?log?n) 的時間內(nèi)求出一個多項式的 ln 了。
-
其中求導和積分都是 O(n)O(n)O(n) ,具體來說就是:(xa)′=axa?1(x^a)'=ax^{a-1}(xa)′=axa?1∫xadx=1a+1xa+1\int{x^a}\ dx=\frac{1}{a+1}x^{a+1}∫xa?dx=a+11?xa+1
-
逐項操作即可求導、積分。
多項式exp
-
給出一個次數(shù)界為 nnn 的多項式 F(x)F(x)F(x) ,要求多項式 G(x)G(x)G(x) 滿足:G(x)≡eF(x)(modxn)G(x)\equiv e^{F(x)}\ (mod\ x^n)G(x)≡eF(x)?(mod?xn)
-
這個需要用到 牛頓迭代法——一種在實數(shù)域和復數(shù)域上近似求解方程的方法。
-
牛頓迭代法的一般式如下(要求的根為 xxx):x′=x?f(x)f′(x)x'=x-\frac{f(x)}{f'(x)}x′=x?f′(x)f(x)?
-
于是我們先將定義式移項:G(x)?eF(x)=0G(x)-e^{F(x)}=0G(x)?eF(x)=0
-
兩邊取對數(shù):lnG(x)?F(x)=0ln\ G(x)-F(x)=0ln?G(x)?F(x)=0
-
把 G(x)G(x)G(x) 當做要求的方程的根,則根據(jù)牛頓迭代法可得:Gt+1(x)=Gt(x)?lnGt(x)?F(x)1Gt(x)=Gt(x)(1+F(x)?lnGt(x))G_{t+1}(x)=G_t(x)-\frac{ln\ G_t(x)-F(x)}{\frac{1}{G_t(x)}}=G_t(x)(1+F(x)-ln\ G_t(x))Gt+1?(x)=Gt?(x)?Gt?(x)1?ln?Gt?(x)?F(x)?=Gt?(x)(1+F(x)?ln?Gt?(x))
-
這樣每次迭代都使精度(次數(shù)界)翻倍,遞歸處理即可,每次遞歸過程中要進行一次多項式求ln 。
-
這樣的時間復雜度是 T(n)=T(n2)+O(nlogn)=O(nlogn)T(n)=T(\frac{n}{2})+O(n\ log\ n)=O(n\ log\ n)T(n)=T(2n?)+O(n?log?n)=O(n?log?n) 的,不過常數(shù)較大。
多項式快速冪
-
給出一個次數(shù)界為 nnn 的多項式 F(x)F(x)F(x) ,要求多項式 G(x)G(x)G(x) 滿足:G(x)≡F(x)k(modxn)G(x)\equiv F(x)^k\ (mod\ x^n)G(x)≡F(x)k?(mod?xn)
-
如果像普通一樣每次直接 NTT 的話,復雜度會很不理想,特別是當 kkk 大到一定程度的時候。
-
我們可以兩邊取對數(shù),則有:lnG(x)=lnF(x)k=klnF(x)ln\ G(x)=ln\ F(x)^k=k\ ln\ F(x)ln?G(x)=ln?F(x)k=k?ln?F(x)
-
右邊多項式系數(shù)直接乘 kkk ,再兩邊exp還原即可。
-
這樣做的話 kkk 的值可以很大,我們只需要知道它模意義下的值即可。
-
這樣多項式快速冪需要用到求ln、exp,但復雜度仍是 O(nlogn)O(n\ log\ n)O(n?log?n) 的。
-
附模板題和我的代碼: 洛谷 P5245 【模板】多項式快速冪
-
值得注意的是這里的 kkk 達到了 1010510^{10^5}10105 數(shù)量級,我們用高精度除單精度求其模意義下的值即可。
Code
#include<cstdio> #include<cstring> #include<algorithm> #include<cctype> using namespace std; typedef long long LL; const int N=1e5+5,L=N<<2,G=3,mo=998244353; int a[L],b[L],c[L]; int f[L],rev[L],wn[L],inv[L]; char s[N]; inline int read() {int X=0,w=0; char ch=0;while(!isdigit(ch)) w|=ch=='-',ch=getchar();while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();return w?-X:X; } inline int ksm(int x,int y) {int s=1;while(y){if(y&1) s=(LL)s*x%mo;x=(LL)x*x%mo;y>>=1;}return s; } inline void NTT(int *y,int len,int ff) {for(int i=0;i<len;i++)if(i<rev[i]) swap(y[i],y[rev[i]]);for(int h=2,d=len>>1;h<=len;h<<=1,d>>=1)for(int i=0,k=h>>1;i<len;i+=h)for(int j=0,cnt=0;j<k;j++,cnt+=d){int u=y[i+j],t=(LL)y[i+j+k]*wn[cnt]%mo;y[i+j]=u+t>=mo?u+t-mo:u+t;y[i+j+k]=u-t<0?u-t+mo:u-t;}if(ff==-1){reverse(y+1,y+len);int invl=inv[len];for(int i=0;i<len;i++) y[i]=(LL)y[i]*invl%mo;} } void getinv(int len,int num) {if(len==1){b[0]=ksm(c[0],mo-2);return;}getinv(len>>1,num-1);for(int i=0;i<len>>1;i++) f[i]=c[i];for(int i=len>>1;i<len;i++) f[i]=0;for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<num-1;int w0=ksm(G,(mo-1)/len);for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;NTT(f,len,1),NTT(b,len,1);for(int i=0;i<len;i++) b[i]=(2-(LL)b[i]*f[i]%mo+mo)*b[i]%mo;NTT(b,len,-1);for(int i=len>>1;i<len;i++) b[i]=0; } inline void getln(int len,int num) {getinv(len,num);for(int i=0;i<len-1;i++) f[i]=(LL)c[i+1]*(i+1)%mo;f[len-1]=0;for(int i=len>>1;i<len;i++) f[i]=b[i]=0;NTT(f,len,1),NTT(b,len,1);for(int i=0;i<len;i++) f[i]=(LL)f[i]*b[i]%mo;NTT(f,len,-1);for(int i=len-1;i>0;i--) f[i]=(LL)f[i-1]*inv[i]%mo;f[0]=0; } void getexp(int len,int num) {if(len==1){c[0]=1;return;}getexp(len>>1,num-1);for(int i=0;i<len;i++) f[i]=b[i]=0;getln(len,num);f[0]=(a[0]+1-f[0]+mo)%mo;for(int i=1;i<len>>1;i++) f[i]=(a[i]-f[i]+mo)%mo;for(int i=len>>1;i<len;i++) f[i]=0;NTT(c,len,1),NTT(f,len,1);for(int i=0;i<len;i++) c[i]=(LL)c[i]*f[i]%mo;NTT(c,len,-1);for(int i=len>>1;i<len;i++) c[i]=0; } int main() {int n=read(),k=0;scanf("%s",s+1);int m=strlen(s+1);for(int i=1;i<=m;i++) k=((LL)k*10+s[i]-'0')%mo;for(int i=0;i<n;i++) c[i]=read();inv[0]=inv[1]=1;for(int i=2;i<L;i++) inv[i]=(LL)(mo-mo/i)*inv[mo%i]%mo;int len=1,num=0;while(len<n<<1) len<<=1,num++;getln(len,num);for(int i=0;i<n;i++) a[i]=(LL)f[i]*k%mo;for(int i=0;i<len;i++) f[i]=b[i]=c[i]=0;getexp(len,num);for(int i=0;i<n;i++) printf("%d ",c[i]);return 0; }多項式開根
-
給出一個次數(shù)界為 nnn 多項式 A(x)A(x)A(x) ,要求一個次數(shù)界為 nnn 的多項式 B(x)B(x)B(x) 滿足:B2(x)≡A(x)(modxn)(?)B^2(x)\equiv A(x)\ (mod\ x^n)\ \ \ \ (*)B2(x)≡A(x)?(mod?xn)????(?)
-
這個在做生成函數(shù)題時用得很多,放在這里是因為這可以用 ln、exp 來處理(需要保證 A0=1A_0=1A0?=1)。
-
兩邊取對數(shù),得:lnB(x)=12lnA(x)ln\ B(x)=\frac{1}{2}\ ln\ A(x)ln?B(x)=21??ln?A(x)
-
再兩邊 exp 還原即可,本質與多項式快速冪是一樣的,時間復雜度 O(nlogn)O(n\ log\ n)O(n?log?n) 。
-
當然開根也有其獨特的解法,跟求逆一樣是倍增來求的。
-
我們已求出 G(x)G(x)G(x) 滿足:G2(x)≡A(x)(modx?n2?)G^2(x)\equiv A(x)\ (mod\ x^{\lfloor\frac{n}{2}\rfloor})G2(x)≡A(x)?(mod?x?2n??)
-
用 (*) 式減上式,得:B2(x)?G2(x)≡0(modx?n2?)B^2(x)-G^2(x)\equiv0\ (mod\ x^{\lfloor\frac{n}{2}\rfloor})B2(x)?G2(x)≡0?(mod?x?2n??)
-
平方差公式:(B(x)+G(x))(B(x)?G(x))≡0(modx?n2?)(B(x)+G(x))(B(x)-G(x))\equiv0\ (mod\ x^{\lfloor\frac{n}{2}\rfloor})(B(x)+G(x))(B(x)?G(x))≡0?(mod?x?2n??)
-
這里的根有兩個,這里考慮其中一個:B(x)?G(x)≡0(modx?n2?)B(x)-G(x)\equiv0\ (mod\ x^{\lfloor\frac{n}{2}\rfloor})B(x)?G(x)≡0?(mod?x?2n??)
-
兩邊平方(次數(shù)界擴大了兩倍):(B(x)?G(x))2≡0(modxn)(B(x)-G(x))^2\equiv0\ (mod\ x^n)(B(x)?G(x))2≡0?(mod?xn)
-
展開來:B2(x)?2B(x)G(x)+G2(x)≡0(modxn)B^2(x)-2B(x)G(x)+G^2(x)\equiv0\ (mod\ x^n)B2(x)?2B(x)G(x)+G2(x)≡0?(mod?xn)
-
消去 B2(x)B^2(x)B2(x) ,即:A(x)?2B(x)G(x)+G2(x)≡0(modxn)A(x)-2B(x)G(x)+G^2(x)\equiv0\ (mod\ x^n)A(x)?2B(x)G(x)+G2(x)≡0?(mod?xn)
-
移項即得:B(x)≡A(x)+G2(x)2G(x)(modxn)B(x)\equiv\frac{A(x)+G^2(x)}{2G(x)}\ (mod\ x^n)B(x)≡2G(x)A(x)+G2(x)??(mod?xn)
-
或者這種形式:B(x)≡12(A(x)G(x)+G(x))(modxn)B(x)\equiv\frac{1}{2}(\frac{A(x)}{G(x)}+G(x))\ (mod\ x^n)B(x)≡21?(G(x)A(x)?+G(x))?(mod?xn)
-
于是我們多項式求逆并計算即可,時間復雜度 O(nlogn)O(n\ log\ n)O(n?log?n) 。
-
注意遞歸底層是 B0=A0B_0=\sqrt {A_0}B0?=A0?? ,若是 0、1 就最好了,否則我們得算一下二次剩余(當然必須得有,不然沒法開根)!
-
如果題目使得這個值是固定的,那么就可以直接枚舉算二次剩余是啥。
-
如果是會變的,也許我們要用用什么 CipollaCipollaCipolla 算法……這里就不做討論了。。
-
有興趣的同學可以看看 a_crazy_czy大神 的 二次剩余Cipolla算法學習小記 。
-
模板題:洛谷 P5205 【模板】多項式開根
Code
#include<cstdio> #include<cstring> #include<algorithm> #include<cctype> using namespace std; typedef long long LL; const int N=1e5+5; int G=3,mo=998244353,inv2=mo+1>>1; int a[N<<2],b[N<<2],c[N<<2]; int f[N<<2],rev[N<<2],wn[N<<2]; inline int read() {int X=0,w=0; char ch=0;while(!isdigit(ch)) w|=ch=='-',ch=getchar();while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();return w?-X:X; } void write(int x) {if(x>9) write(x/10);putchar(x%10+'0'); } inline int ksm(int x,int y) {int s=1;while(y){if(y&1) s=(LL)s*x%mo;x=(LL)x*x%mo;y>>=1;}return s; } inline void NTT(int *y,int len,int ff) {for(int i=0;i<len;i++)if(i<rev[i]) swap(y[i],y[rev[i]]);for(int h=2,d=len>>1;h<=len;h<<=1,d>>=1)for(int i=0,k=h>>1;i<len;i+=h)for(int j=0,cnt=0;j<k;j++,cnt+=d){int u=y[i+j],t=(LL)wn[cnt]*y[i+j+k]%mo;y[i+j]=u+t>=mo?u+t-mo:u+t;y[i+j+k]=u-t<0?u-t+mo:u-t;}if(ff==-1){for(int i=len>>1;i;i--) swap(y[i],y[len-i]);int inv=ksm(len,mo-2);for(int i=0;i<len;i++) y[i]=(LL)y[i]*inv%mo;} } void solve1(int len,int num) {if(len==1){c[0]=ksm(b[0],mo-2);return;}solve1(len>>1,num-1);for(int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|(i&1)<<num-1;int w0=ksm(G,(mo-1)/len);for(int i=wn[0]=1;i<=len;i++) wn[i]=(LL)wn[i-1]*w0%mo;for(int i=0;i<len>>1;i++) f[i]=b[i];for(int i=len>>1;i<len;i++) f[i]=0;NTT(f,len,1),NTT(c,len,1);for(int i=0;i<len;i++) f[i]=(2-(LL)f[i]*c[i]%mo+mo)*c[i]%mo;NTT(f,len,-1);for(int i=0;i<len>>1;i++) c[i]=f[i];for(int i=len>>1;i<len;i++) c[i]=0; } void solve(int len,int num) {if(len==1){b[0]=1;return;}solve(len>>1,num-1);solve1(len,num);for(int i=0;i<len>>1;i++) f[i]=a[i];for(int i=len>>1;i<len;i++) f[i]=0;NTT(f,len,1),NTT(c,len,1);for(int i=0;i<len;i++) f[i]=(LL)f[i]*c[i]%mo;NTT(f,len,-1);for(int i=0;i<len>>1;i++) b[i]=(LL)(f[i]+b[i])*inv2%mo;for(int i=len>>1;i<len;i++) b[i]=0;for(int i=0;i<len;i++) c[i]=0; } int main() {int n=read();for(int i=0;i<n;i++) a[i]=read();int mx=1,ll=0;while(mx<n+n) mx<<=1,ll++;solve(mx,ll);for(int i=0;i<n;i++) write(b[i]),putchar(' ');return 0; }總結
以上是生活随笔為你收集整理的多项式的ln、exp、快速幂和开根学习小记的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: C++ 40行超级加速命令
- 下一篇: 多项式快速插值学习小记