zoukankan      html  css  js  c++  java
  • 多项式全家桶

    开头Orz yww,Orz xfz,Orz dtz,Orz lyy

    部分转载自yww的博客dtz的博客

    一些定义

    一个关于$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$$

    多点求值&快速插值

    咕咕咕

    移到了这篇博文

  • 相关阅读:
    利用python脚本统计和删除redis key
    利用expect交互完成多台linux主机ssh key推送
    iptables -L很慢的原因
    tomcat各个端口的作用
    rabbitmq集群搭建
    ping 没有回icmp reply
    go mod 无法下载依赖问题
    0/1 nodes are available: 1 node(s) had taint
    go 编译:build constraints exclude all Go files in
    k8s单机部署
  • 原文地址:https://www.cnblogs.com/dcdcbigbig/p/9359329.html
Copyright © 2011-2022 走看看