开头Orz yww,Orz xfz,Orz dtz,Orz lyy
一些定义
一个关于$x$的多项式$A(x)$可以表示为如下形式和:
$$A(x)=\sum\limits_{i=0}^{n-1}a_{i}x^{i}$$
其中$a_0,a_1,a_2......a_{n-1}$称为多项式$A(x)$的系数,如果$A(x)$最高次数的一个非零项系数为$a_k$,则称其次数为$k$。任意一个严格大于$k$的整数都可以作为这个多项式的次数界,即对于一个次数界为$n$的多项式,他的次数可以是$[0,n-1]$中的任意一个整数。
多项式加减法
略
多项式乘法
直接FFT or NTT
代码:
FFT:
1 #include<iostream>
2 #include<cstring>
3 #include<cstdio>
4 #include<cmath>
5 #define pw(n) (1<<n)
6 using namespace std; 7 const double pi=acos(-1); 8 struct complex{ 9 double a,b; 10 complex(double _a=0,double _b=0){ 11 a=_a; 12 b=_b; 13 } 14 friend complex operator +(complex x,complex y){return complex(x.a+y.a,x.b+y.b);} 15 friend complex operator -(complex x,complex y){return complex(x.a-y.a,x.b-y.b);} 16 friend complex operator *(complex x,complex y){return complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);} 17 friend complex operator *(complex x,double y){return complex(x.a*y,x.b*y);} 18 friend complex operator /(complex x,double y){return complex(x.a/y,x.b/y);} 19 }a[100001],b[100001]; 20 int n,m,bit,bitnum=0,rev[pw(20)]; 21 void getrev(int l){//Reverse
22 for(int i=0;i<pw(l);i++){ 23 rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1)); 24 } 25 } 26 void FFT(complex *s,int op){ 27 for(int i=0;i<bit;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 28 for(int i=1;i<bit;i<<=1){ 29 complex w(cos(pi/i),op*sin(pi/i)); 30 for(int p=i<<1,j=0;j<bit;j+=p){//Butterfly
31 complex wk(1,0); 32 for(int k=j;k<i+j;k++,wk=wk*w){ 33 complex x=s[k],y=wk*s[k+i]; 34 s[k]=x+y; 35 s[k+i]=x-y; 36 } 37 } 38 } 39 if(op==-1){ 40 for(int i=0;i<=bit;i++){ 41 s[i]=s[i]/(double)bit; 42 } 43 } 44 } 45 int main(){ 46 scanf("%d%d",&n,&m); 47 for(int i=0;i<=n;i++)scanf("%lf",&a[i].a); 48 for(int i=0;i<=m;i++)scanf("%lf",&b[i].a); 49 m+=n; 50 for(bit=1;bit<=m;bit<<=1)bitnum++; 51 getrev(bitnum); 52 FFT(a,1); 53 FFT(b,1); 54 for(int i=0;i<=bit;i++)a[i]=a[i]*b[i]; 55 FFT(a,-1); 56 for(int i=0;i<=m;i++)printf("%d ",(int)(a[i].a+0.5)); 57 return 0; 58 }
NTT:
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #define pw(n) (1<<n)
7 using namespace std; 8 const int N=262144,P=998244353,g=3;//或P=1004535809
9 int n,m,bit,bitnum=0,a[N+5],b[N+5],rev[N+5]; 10 void getrev(int l){ 11 for(int i=0;i<pw(l);i++){ 12 rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1)); 13 } 14 } 15 int fastpow(int a,int b){ 16 int ans=1; 17 for(;b;b>>=1,a=1LL*a*a%P){ 18 if(b&1)ans=1LL*ans*a%P; 19 } 20 return ans; 21 } 22 void NTT(int *s,int op){ 23 for(int i=0;i<bit;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 24 for(int i=1;i<bit;i<<=1){ 25 int w=fastpow(g,(P-1)/(i<<1)); 26 for(int p=i<<1,j=0;j<bit;j+=p){ 27 int wk=1; 28 for(int k=j;k<i+j;k++,wk=1LL*wk*w%P){ 29 int x=s[k],y=1LL*s[k+i]*wk%P; 30 s[k]=(x+y)%P; 31 s[k+i]=(x-y+P)%P; 32 } 33 } 34 } 35 if(op==-1){ 36 reverse(s+1,s+bit); 37 int inv=fastpow(bit,P-2); 38 for(int i=0;i<bit;i++)a[i]=1LL*a[i]*inv%P; 39 } 40 } 41 int main(){ 42 scanf("%d%d",&n,&m); 43 for(int i=0;i<=n;i++)scanf("%d",&a[i]); 44 for(int i=0;i<=m;i++)scanf("%d",&b[i]); 45 m+=n; 46 for(bit=1;bit<=m;bit<<=1)bitnum++; 47 getrev(bitnum); 48 NTT(a,1); 49 NTT(b,1); 50 for(int i=0;i<bit;i++)a[i]=1LL*a[i]*b[i]%P; 51 NTT(a,-1); 52 for(int i=0;i<=m;i++)printf("%d ",a[i]); 53 return 0; 54 }
时间复杂度:$O(nlogn)$
模板:洛谷P3803
Upd:补一个任意模数fft(mtt)
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<vector>
5 #include<cstdio>
6 #include<cmath>
7 #include<queue>
8 #include<stack>
9 #include<map>
10 #include<set>
11 using namespace std; 12 typedef long long ll; 13 const long double pi=acos(-1); 14 struct complex{ 15 long double a,b; 16 complex(long double _a=0,long double _b=0){ 17 a=_a; 18 b=_b; 19 } 20 friend complex operator +(complex x,complex y){return complex(x.a+y.a,x.b+y.b);} 21 friend complex operator -(complex x,complex y){return complex(x.a-y.a,x.b-y.b);} 22 friend complex operator *(complex x,complex y){return complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);} 23 friend complex operator *(complex x,long double y){return complex(x.a*y,x.b*y);} 24 friend complex operator /(complex x,long double y){return complex(x.a/y,x.b/y);} 25 }a[1000001],b[1000001],c[1000001],d[1000001]; 26 int n,m,p,bit=1,bitnum=0,A[1000001],B[1000001],rev[1000001]; 27 void fft(complex *s,int op){ 28 for(int i=0;i<bit;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 29 for(int i=1;i<bit;i<<=1){ 30 complex w(cos(pi/i),op*sin(pi/i)); 31 for(int pp=i<<1,j=0;j<bit;j+=pp){ 32 complex wk(1,0); 33 for(int k=j;k<i+j;k++,wk=wk*w){ 34 complex x=s[k],y=wk*s[k+i]; 35 s[k]=x+y; 36 s[k+i]=x-y; 37 } 38 } 39 } 40 if(op==-1){ 41 for(int i=0;i<bit;i++){ 42 s[i]=s[i]/(double)bit; 43 } 44 } 45 } 46 int main(){ 47 scanf("%d%d%d",&n,&m,&p); 48 for(int i=0;i<=n;i++)scanf("%d",&A[i]); 49 for(int i=0;i<=m;i++)scanf("%d",&B[i]); 50 n+=m+1; 51 for(;bit<=n;bit<<=1)++bitnum; 52 for(int i=0;i<bit;i++){ 53 rev[i]=(rev[i>>1]>>1|((i&1)<<(bitnum-1))); 54 } 55 for(int i=0;i<bit;i++){ 56 a[i].a=A[i]>>15; 57 b[i].a=A[i]&32767; 58 c[i].a=B[i]>>15; 59 d[i].a=B[i]&32767; 60 } 61 fft(a,1); 62 fft(b,1); 63 fft(c,1); 64 fft(d,1); 65 for(int i=0;i<bit;i++){ 66 complex s1=a[i],s2=b[i],s3=c[i],s4=d[i]; 67 a[i]=s1*s3; 68 b[i]=s1*s4+s2*s3; 69 c[i]=s2*s4; 70 } 71 fft(a,-1); 72 fft(b,-1); 73 fft(c,-1); 74 for(int i=0;i<n;i++){ 75 printf("%lld ",(ll)((ll)(a[i].a+0.5)%p*(1LL<<30)%p+(ll)(b[i].a+0.5)%p*(1LL<<15)%p+(ll)(c[i].a+0.5)%p)%p); 76 } 77 return 0; 78 }
模板:洛谷P4245
多项式求逆
如果存在多项式$B(x)$使得:
$$A(x)B(x)≡1(\mod x^n)$$
则$B(x)$称为$A(x)$在模$x^n$意义下的乘法逆元;
一个多项式存在乘法逆元的充要条件是他的常数项存在乘法逆元(不会证)
假设已经求出了$B'(x)$满足:
$$A(x)B'(x)≡1(\mod x^{\lceil \frac{n}{2}\rceil})$$
则
$$A(x)B'(x)-1≡0(\mod x^{\lceil \frac{n}{2}\rceil})$$
$$(A(x)B'(x)-1)^2≡0(\mod x^n)$$
$$A^2(x)B'^2(x)-2A(x)B'(x)+1≡0(\mod x^n)$$
$$2A(x)B'(x)-A^2(x)B'^2(x)≡1(\mod x^n)$$
将此式减去$A(x)B(x)≡1(mod x^n)$,可得
$$2A(x)B'(x)-A^2(x)B'^2(x)-A(x)B(x)≡0(\mod x^n)$$
约去$A(x)$,得
$$2B'(x)-A(x)B'^2(x)-B(x)≡0(\mod x^n)$$
$$2B'(x)-A(x)B'^2(x)≡B(x) (\mod x^n)$$
于是就可以每次算出$2B'(x)-A(x)B'^2(x)$然后赋值给$B(x)$进行下一次迭代,因此时间复杂度为$T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)$,看着像是两个$log$的,实际上只有一个,这里提供一个lyy的证明:
设$n=2^k$
则$T(n)=T(\frac{n}{2})+O(nlogn)=\sum\limits_{i=0}^{k}\frac{n}{2^i}logn=nlogn$
dtz:轻松证明
代码:
1 #include<algorithm> 2 #include<iostream> 3 #include<cstring> 4 #include<vector> 5 #include<cstdio> 6 #include<cmath> 7 #include<queue> 8 #include<stack> 9 #include<map> 10 #include<set> 11 using namespace std; 12 typedef long long ll; 13 const int p=998244353,g=3; 14 int n,m,a[500001],b[500001],rev[500001],tmp[500001]; 15 int fastpow(int x,int y){ 16 int ret=1; 17 for(;y;y>>=1,x=(ll)x*x%p){ 18 if(y&1)ret=(ll)ret*x%p; 19 } 20 return ret; 21 } 22 void ntt(int s[],int n,int op){ 23 for(int i=0;i<n;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 24 for(int i=1;i<n;i<<=1){ 25 int w=fastpow(g,(p-1)/(i<<1)); 26 for(int pp=i<<1,j=0;j<n;j+=pp){ 27 int wk=1; 28 for(int k=j;k<i+j;k++,wk=(ll)wk*w%p){ 29 int x=s[k],y=(ll)s[k+i]*wk%p; 30 s[k]=(x+y)%p; 31 s[k+i]=(x-y+p)%p; 32 } 33 } 34 } 35 if(op==-1){ 36 reverse(s+1,s+n); 37 int inv=fastpow(n,p-2); 38 for(int i=0;i<n;i++){ 39 s[i]=(ll)s[i]*inv%p; 40 } 41 } 42 } 43 void GetInv(int a[],int b[],int l){ 44 if(l==1){ 45 b[0]=fastpow(a[0],p-2); 46 return; 47 } 48 GetInv(a,b,(l+1)/2); 49 int bit=1,bitnum=0; 50 for(;bit<l*2;bit<<=1)bitnum++; 51 for(int i=1;i<bit;i++){ 52 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 53 } 54 for(int i=0;i<l;i++)tmp[i]=a[i]; 55 for(int i=l;i<bit;i++)tmp[i]=0; 56 ntt(tmp,bit,1); 57 ntt(b,bit,1); 58 for(int i=0;i<bit;i++)b[i]=(ll)b[i]*((2-(ll)tmp[i]*b[i]%p+p)%p)%p; 59 ntt(b,bit,-1); 60 for(int i=l;i<bit;i++)b[i]=0; 61 } 62 int main(){ 63 scanf("%d",&n); 64 for(int i=0;i<n;i++)scanf("%d",&a[i]); 65 GetInv(a,b,n); 66 for(int i=0;i<n;i++)printf("%d ",b[i]); 67 return 0; 68 }
模板:洛谷P4238、洛谷P4239(加强版)
多项式除法(取模)
给出多项式$A(x)$,$B(x)$,求多项式$C(x)$,$D(x)$满足
$$A(x)=B(x)C(x)+D(x)$$
乍一看很难搞,于是我们需要一些玄(fa)学操作。
定义反转操作,对于一个次数为$n$的多项式,有
$$A^R(x)=x^{n}A(\frac{1}{n})$$
这个有什么用呢?这样可以把每个原本次数是$k$的项次数变为$n-k$,相当于把整个多项式的系数反转了过来。
设$A(x)$的次数为$n$,$B(x)$的次数为$m$,则$C(x)$的次数为$n-m$,$D(x)$的次数为$m-1$
用$1/x$代替$x$,两边乘以$x^n$,则:
$$x^{n}A(\frac{1}{x})=x^{m}B(\frac{1}{x})x^{n-m}C(\frac{1}{x})+x^{n-m+1}x^{m-1}D(\frac{1}{x})$$
$$A^R(x)=B^R(x)C^R(x)+x^{n-m+1}D^R(x)$$
然后惊人的事情就发生了!我们把整个式子放到模$x^{n-m+1}$意义下,有
$$A^R(x)=B^R(x)C^R(x)(\mod x^{n-m+1})$$
$$C^R(x)=A^R(x)(B^R(x))^{-1}(\mod x^{n-m+1})$$
于是就做完了。。。
时间复杂度:$O(nlogn)$
代码:
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<vector>
5 #include<cstdio>
6 #include<cmath>
7 #include<queue>
8 #include<stack>
9 #include<map>
10 #include<set>
11 using namespace std; 12 typedef long long ll; 13 const int p=998244353,g=3; 14 int n,m,a[500001],b[500001],rev[500001],aa[500001],bb[500001],tmp[500001],invb[500001],c[500001],d[500001]; 15 int fastpow(int x,int y){ 16 int ret=1; 17 for(;y;y>>=1,x=(ll)x*x%p){ 18 if(y&1)ret=(ll)ret*x%p; 19 } 20 return ret; 21 } 22 void ntt(int s[],int n,int op){ 23 for(int i=0;i<n;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 24 for(int i=1;i<n;i<<=1){ 25 int w=fastpow(g,(p-1)/(i<<1)); 26 for(int pp=i<<1,j=0;j<n;j+=pp){ 27 int wk=1; 28 for(int k=j;k<i+j;k++,wk=(ll)wk*w%p){ 29 int x=s[k],y=(ll)s[k+i]*wk%p; 30 s[k]=(x+y)%p; 31 s[k+i]=(x-y+p)%p; 32 } 33 } 34 } 35 if(op==-1){ 36 reverse(s+1,s+n); 37 int inv=fastpow(n,p-2); 38 for(int i=0;i<n;i++){ 39 s[i]=(ll)s[i]*inv%p; 40 } 41 } 42 } 43 void GetInv(int a[],int b[],int l){ 44 if(l==1){ 45 b[0]=fastpow(a[0],p-2); 46 return; 47 } 48 GetInv(a,b,(l+1)/2); 49 int bit=1,bitnum=0; 50 for(;bit<l*2;bit<<=1)bitnum++; 51 for(int i=1;i<bit;i++){ 52 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 53 } 54 for(int i=0;i<l;i++)tmp[i]=a[i]; 55 for(int i=l;i<bit;i++)tmp[i]=0; 56 ntt(tmp,bit,1); 57 ntt(b,bit,1); 58 for(int i=0;i<bit;i++)b[i]=(ll)b[i]*((2-(ll)tmp[i]*b[i]%p+p)%p)%p; 59 ntt(b,bit,-1); 60 for(int i=l;i<bit;i++)b[i]=0; 61 } 62 void GetMul(int a[],int b[],int c[],int n,int m){ 63 int bit=1,bitnum=0; 64 for(;bit<=n+m;bit<<=1)bitnum++; 65 for(int i=0;i<=n;i++)aa[i]=a[i]; 66 for(int i=0;i<=m;i++)bb[i]=b[i]; 67 for(int i=1;i<bit;i++){ 68 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 69 } 70 ntt(aa,bit,1); 71 ntt(bb,bit,1); 72 for(int i=0;i<bit;i++){ 73 c[i]=(ll)aa[i]*bb[i]%p; 74 aa[i]=bb[i]=0; 75 } 76 ntt(c,bit,-1); 77 } 78 int main(){ 79 scanf("%d%d",&n,&m); 80 for(int i=0;i<=n;i++)scanf("%d",&a[i]); 81 for(int i=0;i<=m;i++)scanf("%d",&b[i]); 82 reverse(a,a+n+1); 83 reverse(b,b+m+1); 84 GetInv(b,invb,n-m+1); 85 GetMul(a,invb,c,n-m,n-m); 86 reverse(c,c+n-m+1); 87 for(int i=0;i<=n-m;i++){ 88 printf("%d ",c[i]); 89 } 90 printf("\n"); 91 reverse(a,a+n+1); 92 reverse(b,b+m+1); 93 GetMul(c,b,d,n-m,m); 94 for(int i=0;i<m;i++)d[i]=(a[i]-d[i]+p)%p; 95 for(int i=0;i<m;i++)printf("%d ",d[i]); 96 return 0; 97 }
模板:洛谷P4512
多项式开根
给出多项式$A(x)$,求多项式$B(x)$满足
$$B^2(x)≡A(x)(\mod x^n)$$
类似于求逆,设已经求出了$B'(x)$满足
$$B'^2(x)≡A(x)(\mod x^{\lceil \frac{n}{2}\rceil})$$
则
$$B'^2(x)-A(x)≡0(\mod x^{\lceil \frac{n}{2}\rceil})$$
$$(B'^2(x)-A(x))^2≡0(\mod x^n)$$
$$B'^4(x)-2B'^2(x)A(x)+A^2(x)≡0(\mod x^n)$$
$$B'^4(x)+2B'^2(x)A(x)+A^2(x)≡4B'^2(x)A(x)(\mod x^n)$$
$$(B'^2(x)+A(x))^2≡(2B'^2(x))A(x)(\mod x^n)$$
$$(\frac{B'^2(x)+A(x)}{2B'^2(x)})^2≡A(x)(\mod x^n)$$
因此
$$B(x)=\frac{B'^2(x)+A(x)}{2B'^2(x)}$$
每次倍增向上重新赋值,时间复杂度$O(nlogn)$
代码:
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<vector>
5 #include<cstdio>
6 #include<cmath>
7 #include<queue>
8 #include<stack>
9 #include<map>
10 #include<set>
11 using namespace std; 12 typedef long long ll; 13 const int p=998244353,g=3; 14 int inv_2,n,m,bit,bitnum=0,x,a[500001],b[500001],tmp[500001],c[500001],d[500001],rev[500001]; 15 int fastpow(int x,int y){ 16 int ret=1; 17 for(;y;y>>=1,x=(ll)x*x%p){ 18 if(y&1)ret=(ll)ret*x%p; 19 } 20 return ret; 21 } 22 void ntt(int s[],int n,int op){ 23 for(int i=0;i<n;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 24 for(int i=1;i<n;i<<=1){ 25 int w=fastpow(g,(p-1)/(i<<1)); 26 for(int pp=i<<1,j=0;j<n;j+=pp){ 27 int wk=1; 28 for(int k=j;k<i+j;k++,wk=(ll)wk*w%p){ 29 int x=s[k],y=(ll)s[k+i]*wk%p; 30 s[k]=(x+y)%p; 31 s[k+i]=(x-y+p)%p; 32 } 33 } 34 } 35 if(op==-1){ 36 reverse(s+1,s+n); 37 int inv=fastpow(n,p-2); 38 for(int i=0;i<n;i++){ 39 s[i]=(ll)s[i]*inv%p; 40 } 41 } 42 } 43 void GetInv(int a[],int b[],int l){ 44 if(l==1){ 45 b[0]=fastpow(a[0],p-2); 46 return; 47 } 48 GetInv(a,b,(l+1)/2); 49 int bit=1,bitnum=0; 50 for(;bit<l*2;bit<<=1)bitnum++; 51 for(int i=1;i<bit;i++){ 52 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 53 } 54 for(int i=0;i<l;i++)tmp[i]=a[i]; 55 for(int i=l;i<bit;i++)tmp[i]=0; 56 ntt(tmp,bit,1); 57 ntt(b,bit,1); 58 for(int i=0;i<bit;i++)b[i]=(ll)b[i]*((2-(ll)tmp[i]*b[i]%p+p)%p)%p; 59 ntt(b,bit,-1); 60 for(int i=l;i<bit;i++)b[i]=0; 61 } 62 void GetSqrt(int a[],int b[],int n){ 63 if(n==1){ 64 b[0]=a[0]; 65 return; 66 } 67 GetSqrt(a,b,n/2); 68 for(int i=0;i<=n;i++)c[i]=a[i]; 69 GetInv(b,d,n); 70 int bit=1,bitnum=0; 71 for(;bit<n*2;bit<<=1)bitnum++; 72 for(int i=1;i<bit;i++){ 73 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 74 } 75 ntt(c,n*2,1); 76 ntt(d,n*2,1); 77 for(int i=0;i<n*2;i++){ 78 d[i]=(ll)d[i]*c[i]%p; 79 } 80 ntt(d,n*2,-1); 81 for(int i=0;i<n;i++){ 82 b[i]=(ll)(b[i]+d[i])%p*inv_2%p; 83 } 84 for(int i=0;i<n*2;i++)c[i]=d[i]=0; 85 } 86 int main(){ 87 scanf("%d%d",&n,&m); 88 inv_2=fastpow(2,p-2); 89 for(int i=1;i<=n;i++)scanf("%d",&x),a[x]++; 90 for(bit=1;bit<=m;bit<<=1); 91 for(int i=0;i<bit;i++){ 92 a[i]=(-4*a[i]+p)%p; 93 } 94 a[0]++; 95 GetSqrt(a,b,bit); 96 for(int i=0;i<bit;i++)a[i]=0; 97 b[0]=(b[0]+1)%p; 98 GetInv(b,c,bit); 99 for(int i=0;i<=m;i++)c[i]=(c[i]+c[i])%p; 100 for(int i=1;i<=m;i++)printf("%d\n",c[i]); 101 return 0; 102 }
模板(?):BZOJ3625小朋友和二叉树
多项式对数函数
第一次知道多项式还有相应的初等函数……(连三角函数都有……)
给出(这个我也不知道叫什么)$A(x)=\sum\limits_{i\ge 1}a_{i}x^{i}$
则有定义:
$$ln(1-A(x))=-\sum\limits_{i\ge 1}\frac{A(x)^i}{i}$$
则若有多项式$A(x)=\sum\limits_{i\ge 1}a_{i}x^{i}+1$,要求多项式$B(x)$满足$B(x)=ln(A(x))$,直接求导即可:
$$B(x)=\int\frac{A'(x)}{A(x)}$$
要求个逆元,时间复杂度$O(nlogn)$
代码:
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<vector>
5 #include<cstdio>
6 #include<cmath>
7 #include<queue>
8 #include<stack>
9 #include<map>
10 #include<set>
11 using namespace std; 12 typedef long long ll; 13 const int p=998244353,g=3; 14 int n,m,a[500001],b[500001],rev[500001],s[500001],ss[500001],tmp[500001]; 15 int fastpow(int x,int y){ 16 int ret=1; 17 for(;y;y>>=1,x=(ll)x*x%p){ 18 if(y&1)ret=(ll)ret*x%p; 19 } 20 return ret; 21 } 22 void ntt(int s[],int n,int op){ 23 for(int i=0;i<n;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 24 for(int i=1;i<n;i<<=1){ 25 int w=fastpow(g,(p-1)/(i<<1)); 26 for(int pp=i<<1,j=0;j<n;j+=pp){ 27 int wk=1; 28 for(int k=j;k<i+j;k++,wk=(ll)wk*w%p){ 29 int x=s[k],y=(ll)s[k+i]*wk%p; 30 s[k]=(x+y)%p; 31 s[k+i]=(x-y+p)%p; 32 } 33 } 34 } 35 if(op==-1){ 36 reverse(s+1,s+n); 37 int inv=fastpow(n,p-2); 38 for(int i=0;i<n;i++){ 39 s[i]=(ll)s[i]*inv%p; 40 } 41 } 42 } 43 void GetInv(int a[],int b[],int l){ 44 if(l==1){ 45 b[0]=fastpow(a[0],p-2); 46 return; 47 } 48 GetInv(a,b,(l+1)/2); 49 int bit=1,bitnum=0; 50 for(;bit<l*2;bit<<=1)bitnum++; 51 for(int i=1;i<bit;i++){ 52 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 53 } 54 for(int i=0;i<l;i++)tmp[i]=a[i]; 55 for(int i=l;i<bit;i++)tmp[i]=0; 56 ntt(tmp,bit,1); 57 ntt(b,bit,1); 58 for(int i=0;i<bit;i++)b[i]=(ll)b[i]*((2-(ll)tmp[i]*b[i]%p+p)%p)%p; 59 ntt(b,bit,-1); 60 for(int i=l;i<bit;i++)b[i]=0; 61 } 62 void GetDerivation(int a[],int b[],int n){ 63 for(int i=1;i<n;i++)b[i-1]=(ll)i*a[i]%p; 64 b[n-1]=0; 65 } 66 void GetIntegral(int a[],int b[],int n){ 67 for(int i=1;i<n;i++)b[i]=(ll)a[i-1]*fastpow(i,p-2)%p; 68 b[0]=0; 69 } 70 void Getln(int a[],int b[],int n){ 71 GetDerivation(a,s,n); 72 GetInv(a,ss,n); 73 int bit,bitnum=0; 74 for(bit=1;bit<=n;bit<<=1)bitnum++; 75 for(int i=1;i<bit;i++){ 76 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 77 } 78 ntt(s,bit,1); 79 ntt(ss,bit,1); 80 for(int i=0;i<bit;i++)s[i]=(ll)s[i]*ss[i]%p; 81 ntt(s,bit,-1); 82 GetIntegral(s,b,n); 83 } 84 int main(){ 85 int bit; 86 scanf("%d",&n); 87 for(int i=0;i<n;i++)scanf("%d",&a[i]); 88 for(bit=1;bit<=n;bit<<=1); 89 Getln(a,b,bit); 90 for(int i=0;i<n;i++)printf("%d ",b[i]); 91 return 0; 92 }
模板:洛谷P4725
多项式指数函数
(听说要用牛顿迭代法+泰勒展开?反正我只会背代码)
直接写式子:设$B(x)=exp(A(x))$,则
$$B(x)=B_{0}(x)(1-ln(B_{0}(x))+A(x))$$
其中$B_{0}(x)$表示$B(x)$的前n项,即$B_{0}(x)≡B(x)(\mod x^n)$
时间复杂度:$O(nlogn)$
代码:
1 //未验证正确性
2 #include<algorithm>
3 #include<iostream>
4 #include<cstring>
5 #include<vector>
6 #include<cstdio>
7 #include<cmath>
8 #include<queue>
9 #include<stack>
10 #include<map>
11 #include<set>
12 using namespace std; 13 typedef long long ll; 14 const int p=998244353,g=3; 15 int n,m,a[500001],b[500001],rev[500001],s[500001],ss[500001],tmp[500001],c[500001]; 16 int fastpow(int x,int y){ 17 int ret=1; 18 for(;y;y>>=1,x=(ll)x*x%p){ 19 if(y&1)ret=(ll)ret*x%p; 20 } 21 return ret; 22 } 23 void ntt(int s[],int n,int op){ 24 for(int i=0;i<n;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 25 for(int i=1;i<n;i<<=1){ 26 int w=fastpow(g,(p-1)/(i<<1)); 27 for(int pp=i<<1,j=0;j<n;j+=pp){ 28 int wk=1; 29 for(int k=j;k<i+j;k++,wk=(ll)wk*w%p){ 30 int x=s[k],y=(ll)s[k+i]*wk%p; 31 s[k]=(x+y)%p; 32 s[k+i]=(x-y+p)%p; 33 } 34 } 35 } 36 if(op==-1){ 37 reverse(s+1,s+n); 38 int inv=fastpow(n,p-2); 39 for(int i=0;i<n;i++){ 40 s[i]=(ll)s[i]*inv%p; 41 } 42 } 43 } 44 void GetInv(int a[],int b[],int l){ 45 if(l==1){ 46 b[0]=fastpow(a[0],p-2); 47 return; 48 } 49 GetInv(a,b,(l+1)/2); 50 int bit=1,bitnum=0; 51 for(;bit<l*2;bit<<=1)bitnum++; 52 for(int i=1;i<bit;i++){ 53 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 54 } 55 for(int i=0;i<l;i++)tmp[i]=a[i]; 56 for(int i=l;i<bit;i++)tmp[i]=0; 57 ntt(tmp,bit,1); 58 ntt(b,bit,1); 59 for(int i=0;i<bit;i++)b[i]=(ll)b[i]*((2-(ll)tmp[i]*b[i]%p+p)%p)%p; 60 ntt(b,bit,-1); 61 for(int i=l;i<bit;i++)b[i]=0; 62 } 63 void GetDerivation(int a[],int b[],int n){ 64 for(int i=1;i<n;i++)b[i-1]=(ll)i*a[i]%p; 65 b[n-1]=0; 66 } 67 void GetIntegral(int a[],int b[],int n){ 68 for(int i=1;i<n;i++)b[i]=(ll)a[i-1]*fastpow(i,p-2)%p; 69 b[0]=0; 70 } 71 void Getln(int a[],int b[],int n){ 72 GetDerivation(a,s,n); 73 GetInv(a,ss,n); 74 int bit,bitnum=0; 75 for(bit=1;bit<=n;bit<<=1)bitnum++; 76 for(int i=1;i<bit;i++){ 77 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 78 } 79 ntt(s,bit,1); 80 ntt(ss,bit,1); 81 for(int i=0;i<bit;i++)s[i]=(ll)s[i]*ss[i]%p; 82 ntt(s,bit,-1); 83 GetIntegral(s,b,n); 84 } 85 void Getexp(int a[],int b[],int n){ 86 if(n==1){ 87 b[0]=1; 88 return; 89 } 90 Getexp(a,b,n/2); 91 Getln(b,c,n); 92 for(int i=0;i<n;i++)c[i]=(a[i]-c[i]+p)%p; 93 c[0]++; 94 int bit=1,bitnum=0; 95 for(;bit<n*2;bit<<=1)bitnum++; 96 for(int i=1;i<bit;i++){ 97 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 98 } 99 ntt(b,bit,1); 100 ntt(c,bit,1); 101 for(int i=0;i<bit;i++){ 102 b[i]=(ll)b[i]*c[i]%p; 103 } 104 ntt(b,bit,-1); 105 for(int i=n;i<bit;i++)b[i]=0; 106 } 107 int main(){ 108 int bit; 109 scanf("%d",&n); 110 for(int i=0;i<n;i++)scanf("%d",&a[i]); 111 for(bit=1;bit<=n;bit<<=1); 112 Getexp(a,b,bit); 113 for(int i=0;i<n;i++)printf("%d ",b[i]); 114 return 0; 115 }
模板:洛谷P4726
多项式快速幂
有了上面的几个,就可以乱搞了(要注意先把最低次项$ax^d$提出来)
$$A^k(x)=exp(kln(\frac{A(x)}{ax^d}))(ax^d)^k$$
时间复杂度:$O(nlogn)$(写什么ln和exp肯定是$O(nlognlogk)快速幂啦$)
代码:
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<vector>
5 #include<cstdio>
6 #include<cmath>
7 #include<queue>
8 #include<stack>
9 #include<map>
10 #include<set>
11 using namespace std; 12 typedef long long ll; 13 const int p=1004535809,g=3; 14 int n,m,k,bit=1,bitnum=0,s[100001],tmp[100001],x,_g,rev[100001],blg[100001],inv,G[100001],a[100001],b[100001],c[100001],sum[100001],ans[100001]; 15 int fastpow(int a,int b,int p){ 16 int ret=1; 17 for(;b;b>>=1,a=(ll)a*a%p){ 18 if(b&1)ret=(ll)ret*a%p; 19 } 20 return ret; 21 } 22 bool check(int x,int n){ 23 for(int i=2;i*i<=n;i++){ 24 if((n-1)%i==0&&fastpow(x,(n-1)/i,n)==1)return false; 25 } 26 return true; 27 } 28 int GetG(int n){ 29 if(n==2)return 1; 30 for(int i=2;;i++){ 31 if(check(i,n))return i; 32 } 33 } 34 void ntt(int s[],int n,int op){ 35 for(int i=0;i<n;i++)if(i<rev[i])swap(s[i],s[rev[i]]); 36 for(int i=1;i<n;i<<=1){ 37 int w=fastpow(g,(p-1)/(i<<1),p); 38 for(int pp=i<<1,j=0;j<n;j+=pp){ 39 int wk=1; 40 for(int k=j;k<i+j;k++,wk=(ll)wk*w%p){ 41 int x=s[k],y=(ll)s[k+i]*wk%p; 42 s[k]=(x+y)%p; 43 s[k+i]=(x-y+p)%p; 44 } 45 } 46 } 47 if(op==-1){ 48 reverse(s+1,s+n); 49 int inv=fastpow(n,p-2,p); 50 for(int i=0;i<n;i++){ 51 s[i]=(ll)s[i]*inv%p; 52 } 53 } 54 } 55 void mul(int aa[],int bb[],int cc[]){ 56 for(int i=0;i<bit;i++)b[i]=aa[i],c[i]=bb[i]; 57 ntt(b,bit,1); 58 ntt(c,bit,1); 59 for(int i=0;i<bit;i++)cc[i]=(ll)b[i]*c[i]%p; 60 ntt(cc,bit,-1); 61 for(int i=m-1;i<bit;i++){ 62 cc[i-m+1]=(cc[i-m+1]+cc[i])%p; 63 cc[i]=0; 64 } 65 } 66 void pow(int a[],int n){ 67 ans[0]=1; 68 for(;n;n>>=1,mul(a,a,a)){ 69 if(n&1)mul(a,ans,ans); 70 } 71 } 72 int main(){ 73 scanf("%d%d%d%d",&n,&m,&x,&k); 74 for(int i=1;i<=k;i++)scanf("%d",&s[i]); 75 for(;bit<m*2;bit<<=1)bitnum++; 76 for(int i=0;i<bit;i++){ 77 rev[i]=(rev[i>>1]>>1|(i&1)<<(bitnum-1)); 78 } 79 _g=GetG(m); 80 //printf("%d\n",_g);
81 G[0]=1; 82 for(int i=1;i<m-1;i++){ 83 G[i]=(G[i-1]*_g)%m; 84 blg[G[i]]=i; 85 } 86 for(int i=1;i<=k;i++){ 87 if(s[i])tmp[blg[s[i]]]++; 88 } 89 pow(tmp,n); 90 printf("%d",ans[blg[x]]); 91 return 0; 92 }
模板(?):BZOJ3992 [SDOI2015]序列统计
多项式导数&多项式积分
From yww
设$A(x)=\sum_{i\ge 0}$,则$A(x)$的形式导数为:
$$A'(x)=\sum\limits_{i\ge 1}ia_{i}x^{i-1}$$
$A(x)$同上,则:
$$\int\limits A(x)=\sum\limits_{i\ge 1}\frac{a_{i-1}}{i}x^i$$
多点求值&快速插值
咕咕咕
移到了这篇博文