多项式基础操作 - 学习笔记
生活随笔
收集整理的這篇文章主要介紹了
多项式基础操作 - 学习笔记
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
原文鏈接https://www.cnblogs.com/zhouzhendong/p/polynomial.html?
下載鏈接:
多項式基礎操作
#include <bits/stdc++.h> using namespace std; typedef long long LL; LL read(){LL x=0,f=0;char ch=getchar();while (!isdigit(ch))f|=ch=='-',ch=getchar();while (isdigit(ch))x=(x<<1)+(x<<3)+(ch^48),ch=getchar();return f?-x:x; } const int N=1<<18,mod=998244353; void Add(int &x,int y){if ((x+=y)>=mod)x-=mod; } void Del(int &x,int y){if ((x-=y)<0)x+=mod; } int del(int x,int y){return x-y<0?x-y+mod:x-y; } int Pow(int x,int y){int ans=1;for (;y;y>>=1,x=(LL)x*x%mod)if (y&1)ans=(LL)ans*x%mod;return ans; } int randint(){return ((rand()&65535)<<15)^(rand()&65535); } namespace Rem2{int INIT_TAG=0;int t,w;#define fi first#define se secondvoid init(){INIT_TAG=1;srand('C'+'L'+'Y'+'A'+'K'+'I'+'O'+'I');}pair <int,int> Mul_pii(pair <int,int> A,pair <int,int> B){static int a,b;a=((LL)A.fi*B.fi+(LL)A.se*B.se%mod*w)%mod;b=((LL)A.fi*B.se+(LL)A.se*B.fi)%mod;return make_pair(a,b);}pair <int,int> Pow_pii(pair <int,int> x,int y){pair <int,int> ans=make_pair(1,0);for (;y;y>>=1,x=Mul_pii(x,x))if (y&1)ans=Mul_pii(ans,x);return ans;}int Sqrt(int x){if (!INIT_TAG)init();if (x==0)return 0;if (Pow(x,(mod-1)/2)!=1)return -1;do {t=randint()%(mod-1)+1;w=((LL)t*t+mod-x)%mod;} while (Pow(w,(mod-1)/2)==1);pair <int,int> res=Pow_pii(make_pair(t,1),(mod+1)/2);return min(res.fi,mod-res.fi);} } namespace Polynomial{namespace Fast{const int N=1<<18;int n,Log[N+1],Fac[N+1],InvFac[N+1],Inv[N+1];int ww[N*2],*Ew=ww,*w[N+1];int iww[N*2],*Ei=iww,*iw[N+1];int INIT_TAG=0;void init(int _n){INIT_TAG=1;Log[1]=0,n=_n;for (int i=2;i<=N;i++)Log[i]=Log[i>>1]+1;for (int i=Fac[0]=1;i<=N;i++)Fac[i]=(LL)Fac[i-1]*i%mod;InvFac[N]=Pow(Fac[N],mod-2);for (int i=N;i>=1;i--)InvFac[i-1]=(LL)InvFac[i]*i%mod;for (int i=1;i<=N;i++)Inv[i]=(LL)InvFac[i]*Fac[i-1]%mod;for (int d=0;d<=Log[n];d++){w[d]=Ew,iw[d]=Ei;int n=1<<d;w[d][0]=1,w[d][1]=Pow(3,(mod-1)/n);for (int i=2;i<n;i++)w[d][i]=(LL)w[d][i-1]*w[d][1]%mod;iw[d][0]=1,iw[d][1]=Pow(w[d][1],mod-2);for (int i=2;i<n;i++)iw[d][i]=(LL)iw[d][i-1]*iw[d][1]%mod;Ew+=n,Ei+=n;}}int Rev[N+1],A[N+1],B[N+1];void FFT(int a[],int n,int **w){if (!INIT_TAG)init(N);for (int i=0;i<n;i++)if (Rev[i]<i)swap(a[i],a[Rev[i]]);for (int t=1,d=1;d<n;t++,d<<=1)for (int i=0;i<n;i+=(d<<1))for (int j=0,*W=w[t];j<d;j++){int tmp=(LL)(*W++)*a[i+j+d]%mod;a[i+j+d]=del(a[i+j],tmp);Add(a[i+j],tmp);}}vector <int> Mul(vector <int> &a,vector <int> &b){static vector <int> res;res.clear();LL Br=(LL)a.size()*b.size();LL FF=(a.size()+b.size())*Log[a.size()+b.size()]*10+100;if (Br<=FF){for (int i=0;i<a.size()+b.size();i++)res.push_back(0);for (int i=0;i<a.size();i++)for (int j=0;j<b.size();j++)res[i+j]=((LL)a[i]*b[j]+res[i+j])%mod;}else {int n=1,d=0;for (;n<a.size()+b.size();n<<=1,d++);for (int i=0;i<n;i++)Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(d-1)),A[i]=B[i]=0;for (int i=0;i<a.size();i++)A[i]=a[i];for (int i=0;i<b.size();i++)B[i]=b[i]; // w[0]=1,w[1]=Pow(3,(mod-1)/n); // for (int i=2;i<n;i++) // w[i]=(LL)w[i-1]*w[1]%mod; FFT(A,n,w),FFT(B,n,w);for (int i=0;i<n;i++)A[i]=(LL)A[i]*B[i]%mod; // w[1]=Pow(w[1],mod-2); // for (int i=2;i<n;i++) // w[i]=(LL)w[i-1]*w[1]%mod; FFT(A,n,iw);int inv=Pow(n,mod-2);for (int i=0;i<n;i++)res.push_back((int)((LL)inv*A[i]%mod));}while (!res.empty()&&!res.back())res.pop_back();return res;}vector <int> MulInv(vector <int> &a,vector <int> &b){static vector <int> res;res.clear();int n=1,d=0;for (;n<a.size()*2+b.size();n<<=1,d++);for (int i=0;i<n;i++)Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(d-1)),A[i]=B[i]=0;for (int i=0;i<a.size();i++)A[i]=a[i];for (int i=0;i<b.size();i++)B[i]=b[i]; // w[0]=1,w[1]=Pow(3,(mod-1)/n); // for (int i=2;i<n;i++) // w[i]=(LL)w[i-1]*w[1]%mod; FFT(A,n,w),FFT(B,n,w);for (int i=0;i<n;i++)A[i]=(LL)A[i]*A[i]%mod*B[i]%mod; // w[1]=Pow(w[1],mod-2); // for (int i=2;i<n;i++) // w[i]=(LL)w[i-1]*w[1]%mod; FFT(A,n,iw);int inv=Pow(n,mod-2);for (int i=0;i<n;i++)res.push_back((int)((LL)inv*A[i]%mod));while (!res.empty()&&!res.back())res.pop_back();return res;}}struct Poly{vector <int> v;Poly(){v.clear();}Poly(int x){v.clear();v.push_back(x);}Poly(vector <int> x){v=x;}int operator ()(int x){int ans=0,y=1;for (int i=0;i<v.size();i++)ans=((LL)v[i]*y+ans)%mod,y=(LL)y*x%mod;return ans;}int size(){return v.size();}void print(){for (int i=0;i<v.size();i++)printf("%d ",v[i]);}void print(int x){for (int i=0;i<x;i++)printf("%d ",i>=v.size()?0:v[i]);}void print(string s){print(),cout << s;}void clear(){v.clear();}void push_back(int x){v.push_back(x);}void pop_back(){v.pop_back();}int empty(){return v.empty();}int back(){return v.back();}int &operator [](int x){return v[x];}void operator += (Poly A){while (v.size()<A.size())v.push_back(0);for (int i=0;i<A.size();i++)Add(v[i],A[i]);}void operator -= (Poly &A){while (v.size()<A.size())v.push_back(0);for (int i=0;i<A.size();i++)Del(v[i],A[i]);}void operator *= (Poly &A);void Derivation(){for (int i=0;i<v.size()-1;i++)v[i]=(LL)v[i+1]*(i+1)%mod;v.pop_back();}void Integral(){v.push_back(0);for (int i=v.size()-2;i>=0;i--)v[i+1]=(LL)v[i]*Fast :: Inv[i+1]%mod;v[0]=0;}void operator *= (int x){for (int i=0;i<v.size();i++)v[i]=(LL)v[i]*x%mod;}}pp;//struct Poly end-------------Poly operator + (Poly A,Poly B){pp.clear();for (int i=0;i<max(A.size(),B.size());i++)pp.push_back(0);for (int i=0;i<A.size();i++)Add(pp[i],A[i]);for (int i=0;i<B.size();i++)Add(pp[i],B[i]);return pp;}Poly operator - (Poly A,Poly B){pp.clear();for (int i=0;i<max(A.size(),B.size());i++)pp.push_back(0);for (int i=0;i<A.size();i++)Add(pp[i],A[i]);for (int i=0;i<B.size();i++)Del(pp[i],B[i]);return pp;}Poly operator * (Poly A,Poly B){return Poly(Fast :: Mul(A.v,B.v));}void Poly :: operator *= (Poly &A){v=Fast :: Mul(v,A.v);}Poly operator * (Poly A,int x){pp=A;for (int i=0;i<A.size();i++)pp[i]=(LL)pp[i]*x%mod;return pp;}Poly Inverse(Poly a,int n);Poly operator / (Poly A,Poly B){//Divideint n=A.size(),m=B.size();reverse(A.v.begin(),A.v.end());reverse(B.v.begin(),B.v.end());int k=n-m+1;if (k<0)return Poly(0);while (A.size()>k)A.pop_back();while (B.size()>k)B.pop_back();A=A*Inverse(B,k);while (A.size()>k)A.pop_back();reverse(A.v.begin(),A.v.end());return A;}Poly operator % (Poly A,Poly B){//Modulowhile (!A.empty()&&!A.back())A.pop_back();while (!B.empty()&&!B.back())B.pop_back();A=A-A/B*B;while (A.size()>=B.size())A.pop_back();while (!A.empty()&&!A.back())A.pop_back();return A;}Poly Derivation(Poly A){for (int i=0;i<A.size()-1;i++)A[i]=(LL)A[i+1]*(i+1)%mod;A.pop_back();return A;}Poly Integral(Poly A){A.push_back(0);for (int i=A.size()-2;i>=0;i--)A[i+1]=(LL)A[i]*Fast :: Inv[i+1]%mod;A[0]=0;return A;}Poly Inverse(Poly a,int n){static Poly A,B;while (!a.empty()&&!a.back())a.pop_back();if (a.empty())return a;A.clear(),B.clear();B.push_back(a[0]);A.push_back(Pow(B[0],mod-2));for (int t=1;t<n;){for (int i=t;i<min(a.size(),(t<<1));i++)B.push_back(a[i]);t<<=1;A=A*2-Poly(Fast :: MulInv(A.v,B.v));while (A.size()>t)A.pop_back();}while (A.size()>n)A.pop_back();return A;}Poly Sqrt(Poly a,int n){static Poly A,B;while (!a.empty()&&!a.back())a.pop_back();if (a.empty())return a;A.clear(),B.clear();B.push_back(a[0]);A.push_back(Rem2 :: Sqrt(B[0]));for (int t=1;t<n;){for (int i=t;i<min(a.size(),(t<<1));i++)B.push_back(a[i]);t<<=1;A+=B*Inverse(A,t);while (A.size()>t)A.pop_back();A*=499122177;}if (A[0]>mod-A[0])for (int i=0;i<A.size();i++)A[i]=(mod-A[i])%mod;while (A.size()>n)A.pop_back();return A;}Poly Ln(Poly a,int n){while (!a.empty()&&!a.back())a.pop_back();if (a.empty()||a[0]!=1)return a;a=Integral(Derivation(a)*Inverse(a,n));while (a.size()>n)a.pop_back();return a;}Poly Exp(Poly a,int n){static Poly A,B;while (!a.empty()&&!a.back())a.pop_back();if (a.empty())return Poly(1);if (a[0]!=0)return a;A.clear(),B.clear();B.push_back(1);A.push_back(a[0]);for (int t=1;t<n;){for (int i=t;i<min(a.size(),(t<<1));i++)A.push_back(a[i]);t<<=1;B=B*(Poly(1)+A-Ln(B,t));while (B.size()>t)B.pop_back();}while (B.size()>n)B.pop_back();return B;}Poly PolyPow(Poly x,int y,int n){static Poly A,B;int k0=0,kc,ivkc;while (!x.empty()&&!x.back())x.pop_back();if (x.empty())return x;while (k0<x.size()&&x[k0]==0)k0++;kc=x[k0],ivkc=Pow(kc,mod-2);A.clear();for (int i=k0;i<x.size();i++)A.push_back((int)((LL)x[i]*ivkc%mod));A=Exp(Ln(A,n)*y,n);B.clear();if ((LL)k0*y>=n)return B;kc=Pow(kc,y),k0*=y;for (int i=0;i<k0;i++)B.push_back(0);for (int i=0;i<min(A.size(),n-k0);i++)B.push_back((int)((LL)A[i]*kc%mod));while (B.size()>n)B.pop_back();return B;}namespace Qiuzhi{Poly P[N<<2],f[N<<2],M;vector <int> x,y;int n;void GetP(int rt,int L,int R){if (L==R){P[rt].clear();P[rt].push_back((mod-x[L])%mod);P[rt].push_back(1);return;}int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;GetP(ls,L,mid);GetP(rs,mid+1,R);P[rt]=P[ls]*P[rs];}void qiuzhi(int rt,int L,int R){if (f[rt].empty())f[rt].push_back(0);if (L==R)return (void)(y[L]=f[rt][0]);int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;f[ls]=f[rt]%P[ls];f[rs]=f[rt]%P[rs];qiuzhi(ls,L,mid);qiuzhi(rs,mid+1,R);}vector <int> Get_Val(vector <int> A,Poly F){n=A.size();x.clear(),y.clear();for (int i=0;i<n;i++){x.push_back(A[i]);y.push_back(0);}GetP(1,0,n-1);f[1]=F;qiuzhi(1,0,n-1);return y;}}namespace Chazhi{Poly P[N<<2],M;vector <int> x,y;int n;void GetP(int rt,int L,int R){if (L==R){P[rt].clear();P[rt].push_back((mod-x[L])%mod);P[rt].push_back(1);return;}int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;GetP(ls,L,mid);GetP(rs,mid+1,R);P[rt]=P[ls]*P[rs];}Poly chazhi(int rt,int L,int R){if (L==R)return Poly(y[L]);int mid=(L+R)>>1,ls=rt<<1,rs=ls|1;return chazhi(ls,L,mid)*P[rs]+chazhi(rs,mid+1,R)*P[ls];}Poly Get_Poly(vector <int> A,vector <int> B){n=A.size();x=A;int Product=1;GetP(1,0,n-1);M=Derivation(P[1]);y=Qiuzhi :: Get_Val(A,M);for (int i=0;i<y.size();i++)y[i]=(LL)B[i]*Pow(y[i],mod-2)%mod;return chazhi(1,0,n-1);}} }// be careful about init!!!!!! using namespace Polynomial; Poly A,B; vector <int> x,y; int main(){int n=read()+1,k=read();A.clear();for (int i=0;i<n;i++)A.push_back(read());B=Exp(Integral(Inverse(Sqrt(A,n),n)),n);B=Poly(1)+Ln(Poly(2)+A-A(0)-B,n);B=Derivation(PolyPow(B,k,n));n--;B.print(n);return 0; } Polynomial?
UPD(2019-03-12): 感謝神仙cly的指出。已更。
UPD(2019-04-21): 補一個相對好寫一些的LOJ#150板子。
#include <bits/stdc++.h> #define clr(x) memset(x,0,sizeof (x)) #define clrint(x,n) memset(x,0,(n)<<2) #define cpyint(a,b,n) memcpy(a,b,(n)<<2) #define For(i,a,b) for (int i=a;i<=b;i++) #define Fod(i,b,a) for (int i=b;i>=a;i--) #define pb(x) push_back(x) #define mp(x,y) make_pair(x,y) #define fi first #define se second #define real __zzd001 #define _SEED_ ('C'+'L'+'Y'+'A'+'K'+'I'+'O'+'I') #define outval(x) printf(#x" = %d\n",x) #define outvec(x) printf("vec "#x" = ");for (auto _v : x)printf("%d ",_v);puts("") #define outtag(x) puts("----------"#x"----------") #define outarr(a,L,R) printf(#a"[%d...%d] = ",L,R);\For(_v2,L,R)printf("%d ",a[_v2]);puts(""); using namespace std; typedef long long LL; typedef unsigned long long ULL; typedef vector <int> vi; LL read(){LL x=0,f=0;char ch=getchar();while (!isdigit(ch))f|=ch=='-',ch=getchar();while (isdigit(ch))x=(x<<1)+(x<<3)+(ch^48),ch=getchar();return f?-x:x; } const int N=1<<19,mod=998244353,inv2=(mod+1)>>1; const int YG=3; int Pow(int x,int y){int ans=1;for (;y;y>>=1,x=(LL)x*x%mod)if (y&1)ans=(LL)ans*x%mod;return ans; } void Add(int &x,int y){if ((x+=y)>=mod)x-=mod; } void Del(int &x,int y){if ((x-=y)<0)x+=mod; } int Add(int x){return x>=mod?x-mod:x; } int Del(int x){return x<0?x+mod:x; } namespace Math{int Iv[N];void prework(){int n=N-1;Iv[1]=1;For(i,2,n)Iv[i]=(LL)(mod-mod/i)*Iv[mod%i]%mod;}map <int,int> Map;int ind(int x){static int M,bas;if (Map.empty()){M=max((int)sqrt(mod),1);bas=Pow(YG,M);for (int i=1,v=YG;i<=M;i++,v=(LL)v*YG%mod)Map[v]=i;}for (int i=M,v=(LL)bas*Pow(x,mod-2)%mod;i<=mod-1+M;i+=M,v=(LL)v*bas%mod)if (Map[v])return i-Map[v];return -1;} } namespace fft{int w[N],R[N];int Log[N+1];void init(int n){if (!Log[2]){For(i,2,N)Log[i]=Log[i>>1]+1;}int d=Log[n];assert(n==(1<<d));For(i,0,n-1)R[i]=(R[i>>1]>>1)|((i&1)<<(d-1));w[0]=1,w[1]=Pow(YG,(mod-1)/n);For(i,2,n-1)w[i]=(LL)w[i-1]*w[1]%mod;}void FFT(int *a,int n,int flag){if (flag<0)reverse(w+1,w+n);For(i,0,n-1)if (i<R[i])swap(a[i],a[R[i]]);for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)for (int i=0;i<n;i+=d<<1)for (int j=0;j<d;j++){int tmp=(LL)w[t*j]*a[i+j+d]%mod;a[i+j+d]=Del(a[i+j]-tmp);Add(a[i+j],tmp);}if (flag<0){reverse(w+1,w+n);int inv=Pow(n,mod-2);For(i,0,n-1)a[i]=(LL)a[i]*inv%mod;}}void CirMul(int *a,int *b,int *c,int n){init(n),FFT(a,n,1),FFT(b,n,1);For(i,0,n-1)c[i]=(LL)a[i]*b[i]%mod;FFT(c,n,-1);} } using fft::FFT; using fft::CirMul; int calc_up(int x){int n=1;while (n<=x)n<<=1;return n; } void Inv(int *a,int *b,int n){static int f[N],g[N];b[0]=Pow(a[0],mod-2);int now=1;while (now<n){int len=now<<2;For(i,0,len-1)f[i]=g[i]=0;cpyint(g,b,now),now<<=1,cpyint(f,a,min(n,now));fft::init(len);FFT(f,len,1),FFT(g,len,1);For(i,0,len-1)g[i]=(2LL*g[i]-(LL)f[i]*g[i]%mod*g[i]%mod+mod)%mod;FFT(g,len,-1);cpyint(b,g,min(n,now));} } int Sqrt(int a){int k=Math::ind(a);assert(~k&1);k=Pow(YG,k>>1);return min(k,mod-k); } void Sqrt(int *a,int *b,int n){static int f[N],g[N],h[N];b[0]=Sqrt(a[0]);int now=1;while (now<n){int len=now<<2;For(i,0,len-1)f[i]=g[i]=h[i]=0;cpyint(f,b,now),now<<=1,Inv(f,h,now),cpyint(g,a,min(n,now));CirMul(g,h,g,len);For(i,0,len-1)f[i]=((g[i]+f[i])&1)?Add(((LL)g[i]+f[i]+mod)>>1):((g[i]+f[i])>>1);cpyint(b,f,min(n,now));} } void Der(int *a,int n){For(i,0,n-2)a[i]=(LL)a[i+1]*(i+1)%mod;a[n-1]=0; } void Int(int *a,int n){if (!Math::Iv[1])Math::prework();Fod(i,n,1)a[i]=(LL)a[i-1]*Math::Iv[i]%mod;a[0]=0; } void Ln(int *a,int *b,int n){static int f[N],g[N];int len=calc_up(n*2);For(i,0,len-1)f[i]=g[i]=0;cpyint(f,a,n),Inv(f,g,n),Der(f,n);CirMul(f,g,f,len);Int(f,n),cpyint(b,f,n); } void Exp(int *a,int *b,int n){static int f[N],g[N],h[N];b[0]=1;int now=1;while (now<n){int len=now<<2;For(i,0,len-1)f[i]=g[i]=h[i]=0;cpyint(f,b,now),now<<=1,Ln(f,g,now),cpyint(h,a,min(n,now));For(i,0,now-1)g[i]=Del(h[i]-g[i]);Add(g[0],1);CirMul(f,g,f,len),cpyint(b,f,min(n,now));} } void Pow(int *a,int *b,int n,int k){static int f[N];clrint(b,n);if (k==0)return (void)(b[0]=1);int fir=0;for (;fir<n&&!a[fir];fir++);if ((LL)fir*k>=n)return;int m=n-fir*k;cpyint(f,a+fir,m);int t=Pow(f[0],k),it=Pow(f[0],mod-2);For(i,0,m-1)f[i]=(LL)f[i]*it%mod;Ln(f,f,m);For(i,0,m-1)f[i]=(LL)f[i]*k%mod;Exp(f,b+fir*k,m);For(i,fir*k,n-1)b[i]=(LL)b[i]*t%mod; } int n,k; int f[N],g[N],h[N]; int main(){n=read()+1,k=read();For(i,0,n-1)f[i]=read();Sqrt(f,g,n);Inv(g,h,n);Int(h,n);Exp(h,g,n);For(i,0,n-1)g[i]=Del(f[i]-g[i]);Add(g[0],2),Del(g[0],f[0]);Ln(g,h,n);Add(h[0],1);Pow(h,g,n,k);Der(g,n);For(i,0,n-2)printf("%d ",g[i]);puts("");return 0; } View CodeUPD(2019-04-21晚): 再加一個歷時56min Rush出來的板子。感覺看起來比上面那個更好看。相比之下唯一的缺點就是要手動調用一個 prework() 函數。
#include <bits/stdc++.h> #define clr(x) memset(x,0,sizeof x) #define cpyint(a,b,n) memcpy(a,b,(n)<<2) #define For(i,a,b) for (int i=a;i<=b;i++) #define Fod(i,b,a) for (int i=b;i>=a;i--) #define outval(x) printf(#x" = %d\n",x) #define outtag(x) puts("-------------"#x"-------------"); #define outarr(a,L,R) printf(#a"[%d..%d] = ",L,R);\For(_x,L,R)printf("%d ",a[_x]);puts("") using namespace std; typedef long long LL; LL read(){LL x=0,f=0;char ch=getchar();while (!isdigit(ch))f|=ch=='-',ch=getchar();while (isdigit(ch))x=(x<<1)+(x<<3)+(ch^48),ch=getchar();return f?-x:x; } const int N=1<<19,mod=998244353; int Pow(int x,int y){int ans=1;for (;y;y>>=1,x=(LL)x*x%mod)if (y&1)ans=(LL)ans*x%mod;return ans; } void Add(int &x,int y){if ((x+=y)>=mod)x-=mod; } void Del(int &x,int y){if ((x-=y)<0)x+=mod; } int Add(int x){return x>=mod?x-mod:x; } int Del(int x){return x<0?x+mod:x; } namespace fft{int w[N],R[N];void init(int n){int d=0;while ((1<<d)<n)d++;For(i,0,n-1)R[i]=(R[i>>1]>>1)|((i&1)<<(d-1));w[0]=1,w[1]=Pow(3,(mod-1)/n);For(i,2,n-1)w[i]=(LL)w[i-1]*w[1]%mod;}void FFT(int *a,int n,int flag){if (flag<0)reverse(w+1,w+n);For(i,0,n-1)if (i<R[i])swap(a[i],a[R[i]]);for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)for (int i=0;i<n;i+=d<<1)for (int j=0;j<d;j++){int tmp=(LL)w[t*j]*a[i+j+d]%mod;a[i+j+d]=Del(a[i+j]-tmp);Add(a[i+j],tmp);}if (flag<0){reverse(w+1,w+n);int inv=Pow(n,mod-2);For(i,0,n-1)a[i]=(LL)a[i]*inv%mod;}}void Mul(int *a,int *b,int *c,int n){init(n),FFT(a,n,1),FFT(b,n,1);For(i,0,n-1)c[i]=(LL)a[i]*b[i]%mod;FFT(c,n,-1);} } using fft::FFT; using fft::Mul; int Iv[N]; void prework(){int n=N-1;Iv[1]=1;For(i,2,n)Iv[i]=(LL)(mod-mod/i)*Iv[mod%i]%mod; } int Ind(int x){static map <int,int> Map;static int M,bas;if (Map.empty()){M=(int)sqrt(mod);bas=Pow(3,M);for (int i=1,v=3;i<=M;i++,v=(LL)v*3%mod)Map[v]=i;}for (int i=M,v=(LL)bas*Pow(x,mod-2)%mod;i<=mod+M;i+=M,v=(LL)v*bas%mod)if (Map[v])return i-Map[v];return -1; } int calup(int x){int n=1;while (n<x)n<<=1;return n; } void Inv(int *a,int *b,int n){static int f[N],g[N];b[0]=Pow(a[0],mod-2);int now=1;while (now<n){int len=now<<2;For(i,0,len-1)f[i]=g[i]=0;cpyint(g,b,now),now<<=1,cpyint(f,a,min(n,now));fft::init(len),FFT(f,len,1),FFT(g,len,1);For(i,0,len-1)f[i]=(2LL*g[i]-(LL)f[i]*g[i]%mod*g[i]%mod+mod)%mod;FFT(f,len,-1);cpyint(b,f,min(n,now));} } int Sqrt(int x){int k=Ind(x);assert(~k&1);k=Pow(3,k>>1);return min(k,mod-k); } void Sqrt(int *a,int *b,int n){static int f[N],g[N],h[N];b[0]=Sqrt(a[0]);int now=1;while (now<n){int len=now<<2;For(i,0,len-1)f[i]=g[i]=h[i]=0;cpyint(f,b,now),now<<=1,Inv(f,g,now),cpyint(h,a,min(n,now));Mul(g,h,g,len);For(i,0,len-1)f[i]=Add(f[i]+g[i]),f[i]=(f[i]&1)?((f[i]+mod)>>1):(f[i]>>1);cpyint(b,f,min(n,now));} } void Der(int *a,int n){For(i,0,n-2)a[i]=(LL)a[i+1]*(i+1)%mod;a[n-1]=0; } void Int(int *a,int n){Fod(i,n,1)a[i]=(LL)a[i-1]*Iv[i]%mod;a[0]=0; } void Ln(int *a,int *b,int n){static int f[N],g[N];int len=calup(n<<1);For(i,0,len-1)f[i]=g[i]=0;cpyint(f,a,n),Inv(f,g,n),Der(f,n);Mul(f,g,f,len);Int(f,n),cpyint(b,f,n); } void Exp(int *a,int *b,int n){static int f[N],g[N],h[N];b[0]=1;int now=1;while (now<n){int len=now<<2;For(i,0,len-1)f[i]=g[i]=h[i]=0;cpyint(g,b,now),now<<=1,Ln(g,h,now),cpyint(f,a,min(n,now));For(i,0,now-1)Del(f[i],h[i]);Add(f[0],1);Mul(g,f,g,len);cpyint(b,g,min(n,now));} } void Pow(int *a,int *b,int n,int k){static int f[N],g[N];memset(b,0,sizeof(int)*n);if (k==0)return (void)(b[0]=1);int fir=0;for (;fir<n&&!a[fir];fir++);if ((LL)fir*k>=n)return;int m=n-fir*k;cpyint(f,a+fir,m);int t=Pow(f[0],k),it=Pow(f[0],mod-2);For(i,0,m-1)f[i]=(LL)f[i]*it%mod;Ln(f,g,m);For(i,0,m-1)g[i]=(LL)g[i]*k%mod;Exp(g,f,m);For(i,0,m-1)f[i]=(LL)f[i]*t%mod;cpyint(b+fir*k,f,m); } int n,k; int a[N],b[N],c[N],d[N]; int main(){prework();n=read()+1,k=read();For(i,0,n-1)a[i]=read();Sqrt(a,b,n);Inv(b,c,n);Int(c,n);Exp(c,b,n);For(i,0,n-1)b[i]=Del(a[i]-b[i]);Add(b[0],2),Del(b[0],a[0]);Ln(b,c,n);Add(c[0],1);Pow(c,b,n,k);Der(b,n);n--;For(i,0,n-1)printf("%d ",b[i]);puts("");return 0; } View Code?
轉載于:https://www.cnblogs.com/zhouzhendong/p/Polynomial.html
總結
以上是生活随笔為你收集整理的多项式基础操作 - 学习笔记的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: [转帖]2019 简易Web开发指南
- 下一篇: Cocoapods ----- pod