天天看點

多項式全家桶一些定義多項式加減法多項式乘法多項式求逆多項式除法(取模)多項式開根多項式對數函數多項式指數函數多項式快速幂多項式導數&多項式積分多點求值&快速插值

開頭Orz yww,Orz xfz,Orz dtz,Orz lyy

部分轉載自yww的部落格、dtz的部落格

一些定義

一個關于$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$$

多點求值&快速插值

咕咕咕

移到了這篇博文

轉載于:https://www.cnblogs.com/dcdcbigbig/p/9359329.html