多项式全家桶


开头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$$

多点求值&快速插值

咕咕咕

移到了这篇博文


免责声明!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系本站邮箱yoyou2525@163.com删除。



 
粤ICP备18138465号  © 2018-2025 CODEPRJ.COM