【学习笔记】之多项式使人头秃
真的自閉= =
多項式是什么鬼哦
?
首先 介紹 FFT
我才不想寫那么多柿子呢
大體說一下FFT干了啥
我們對兩個多項式進行卷積(即多項式乘法)
?也就是?
暴力計算的話是n^2的
我們考慮把它變成點值[即(x,y)表示f(x)=y] 點值相乘就快了嘛 但是變成點值了以后咋變回來呢
有個叫傅里葉的nb的人 他發(fā)明了一個nb的東西叫傅里葉變換= =
也就是通過 虛數(shù)中的單位根 來計算就可以變回來了
單位根是什么東西呢 就是在復平面上的一個單位圓 將其弧等分成若干份 第一個點位于(0,1)的n個點 把這些數(shù)帶進去就可以做啦
說起來很奇特對不對 他其實就很奇特= =
具體詳細的東西右轉(zhuǎn)百度吧 我實在是懶得寫QAQ
?實現(xiàn)上可以直接使用模板庫里的complex(雖然我用起來非常不習慣)
扔個代碼跑路。
#include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define maxn 3000010 using namespace std;const double Pi = acos(-1.0);struct complex {double x,y;complex(double xx=0.0,double yy=0.0){x=xx,y=yy;} }A[maxn],B[maxn]; int l,r[maxn],limit=1; complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);} complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);} complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}void FFT(complex *a,int type) {for(int i=0;i<limit;i++)if(i<r[i]) swap(a[i],a[r[i]]);for(int mid=1;mid<limit;mid<<=1){complex Wn=complex(cos(Pi/mid),type*sin(Pi/mid));for(int R=mid<<1,j=0;j<limit;j+=R){complex w=complex(1,0);for(int k=0;k<mid;k++,w=w*Wn){complex x=a[j+k],y=w*a[j+mid+k];a[j+k]=x+y;a[j+mid+k]=x-y;}}} }int main() {int N,M,i,j;scanf("%d%d",&N,&M);for(i=0;i<=N;i++) scanf("%lf",&A[i].x);for(i=0;i<=M;i++) scanf("%lf",&B[i].x);while(limit<=(N+M)) limit<<=1,l++;for(i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));FFT(A,1);FFT(B,1);for(i=0;i<=limit;i++)A[i]=A[i]*B[i];FFT(A,-1);for(i=0;i<=N+M;i++)printf("%d ",(int)(A[i].x/limit+0.5));return 0; }?
然后我們就遇到了一個神奇的模數(shù) 998244353 才不是1XXXXXX7
為什么是這個模數(shù)呢 因為他是一個2^x* ...+1的一個素數(shù) 具有一些優(yōu)美的性質(zhì)
我們就可以進行NTT[快速數(shù)論變換] /斜眼笑/
我們剛剛FFT中用的復平面中的單位根 所以是有小數(shù)的
這個樣子可不大好因為我們要取模 所以我們有了一個很nb的東西叫做 原根
原根有一些優(yōu)美的性質(zhì) 就跟單位根一樣 G ^ ((p-1)/i) 就是可以當成單位根使用的數(shù)啦 [i|(p-1)]這也就是p為什么要是 2^x *... +1的原因啦?
小姿勢:998244353的原根是3
其他詳細的細節(jié)還是右轉(zhuǎn)百度吧【大霧
扔個代碼。
#include<cstdio> #include<algorithm> #include<cstring> #define maxn 300005 #define modn 998244353 #define G 3 #define ll long long using namespace std;int q_pow(ll base,ll pow) {ll ans=1;while(pow){if(pow&1){ans*=base;ans%=modn;}base*=base;base%=modn;pow>>=1;}return (int)ans; }int A[maxn],B[maxn],C[maxn]; int l,r[maxn],limit,inv; void FFT(int *a,int type) {int i,j,k;for(i=0;i<limit;i++)if(r[i]>i)swap(a[r[i]],a[i]);for(i=2;i<=limit;i<<=1){int mid=i>>1;int Wn=q_pow(G,(modn-1)/i);if(type) Wn=q_pow(Wn,(modn-2));for(j=0;j<limit;j+=i){int w=1;for(k=0;k<mid;k++){int x=a[j+k],y=a[j+mid+k];a[j+k]=x+(ll)w*y%modn;if(a[j+k]>=modn) a[j+k]-=modn;a[j+k+mid]=x-(ll)w*y%modn;if(a[j+mid+k]<0) a[j+mid+k]+=modn;w=(ll)w*Wn%modn;}}}if(type){for(i=0;i<limit;i++)a[i]=(ll)a[i]*inv%modn;} }int main() {int n,i,j,k,m,s;scanf("%d%d",&n,&m);for(i=0;i<=n;i++) scanf("%d",&A[i]);for(i=0;i<=m;i++) scanf("%d",&B[i]);limit=1;while(limit<=(n+m)) limit<<=1,l++;for(i=0;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));inv=q_pow(limit,modn-2);FFT(A,0);FFT(B,0);for(i=0;i<=limit;i++)C[i]=(ll)A[i]*B[i]%modn;FFT(C,1);for(i=0;i<=(n+m);i++)printf("%d\n",C[i]);return 0; }?
于是我的姿勢還只停留在 FFT/NTT 只是能求個卷積【大霧
?
然后就有一些奇奇怪怪的東西了=.=+
?
【奇奇怪怪的東西一】多項式求逆
我們現(xiàn)在有一個多項式f 我們要求一個多項式g滿足?
這玩意看上去是不是非常奇怪啊【明明就是非常奇怪!
假設我們現(xiàn)在已知一個多項式h 滿足?
我們可以得到
平方
卷上f?
移項
遞歸求解就好啦
邊界當然是n=1的時候 g直接取f的常數(shù)項的逆元啦qwq。
附代碼。
#include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define inf 20021225 #define ll long long #define mdn 998244353 #define G 3 #define mxn 300100 using namespace std;int rev[mxn];int ksm(int bs,int mi) {int ans=1;while(mi){if(mi&1) ans=(ll)ans*bs%mdn;bs=(ll)bs*bs%mdn;mi>>=1;}return ans; } int inv; void NTT(int *a,int lim,int f) {for(int i=0;i<lim;i++)if(rev[i]>i) swap(a[i],a[rev[i]]);for(int k=2;k<=lim;k<<=1){int Wn=ksm(G,(mdn-1)/k),mid=k>>1;if(f) Wn=ksm(Wn,mdn-2);for(int w=1,i=0;i<lim;i+=k,w=1){for(int j=0;j<mid;j++,w=(ll)w*Wn%mdn){int x=a[i+j],y=(ll)w*a[i+j+mid]%mdn;a[i+j]=(x+y)%mdn;a[i+j+mid]=(x-y+mdn)%mdn;}}}if(f) for(int i=0;i<lim;i++) a[i]=(ll)a[i]*inv%mdn; } int g[mxn],h[mxn],f[mxn];void poly_inv(int n) {if(n==1){g[0]=ksm(f[0],mdn-2);//printf("%d\n",g[0]);return;}poly_inv((n+1)>>1);int lim=1,l=0;while(lim<(n<<1)) lim<<=1,l++;for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));inv=ksm(lim,mdn-2);for(int i=0;i<n;i++) h[i]=f[i];for(int i=n;i<lim;i++) h[i]=0;NTT(h,lim,0);NTT(g,lim,0);for(int i=0;i<lim;i++)g[i]=(ll)(2ll-(ll)g[i]*h[i]%mdn+mdn)%mdn*g[i]%mdn;NTT(g,lim,1);for(int i=n;i<lim;i++) g[i]=0; }int main() {int n;scanf("%d",&n);for(int i=0;i<n;i++) scanf("%d",&f[i]);poly_inv(n);for(int i=0;i<n;i++) printf("%d ",g[i]);return 0; }?
【奇奇怪怪的東西二】多項式對數(shù)函數(shù)
看上去是不是很高大上!【實則蠢得一批
對于多項式 f 求g=ln f
這之前我們科普一點求導和積分的小姿勢
對于一個普通多項式
求導
積分
兩個過程都很像噠 就是反過來做而已233
ln x求導是 1/x
復合函數(shù)求導
然后ln f的求導
直接多項式求逆然后求導再積分回去就好啦qwq
代碼等我一哈【不咕不咕必定不可能咕
update:真的沒有咕!
#include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define inf 20021225 #define ll long long #define mdn 998244353 #define G 3 #define mxn 300100 using namespace std;int ksm(int bs,int mi) {int ans=1;while(mi){if(mi&1) ans=(ll)ans*bs%mdn;bs=(ll)bs*bs%mdn;mi>>=1;}return ans; } int inv,rev[mxn]; void NTT(int *a,int lim,int f) {for(int i=0;i<lim;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);for(int k=2;k<=lim;k<<=1){int mid=k>>1,Wn=ksm(G,(mdn-1)/k);if(f) Wn=ksm(Wn,mdn-2);for(int w=1,i=0;i<lim;i+=k,w=1)for(int j=0;j<mid;j++,w=(ll)w*Wn%mdn){int x=a[i+j],y=(ll)w*a[i+j+mid]%mdn;a[i+j]=(x+y)%mdn;a[i+j+mid]=(x-y+mdn)%mdn;}}if(f) for(int i=0;i<lim;i++) a[i]=(ll)a[i]*inv%mdn; }void init(int lim,int l) {for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);inv=ksm(lim,mdn-2); }int h[mxn],g[mxn],f[mxn];void poly_inv(int n) {if(n==1){g[0]=ksm(f[0],mdn-2);return;}poly_inv((n+1)>>1);int lim=1,l=0;while(lim<(n<<1)) lim<<=1,l++;for(int i=0;i<n;i++) h[i]=f[i];for(int i=n;i<lim;i++) h[i]=0;init(lim,l);NTT(h,lim,0);NTT(g,lim,0);for(int i=0;i<lim;i++)g[i]=(2ll -(ll)g[i] *h[i] %mdn +mdn) %mdn *g[i]%mdn;NTT(g,lim,1);for(int i=n;i<lim;i++) g[i]=0; } int d[mxn]; void poly_ln(int n) {poly_inv(n);for(int i=0;i<n;i++) d[i] =(ll)f[i+1] * (i+1) %mdn;int l=0,lim=1;while(lim<(n<<1)) lim<<=1,l++;init(lim,l);NTT(d,lim,0);NTT(g,lim,0);for(int i=0;i<lim;i++) g[i]=(ll)g[i]*d[i]%mdn;NTT(g,lim,-1);for(int i=0;i<n;i++) d[i+1]=(ll)g[i]*ksm(i+1,mdn-2)%mdn;d[0]=0; }int main() {int n;scanf("%d",&n);for(int i=0;i<n;i++) scanf("%d",&f[i]);poly_ln(n);for(int i=0;i<n;i++) printf("%d ",d[i]);return 0; }【奇奇怪怪的東西三】多項式指數(shù)函數(shù)
對于給定 f(x) 求h(x)= e^(f(x)) mod(x^n)
有點鬼畜是不是= =||?
還是先講一些前綴姿勢
1. 泰勒展開
對于一個函數(shù) f(x) 我們可以使用高階導數(shù)對其無限逼近
2.牛頓迭代
我們現(xiàn)在要求?其中g(shù)已知
假設我們知道原式的答案為
根據(jù)泰勒展開
?
顯然所以從第三項開始都是模x^n 為0的 我們可以不用考慮
那么就是
拆開移項得
那么我們就可以遞歸求解啦= =+
?
誒你突然發(fā)現(xiàn)問題了對不對= =+
g我們好像還不知道是什么呢
g當然是我們自己構(gòu)造的啦 由于兩邊同時取ln
得到
根據(jù)牛頓迭代的柿子??
所以把式子帶進去
這次就沒問題啦= =+
【喝了口開水差點被燙死】
#include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define inf 20021225 #define ll long long #define mxn 400100 #define mdn 998244353 #define G 3 using namespace std;int ksm(int bs,int mi) {int ans=1;while(mi){if(mi&1) ans=(ll)ans*bs%mdn;bs=(ll)bs*bs%mdn;mi>>=1;}return ans; } int rev[mxn],inv; int init(int n) {int lim=1,l=0;while(lim<n) lim<<=1,l++;for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);inv = ksm(lim,mdn-2);return lim; } void NTT(int *a,int lim,int f) {for(int i=0;i<lim;i++) if(rev[i]>i) swap(a[rev[i]],a[i]);for(int k=2;k<=lim;k<<=1){int mid=k>>1,Wn=ksm(G,(mdn-1)/k);if(f) Wn=ksm(Wn,mdn-2);for(int w=1,i=0;i<lim;i+=k,w=1)for(int j=0;j<mid;j++,w=(ll)w*Wn%mdn){int x=a[i+j],y=(ll)w*a[i+mid+j]%mdn;a[i+j]=(x+y)%mdn;a[i+mid+j]=(x-y+mdn)%mdn;}}if(f) for(int i=0;i<lim;i++) a[i]=(ll)a[i]*inv%mdn; } int g[mxn],h[mxn]; void poly_inv(int *a,int n) {if(n==1){g[0]=ksm(a[0],mdn-2);return;}poly_inv(a,n+1>>1);int lim=init(n<<1);for(int i=0;i<n;i++) h[i]=a[i];for(int i=n;i<lim;i++) h[i]=0;NTT(h,lim,0);NTT(g,lim,0);for(int i=0;i<lim;i++)g[i]=(2ll-(ll)g[i]*h[i]%mdn+mdn)%mdn*g[i]%mdn;NTT(g,lim,1);for(int i=n;i<lim;i++) g[i]=0; } int d[mxn],e[mxn]; int get_inv(int x){return ksm(x,mdn-2);} void get_d(int *a,int n) {d[n-1]=0;for(int i=0;i<n;i++) d[i]=(ll)a[i+1]*(i+1)%mdn; } void get_e(int *a,int n) {e[0]=0;for(int i=1;i<n;i++) e[i]=(ll)a[i-1]*get_inv(i)%mdn; } void poly_ln(int *a,int n) {get_d(a,n);poly_inv(a,n);int lim=init(n<<1);NTT(d,lim,0);NTT(g,lim,0);for(int i=0;i<lim;i++) g[i]=(ll)d[i]*g[i]%mdn;NTT(g,lim,1);get_e(g,n); } int s[mxn],tmp[mxn]; void poly_exp(int *a,int n) {if(n==1){s[0]=1;return;}int mid=n+1>>1;poly_exp(a,mid);for(int i=0;i<(n<<1);i++) e[i]=g[i]=0;poly_ln(s,n);for(int i=0;i<n;i++) tmp[i]=a[i];int lim=init(n<<1);NTT(e,lim,0);NTT(tmp,lim,0);NTT(s,lim,0);for(int i=0;i<lim;i++)s[i]=((ll)tmp[i]-e[i]+1+mdn)%mdn*s[i]%mdn;NTT(s,lim,1);for(int i=n;i<lim;i++) s[i]=tmp[i]=0; } int f[mxn]; int main() {int n;scanf("%d",&n);for(int i=0;i<n;i++) scanf("%d",&f[i]);poly_exp(f,n);for(int i=0;i<n;i++) printf("%d ",s[i]);return 0; }【奇奇怪怪的東西五】多項式除法
聽起來是不是十分酷炫 實際上用到一個神奇的轉(zhuǎn)化就可以輕輕松松松松松(大霧)的通過啦
我們對于給定 n次多項式f(x) 和 m次多項式g(x) 求中的 d 和 r 沒錯就是平常見過的多項式除法 其中(n>m)
這個玩意如何實現(xiàn)呢? 直接做的話肯定是O(nm) 十分不優(yōu)秀 況且既然有了 FFT這種nb的東西 怎么能讓這么不優(yōu)美的復雜度存在呢
下面來介紹這個神奇的操作
我們將一個多項式的系數(shù)翻轉(zhuǎn)【沒錯 就是 前后倒過來那個翻轉(zhuǎn)】
這個應該怎么在數(shù)學中實現(xiàn)呢 就是這個樣子
然后我們來干一些有趣的事情 把前面的所有多項式中的x替換成然后等式兩邊同時乘?于是就有了
然后我們發(fā)現(xiàn)這個玩意很優(yōu)美 可以化成
我們發(fā)現(xiàn) 余數(shù)項的最低次都是n-m+1 所以我們可以讓整個柿子對取模來消除r對答案的影響。
然后我們就可以魚塊的多項式求逆來求出再翻轉(zhuǎn)回來得到 帶回原式把求出來就做完啦
所以其實比前面的還要好寫= =+
代碼。
#include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define inf 20021225 #define ll long long #define mdn 998244353 #define mxn 400100 #define G 3 using namespace std;int ksm(int bs,int mi) {int ans=1;while(mi){if(mi&1) ans=(ll)ans*bs%mdn;bs=(ll)bs*bs%mdn;mi>>=1;}return ans; } int inv,rev[mxn]; void NTT(int *a,int lim,int f) {for(int i=1;i<lim;i++) if(rev[i]>i) swap(a[rev[i]],a[i]);for(int k=2;k<=lim;k<<=1){int mid=k>>1,Wn=ksm(G,(mdn-1)/k);if(f) Wn=ksm(Wn,mdn-2);for(int i=0,w=1;i<lim;w=1,i+=k)for(int j=0;j<mid;j++,w=(ll)w*Wn%mdn){int x=a[i+j],y=(ll)w*a[i+mid+j]%mdn;a[i+j]=(x+y)%mdn;a[i+mid+j]=(x-y+mdn)%mdn;}}if(f) for(int i=0;i<lim;i++) a[i]=(ll)a[i]*inv%mdn; } int g[mxn],h[mxn]; int init(int n) {int lim=1,l=0;while(lim<n) lim<<=1,l++;for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);inv=ksm(lim,mdn-2);return lim; } void poly_inv(int *a,int n)// g {if(n==1){g[0]=ksm(a[0],mdn-2);return;}int mid=(n+1)>>1;poly_inv(a,mid);int lim=init(n<<1);for(int i=0;i<n;i++) h[i]=a[i];for(int i=n;i<lim;i++) h[i]=0;NTT(h,lim,0);NTT(g,lim,0);for(int i=0;i<lim;i++)g[i]=(2ll -(ll)h[i]*g[i]%mdn +mdn)%mdn*g[i]%mdn;NTT(g,lim,1);for(int i=n;i<lim;i++) g[i]=0; } int ff[mxn],f[mxn],d[mxn],fd[mxn],r[mxn]; void reverse(int *a,int *b,int n) {for(int i=0;i<n;i++)b[i]=a[n-i-1]; } void poly_div(int n,int m) {reverse(f,ff,n);reverse(d,fd,m);int nn=n-m+1;poly_inv(fd,nn);//for(int i=0;i<nn/2;i++) swap(g[i],g[nn-i-1]);int lim=init(n<<1);// 2n-m+1NTT(g,lim,0);NTT(ff,lim,0);for(int i=0;i<lim;i++) g[i]=(ll)g[i]*ff[i]%mdn;NTT(g,lim,1);for(int i=nn;i<lim;i++) g[i]=0;for(int i=0;i<nn/2;i++) swap(g[i],g[nn-i-1]);NTT(g,lim,0);NTT(d,lim,0);NTT(f,lim,0);for(int i=0;i<lim;i++) d[i]=(ll)d[i]*g[i]%mdn;for(int i=0;i<lim;i++) r[i]=(f[i]-d[i]+mdn)%mdn;NTT(r,lim,1);NTT(g,lim,1);for(int i=0;i<nn;i++) printf("%d ",g[i]);printf("\n");for(int i=0;i<m-1;i++) printf("%d ",r[i]);printf("\n"); } int main() {int n,m;scanf("%d%d",&n,&m);n++;m++;for(int i=0;i<n;i++) scanf("%d",&f[i]);for(int i=0;i<m;i++) scanf("%d",&d[i]);poly_div(n,m);return 0; }持續(xù)更新!= =+
多項式全家桶大概到這里就告一段落啦~
或許什么時候我會寫多點求值,但那也是要先學完插值的啦。
所以到這里
完結(jié)撒花?(? ? ??)
我竟然沒有鴿
轉(zhuǎn)載于:https://www.cnblogs.com/hanyuweining/p/10321933.html
總結(jié)
以上是生活随笔為你收集整理的【学习笔记】之多项式使人头秃的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 浅拷贝和深拷贝的应用
- 下一篇: python的递归