开头Orz yww,Orz xfz,Orz dtz,Orz lyy
一些定义
一个关于$x$的多项式$A(x)$可以表示为如下形式和:
$$A(x)=sumlimits_{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} ceil})$$
则
$$A(x)B'(x)-1≡0(mod x^{lceil frac{n}{2}
ceil})$$
$$(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)=sumlimits_{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("
");
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} ceil})$$
则
$$B'^2(x)-A(x)≡0(mod x^{lceil frac{n}{2} ceil})$$
$$(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
",c[i]);
101 return 0;
102 }
模板(?):BZOJ3625小朋友和二叉树
多项式对数函数
第一次知道多项式还有相应的初等函数……(连三角函数都有……)
给出(这个我也不知道叫什么)$A(x)=sumlimits_{ige 1}a_{i}x^{i}$
则有定义:
$$ln(1-A(x))=-sumlimits_{ige 1}frac{A(x)^i}{i}$$
则若有多项式$A(x)=sumlimits_{ige 1}a_{i}x^{i}+1$,要求多项式$B(x)$满足$B(x)=ln(A(x))$,直接求导即可:
$$B(x)=intfrac{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
",_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_{ige 0}$,则$A(x)$的形式导数为:
$$A'(x)=sumlimits_{ige 1}ia_{i}x^{i-1}$$
$A(x)$同上,则:
$$intlimits A(x)=sumlimits_{ige 1}frac{a_{i-1}}{i}x^i$$
多点求值&快速插值
咕咕咕
移到了这篇博文