码迷,mamicode.com
首页 > 编程语言 > 详细

【算法总结】多项式相关

时间:2018-04-20 00:01:57      阅读:200      评论:0      收藏:0      [点我收藏+]

标签:better   sig   number   ring   sed   cal   swa   include   fft   

【快速傅里叶变换】

〖相关资料

《虚数的图解》        《虚数的意义》        《FFT学习笔记》

《从多项式乘法到快速傅里叶变换》        Fast Fourier Transform》

《【快速傅里叶变换】【FFT】【WikiOI】【P3132】【高精度练习之超大整数乘法】》

〖模板代码

  [FFT]

技术分享图片
 1 const double pi=acos(-1);
 2 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};};
 3 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
 4 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
 5 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
 6 void fft(cpx *a,int n,int f)
 7 {
 8     int k=0;while((1<<k)<n)k++;
 9     for(int i=0;i<n;i++)
10     {
11         int t=0;
12         for(int j=0;j<k;j++)
13             if(i&(1<<j))t|=(1<<(k-j-1));
14         if(i<t)swap(a[i],a[t]);
15     }
16     for(int l=2;l<=n;l<<=1)
17     {
18         int m=l>>1;
19         cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
20         for(cpx *p=a;p!=a+n;p+=l)
21         {
22             cpx w=cpx(1,0);
23             for(int i=0;i<m;i++)
24             {
25                 cpx t=p[m+i]*w;w=w*nw;
26                 p[m+i]=p[i]-t;p[i]=p[i]+t;
27             }
28         }
29     }
30 }
31 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
32 {
33     int n=1;while(n<n1+n2)n<<=1;
34     fft(c1,n,1);fft(c2,n,1);
35     for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
36     fft(c1,n,-1);
37     for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
38 }
View Code

  [MTT]

技术分享图片
 1 const double pi=acos(-1);
 2 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}
 3 w,a1[N],b1[N],a2[N],b2[N],c1[N],c2[N],c3[N],omg[N],omgi[N];
 4 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
 5 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
 6 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
 7 void fft(cpx *a,int n,int f)
 8 {
 9     int k=0;while((1<<k)<n)k++;
10     for(int i=0;i<n;i++)
11     {
12         int t=0;
13         for(int j=0;j<k;j++)
14             if(i&(1<<j))t|=(1<<(k-j-1));
15         if(i<t)swap(a[i],a[t]);
16     }
17     for(int l=2;l<=n;l<<=1)
18     {
19         int m=l>>1;
20         for(cpx *p=a;p!=a+n;p+=l)
21         {
22             for(int i=0;i<m;i++)
23             {
24                 if(f==1)w=omg[n/l*i];
25                 else w=omgi[n/l*i];
26                 cpx t=p[m+i]*w;
27                 p[m+i]=p[i]-t;p[i]=p[i]+t;
28             }
29         }
30     }
31     if(f==-1)for(int i=0;i<n;i++)a[i].r/=n;
32 }
33 struct mtt{LL s[N];}zero,st,ans;
34 int power(int a,int b,int mod)
35 {
36     int ans=1;
37     while(b)
38     {
39         if(b&1)ans=1ll*ans*a%mod;
40         a=1ll*a*a%mod;b>>=1;
41     }
42     return ans;
43 }
44 mtt operator * (mtt x,mtt y)
45 {
46     for(int i=0;i<n;i++)
47     {
48         a1[i].r=x.s[i]/qmod;a1[i].i=0;
49         a2[i].r=x.s[i]%qmod;a2[i].i=0;
50         b1[i].r=y.s[i]/qmod;b1[i].i=0;
51         b2[i].r=y.s[i]%qmod;b2[i].i=0;
52     }
53     fft(a1,n,1);fft(a2,n,1);
54     fft(b1,n,1);fft(b2,n,1);
55     tim++;
56     for(int i=0;i<n;i++)
57     {
58         c1[i]=a1[i]*b1[i];
59         c2[i]=a2[i]*b1[i]+a1[i]*b2[i];
60         c3[i]=a2[i]*b2[i];
61     }
62     fft(c1,n,-1);fft(c2,n,-1);fft(c3,n,-1);
63     mtt ret=zero;
64     for(int i=0;i<n;i++)
65     {
66         ret.s[i]=((LL)(c1[i].r+0.5)%mod*qmod%mod*qmod%mod
67                 +(LL)(c2[i].r+0.5)%mod*qmod%mod
68                 +(LL)(c3[i].r+0.5)%mod)%mod;
69     }
70     return ret;
71 }
72 int main()
73 {
74     n=1;while(n<m+m)n<<=1;
75     for(int i=0;i<n;i++)
76     {
77         omg[i]=cpx(cos(2*pi/n*i),sin(2*pi/n*i));
78         omgi[i]=cpx(cos(2*pi/n*i),sin(-2*pi/n*i));
79     }
80     return 0;
81 }
View Code

〖相关题目

1.【bzoj2179】FFT快速傅立叶

题意:给出两个n位10进制整数x和y,输出x*y。

分析:FFT裸题

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cmath>
 4 #include<cstdlib>
 5 #include<cstring>
 6 #include<iostream>
 7 using namespace std;
 8 const double pi=acos(-1);
 9 const int N=131072+5;
10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
11 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
12 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
13 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
14 int n,ans[N];
15 char ch[N];
16 void fft(cpx *a,int n,int f)
17 {
18     int k=0;while((1<<k)<n)k++;
19     for(int i=0;i<n;i++)
20     {
21         int t=0;
22         for(int j=0;j<k;j++)
23             if(i&(1<<j))t|=(1<<(k-j-1));
24         if(i<t)swap(a[i],a[t]);
25     }
26     for(int l=2;l<=n;l<<=1)
27     {
28         int m=l>>1;
29         cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
30         for(cpx *p=a;p!=a+n;p+=l)
31         {
32             cpx w=cpx(1,0);
33             for(int i=0;i<m;i++)
34             {
35                 cpx t=p[m+i]*w;w=w*nw;
36                 p[m+i]=p[i]-t;p[i]=p[i]+t;
37             }
38         }
39     }
40 }
41 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
42 {
43     int n=1;while(n<n1+n2)n<<=1;
44     fft(c1,n,1);fft(c2,n,1);
45     for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
46     fft(c1,n,-1);
47     for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
48 }
49 int main()
50 {
51     scanf("%d",&n);
52     scanf("%s",ch);
53     for(int i=0;i<n;i++)a[i].r=ch[n-i-1]-0;
54     scanf("%s",ch);
55     for(int i=0;i<n;i++)b[i].r=ch[n-i-1]-0;
56     mul(a,n,b,n,ans);int m=n*2-2;
57     for(int i=0;i<=m;i++)
58         if(ans[i]>=10)
59         {
60             ans[i+1]+=ans[i]/10;ans[i]%=10;
61             if(i==m)m++;
62         }
63     for(int i=m;i>=0;i--)printf("%d",ans[i]);
64     return 0;
65 }
66 
View Code

2.【bzoj2194】快速傅立叶之二

题意:计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n。

分析:Heret1cの博客

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cmath>
 4 #include<cstdlib>
 5 #include<cstring>
 6 #include<iostream>
 7 using namespace std;
 8 const double pi=acos(-1);
 9 const int N=262144+5;
10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
11 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
12 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
13 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
14 int n,ans[N];
15 int read()
16 {
17     int x=0,f=1;char c=getchar();
18     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
19     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
20     return x*f;
21 }
22 void fft(cpx *a,int n,int f)
23 {
24     int k=0;while((1<<k)<n)k++;
25     for(int i=0;i<n;i++)
26     {
27         int t=0;
28         for(int j=0;j<k;j++)
29             if(i&(1<<j))t|=(1<<(k-j-1));
30         if(i<t)swap(a[i],a[t]);
31     }
32     for(int l=2;l<=n;l<<=1)
33     {
34         int m=l>>1;
35         cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
36         for(cpx *p=a;p!=a+n;p+=l)
37         {
38             cpx w=cpx(1,0);
39             for(int i=0;i<m;i++)
40             {
41                 cpx t=p[m+i]*w;w=w*nw;
42                 p[m+i]=p[i]-t;p[i]=p[i]+t;
43             }
44         }
45     }
46 }
47 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
48 {
49     int n=1;while(n<n1+n2)n<<=1;
50     fft(c1,n,1);fft(c2,n,1);
51     for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
52     fft(c1,n,-1);
53     for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
54 }
55 int main()
56 {
57     n=read();
58     for(int i=0;i<n;i++)a[i]=read(),b[n-i-1]=read();
59     mul(a,n,b,n,ans);
60     for(int i=n-1;i<n*2-1;i++)printf("%d\n",ans[i]);
61     return 0;
62 }
View Code

3.【bzoj3527】[Zjoi2014]力

题意:见原题。

分析:Mychaelの博客

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cmath>
 4 #include<cstdlib>
 5 #include<cstring>
 6 #include<iostream>
 7 using namespace std;
 8 const double pi=acos(-1);
 9 const int N=262144+5;
10 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}f1[N],f2[N],g[N];
11 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
12 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
13 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
14 int n;double x,ans1[N],ans2[N];
15 void fft(cpx *a,int n,int f)
16 {
17     int k=0;while((1<<k)<n)k++;
18     for(int i=0;i<n;i++)
19     {
20         int t=0;
21         for(int j=0;j<k;j++)
22             if(i&(1<<j))t|=(1<<(k-j-1));
23         if(i<t)swap(a[i],a[t]);
24     }
25     for(int l=2;l<=n;l<<=1)
26     {
27         int m=l>>1;
28         cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
29         for(cpx *p=a;p!=a+n;p+=l)
30         {
31             cpx w=cpx(1,0);
32             for(int i=0;i<m;i++)
33             {
34                 cpx t=p[m+i]*w;w=w*nw;
35                 p[m+i]=p[i]-t;p[i]=p[i]+t;
36             }
37         }
38     }
39 }
40 int main()
41 {
42     scanf("%d",&n);
43     for(int i=0;i<n;i++)
44     {
45         scanf("%lf",&x);
46         f1[i].r=f2[n-i-1].r=x;
47     }
48     for(int i=1;i<n;i++)g[i].r=1.0/i/i;
49     int mx=1;while(mx<n*2)mx<<=1;
50     fft(f1,mx,1);fft(f2,mx,1);fft(g,mx,1);
51     for(int i=0;i<mx;i++)f1[i]=f1[i]*g[i];
52     for(int i=0;i<mx;i++)f2[i]=f2[i]*g[i];
53     fft(f1,mx,-1);fft(f2,mx,-1);
54     for(int i=0;i<n*2-1;i++)ans1[i]=f1[i].r/mx;
55     for(int i=0;i<n*2-1;i++)ans2[i]=f2[i].r/mx;
56     for(int i=0;i<n;i++)
57         printf("%.3lf\n",ans1[i]-ans2[n-i-1]);
58     return 0;
59 }
View Code

4.【bzoj3160】万径人踪灭

题意:给定一个由‘a‘和‘b‘构成的字符串,求不连续回文子序列的个数。

分析:PoPoQQQの博客

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=262144;
 8 const int M=1e6+5;
 9 const int mod=1e9+7;
10 const double pi=acos(-1);
11 int len,ans;
12 int bin[M],tmp[N],p[N];
13 char ch[M],s[N];
14 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
15 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
16 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
17 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
18 void fft(cpx *a,int n,int f)
19 {
20     int k=0;while((1<<k)<n)k++;
21     for(int i=0;i<n;i++)
22     {
23         int t=0;
24         for(int j=0;j<k;j++)
25             if(i&(1<<j))t|=(1<<(k-j-1));
26         if(i<t)swap(a[i],a[t]);
27     }
28     for(int l=2;l<=n;l<<=1)
29     {
30         int m=l>>1;
31         cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
32         for(cpx *p=a;p!=a+n;p+=l)
33         {
34             cpx w=cpx(1,0);
35             for(int i=0;i<m;i++)
36             {
37                 cpx t=p[m+i]*w;w=w*nw;
38                 p[m+i]=p[i]-t;p[i]=p[i]+t;
39             }
40         }
41     }
42 }
43 int manacher()
44 {
45     int n=0;s[++n]=#;
46     for(int i=1;i<=len;i++)s[++n]=ch[i],s[++n]=#;
47     int mxr=0,id=0,sum=0;
48     for(int i=1;i<=n;i++)
49     {
50         if(i<mxr)p[i]=min(p[2*id-i],mxr-i);else p[i]=1;
51         for(;i+p[i]<=n&&s[i-p[i]]==s[i+p[i]];p[i]++);
52         if(i+p[i]>mxr)mxr=i+p[i],id=i;
53         sum=(sum+p[i]/2)%mod;
54     }
55     return sum;
56 }
57 int main()
58 {
59     scanf("%s",ch+1);len=strlen(ch+1);
60     bin[0]=1;for(int i=1;i<=M-5;i++)bin[i]=bin[i-1]*2%mod;
61     int n=1;while(n<len*2)n<<=1;
62     for(int i=1;i<=len;i++)a[i].r=(ch[i]==a);
63     fft(a,n,1);for(int i=0;i<n;i++)b[i]=a[i]*a[i];
64     for(int i=0;i<n;i++)a[i].r=a[i].i=0;
65     for(int i=1;i<=len;i++)a[i].r=(ch[i]==b);
66     fft(a,n,1);for(int i=0;i<n;i++)b[i]=b[i]+a[i]*a[i];
67     fft(b,n,-1);for(int i=0;i<n;i++)tmp[i]=(LL)(b[i].r/n+0.5);
68     for(int i=0;i<n;i++)ans=(ans+bin[(tmp[i]+1)>>1]-1)%mod;
69     printf("%d",(ans-manacher()+mod)%mod);
70     return 0;
71 }
View Code

5.【bzoj3992】[SDOI2015]序列统计

题意:见原题。

分析:L_0_Forever_LFの博客

技术分享图片
  1 #include<cstdio>
  2 #include<cstring>
  3 #include<cmath>
  4 #include<cstdlib>
  5 #include<algorithm>
  6 #define LL long long
  7 using namespace std;
  8 const double pi=acos(-1);
  9 const int N=2e4+5;
 10 const int mod=1004535809;
 11 int n,ci,m,x,num,tmp,qmod;
 12 int lg[N],b[N];
 13 int read()
 14 {
 15     int x=0,f=1;char c=getchar();
 16     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
 17     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
 18     return x*f;
 19 }
 20 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}
 21 w,a1[N],b1[N],a2[N],b2[N],c1[N],c2[N],c3[N],omg[N],omgi[N];
 22 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
 23 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
 24 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
 25 void fft(cpx *a,int n,int f)
 26 {
 27     int k=0;while((1<<k)<n)k++;
 28     for(int i=0;i<n;i++)
 29     {
 30         int t=0;
 31         for(int j=0;j<k;j++)
 32             if(i&(1<<j))t|=(1<<(k-j-1));
 33         if(i<t)swap(a[i],a[t]);
 34     }
 35     for(int l=2;l<=n;l<<=1)
 36     {
 37         int m=l>>1;
 38         for(cpx *p=a;p!=a+n;p+=l)
 39         {
 40             for(int i=0;i<m;i++)
 41             {
 42                 if(f==1)w=omg[n/l*i];
 43                 else w=omgi[n/l*i];
 44                 cpx t=p[m+i]*w;
 45                 p[m+i]=p[i]-t;p[i]=p[i]+t;
 46             }
 47         }
 48     }
 49     if(f==-1)for(int i=0;i<n;i++)a[i].r/=n;
 50 }
 51 struct mtt{LL s[N];}zero,st,ans;
 52 int power(int a,int b,int mod)
 53 {
 54     int ans=1;
 55     while(b)
 56     {
 57         if(b&1)ans=1ll*ans*a%mod;
 58         a=1ll*a*a%mod;b>>=1;
 59     }
 60     return ans;
 61 }
 62 int find(int n)
 63 {
 64     int mx=(int)(sqrt(n)+0.5),num=n-1,tot=0;
 65     for(int i=2;i<=mx;i++)
 66         if(num%i==0)
 67         {
 68             b[++tot]=i;
 69             while(num%i==0)num/=i;
 70         }
 71     if(num!=1)b[++tot]=num;
 72     for(int i=2;i<=n;i++)
 73     {
 74         bool ok=true;
 75         for(int j=1;j<=tot;j++)
 76             if(power(i,(n-1)/b[j],n)==1){ok=false;break;}
 77         if(ok)return i;
 78     }
 79     return 0;
 80 }
 81 void pre()
 82 {
 83     int g=find(m),sum=1;
 84     for(int i=0;i<m-1;i++)
 85     {
 86         lg[sum]=i;
 87         sum=1ll*sum*g%m;
 88     }
 89 }
 90 int tim=0; 
 91 mtt operator * (mtt x,mtt y)
 92 {
 93     for(int i=0;i<n;i++)
 94     {
 95         a1[i].r=x.s[i]/qmod;a1[i].i=0;
 96         a2[i].r=x.s[i]%qmod;a2[i].i=0;
 97         b1[i].r=y.s[i]/qmod;b1[i].i=0;
 98         b2[i].r=y.s[i]%qmod;b2[i].i=0;
 99     }
100     fft(a1,n,1);fft(a2,n,1);
101     fft(b1,n,1);fft(b2,n,1);
102 //  printf("-----------%d-------\n",n);
103     tim++;
104     for(int i=0;i<n;i++)
105     {
106         c1[i]=a1[i]*b1[i];
107         c2[i]=a2[i]*b1[i]+a1[i]*b2[i];
108         c3[i]=a2[i]*b2[i];
109 //      if(tim<=10)printf("%.3lf %.3lf %.3lf\n",c1[i].r,c2[i].r,c3[i].r);
110     }
111     fft(c1,n,-1);fft(c2,n,-1);fft(c3,n,-1);
112     mtt ret=zero;
113     for(int i=0;i<n;i++)
114     {
115         LL temp=((LL)(c1[i].r+0.5)%mod*qmod%mod*qmod%mod
116                 +(LL)(c2[i].r+0.5)%mod*qmod%mod
117                 +(LL)(c3[i].r+0.5)%mod)%mod;
118         ret.s[i%(m-1)]=(ret.s[i%(m-1)]+temp)%mod;
119     }
120     return ret;
121 }
122 void work(mtt a,int b)
123 {
124     bool flag=false;
125     while(b)
126     {
127         if(b&1)
128         {
129             if(!flag)ans=a,flag=true;
130             else ans=ans*a;
131         }
132         a=a*a;b>>=1;
133 //      puts("");
134 //      for(int i=0;i<n;i++)printf("%lld ",a.s[i]);
135 //      puts("");
136     }
137 }
138 int main()
139 {
140     ci=read();m=read();x=read();num=read();
141     pre();n=1;while(n<m+m)n<<=1;
142     for(int i=0;i<n;i++)
143     {
144         omg[i]=cpx(cos(2*pi/n*i),sin(2*pi/n*i));
145         omgi[i]=cpx(cos(2*pi/n*i),sin(-2*pi/n*i));
146     }
147     for(int i=0;i<n;i++)zero.s[i]=0;
148     st=zero;
149     for(int i=1;i<=num;i++)
150     {
151         tmp=read();
152         if(tmp)st.s[lg[tmp]]=1;
153     }
154 //  for(int i=0;i<n;i++)printf("%d %d\n",i,lg[i]);
155 //  for(int i=0;i<n;i++)printf("%lld ",st.s[i]);
156 //  puts("");
157     qmod=sqrt(1.0*mod);
158     work(st,ci);
159     printf("%lld",(ans.s[lg[x]]%mod+mod)%mod);
160     return 0;
161 }
View Code

6.【bzoj4827】[Hnoi2017]礼物

题意:见原题。

分析:cdcqの博客

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm> 
 3 #include<cstring>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=262144+5;
 8 const double pi=acos(-1);
 9 int n,m,anss,sum,x[N],y[N],ans[N];
10 int read()
11 {
12     int x=0,f=1;char c=getchar();
13     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
14     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
15     return x*f;
16 }
17 struct cpx{double r,i;cpx(double _r=0,double _i=0):r(_r),i(_i){};}a[N],b[N];
18 cpx operator + (cpx a,cpx b){return cpx(a.r+b.r,a.i+b.i);}
19 cpx operator - (cpx a,cpx b){return cpx(a.r-b.r,a.i-b.i);}
20 cpx operator * (cpx a,cpx b){return cpx(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
21 void fft(cpx *a,int n,int f)
22 {
23     int k=0;while((1<<k)<n)k++;
24     for(int i=0;i<n;i++)
25     {
26         int t=0;
27         for(int j=0;j<k;j++)
28             if(i&(1<<j))t|=(1<<(k-j-1));
29         if(i<t)swap(a[i],a[t]);
30     }
31     for(int l=2;l<=n;l<<=1)
32     {
33         int m=l>>1;
34         cpx nw=cpx(cos(2*pi/l),sin(2*pi/l)*f);
35         for(cpx *p=a;p!=a+n;p+=l)
36         {
37             cpx w=cpx(1,0);
38             for(int i=0;i<m;i++)
39             {
40                 cpx t=p[m+i]*w;w=w*nw;
41                 p[m+i]=p[i]-t;p[i]=p[i]+t;
42             }
43         }
44     }
45 }
46 void mul(cpx *c1,int n1,cpx *c2,int n2,int *ans)
47 {
48     int n=1;while(n<n1+n2)n<<=1;
49     fft(c1,n,1);fft(c2,n,1);
50     for(int i=0;i<n;i++)c1[i]=c1[i]*c2[i];
51     fft(c1,n,-1);
52     for(int i=0;i<n1+n2-1;i++)ans[i]=(int)(c1[i].r/n+0.5);
53 }
54 int main()
55 {
56     n=read();m=read();
57     for(int i=0;i<n;i++)x[i]=read();
58     for(int i=0;i<n;i++)a[i].r=x[n-i-1];
59     for(int i=0;i<n;i++)y[i]=read(),b[i].r=b[i+n].r=y[i],sum+=x[i]-y[i];
60     mul(a,n,b,n*2-1,ans);
61     for(int i=n-1;i<=2*n-2;i++)anss=max(anss,ans[i]);
62     anss=-anss*2;int mx=0x3f3f3f3f;
63     for(int c=0;c<=m;c++)mx=min(mx,sum*c*2+n*c*c);
64     anss+=mx;
65     for(int i=0;i<n;i++)anss+=x[i]*x[i]+y[i]*y[i];
66     printf("%d",anss);
67     return 0;
68 }
View Code

【生成函数】

〖相关题目

1.【bzoj3028】食物

题意:给出n种重量为1的食物以及它们的限制,求出食物个数总和恰好为n的方案数。

分析:maijingの博客

技术分享图片
 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 using namespace std;
 5 const int mod=1e4+7;
 6 int n,ans;
 7 char ch[505];
 8 int main()
 9 {
10     scanf("%s",ch+1);
11     n=strlen(ch+1);
12     for(int i=1;i<=n;i++)
13         ans=(ans*10+ch[i]-0)%mod;
14     printf("%d\n",ans*(ans+1)%mod*(ans+2)%mod*1668%mod);
15     return 0;
16 }
View Code

【快速数论变换】

〖模板代码

  [NTT]

技术分享图片
 1 const int N=131072+5;
 2 const int mod=998244353;
 3 int qmod(int a,int b)
 4 {
 5     int ans=1;
 6     while(b)
 7     {
 8         if(b&1)ans=1ll*ans*a%mod;
 9         a=1ll*a*a%mod;b>>=1;
10     }
11     return ans;
12 }
13 void ntt(int *a,int n,int f)
14 {
15      int k=0;while((1<<k)<n)k++;
16      for(int i=0;i<n;i++)
17      {
18          int t=0;
19          for(int j=0;j<k;j++)
20              if(i&(1<<j))t|=(1<<(k-j-1));
21          if(i<t)swap(a[i],a[t]);
22      }
23      for(int l=2;l<=n;l<<=1)
24      {
25          int m=l>>1,nw=qmod(3,(mod-1)/l);
26          if(f==-1)nw=qmod(nw,mod-2);
27          for(int *p=a;p!=a+n;p+=l)
28          {
29              int w=1;
30              for(int i=0;i<m;i++)
31              {
32                  int t=1ll*p[m+i]*w%mod;
33                  p[m+i]=(p[i]-t+mod)%mod;
34                  p[i]=(p[i]+t)%mod;
35                  w=1ll*w*nw%mod;
36              }
37          }
38      }
39      if(f==-1)
40      {
41          int inv=qmod(n,mod-2);
42          for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
43      }
44 }
45 int main()
46 {
47     n=1;while(n<len1+len2)n<<=1;
48     ntt(a,n,1);ntt(b,n,1);
49     for(int i=0;i<n;i++)c[i]=1ll*a[i]*b[i]%mod;
50     ntt(c,n,-1);
51     return 0;
52 }
View Code

  [多项式求逆]

技术分享图片
 1 int power(int a,int b)
 2 {
 3     int ans=1;
 4     while(b)
 5     {
 6         if(b&1)ans=1ll*ans*a%mod;
 7         a=1ll*a*a%mod;b>>=1;
 8     }
 9     return ans;
10 }
11 void ntt(int *a,int n,int f)
12 {
13     int k=0;while((1<<k)<n)k++;
14     for(int i=0;i<n;i++)
15     {
16         int t=0;
17         for(int j=0;j<k;j++)
18             if(i&(1<<j))t|=(1<<(k-j-1));
19         if(i<t)swap(a[i],a[t]);
20     }
21     for(int l=2;l<=n;l<<=1)
22     {
23         int m=l>>1,nw=power(3,(mod-1)/l);
24         if(f==-1)nw=power(nw,mod-2);
25         for(int *p=a;p!=a+n;p+=l)
26         {
27             int w=1;
28             for(int i=0;i<m;i++)
29             {
30                 int t=1ll*p[m+i]*w%mod;
31                 p[m+i]=(p[i]-t+mod)%mod;
32                 p[i]=(p[i]+t)%mod;
33                 w=1ll*w*nw%mod;
34             }
35         }
36     }
37     if(f==-1)
38     {
39         int inv=power(n,mod-2);
40         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
41     }
42 }
43 void pinv(int *a,int *b,int n)
44 {
45     if(n==1){b[0]=power(a[0],mod-2);return;}
46     pinv(a,b,n>>1);n<<=1;
47     for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
48     ntt(c,n,1);ntt(b,n,1);
49     for(int i=0;i<n;i++)b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
50     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
51 }
View Code

  [多项式开根]

技术分享图片
 1 int power(int a,int b)
 2 {
 3     int ans=1;
 4     while(b)
 5     {
 6         if(b&1)ans=1ll*ans*a%mod;
 7         a=1ll*a*a%mod;b>>=1;
 8     }
 9     return ans;
10 }
11 void ntt(int *a,int n,int f)
12 {
13     int k=0;while((1<<k)<n)k++;
14     for(int i=0;i<n;i++)
15     {
16         int t=0;
17         for(int j=0;j<k;j++)
18             if(i&(1<<j))t|=(1<<(k-j-1));
19         if(i<t)swap(a[i],a[t]);
20     }
21     for(int l=2;l<=n;l<<=1)
22     {
23         int m=l>>1,nw=power(3,(mod-1)/l);
24         if(f==-1)nw=power(nw,mod-2);
25         for(int *p=a;p!=a+n;p+=l)
26         {
27             int w=1;
28             for(int i=0;i<m;i++)
29             {
30                 int t=1ll*p[m+i]*w%mod;
31                 p[m+i]=(p[i]-t+mod)%mod;
32                 p[i]=(p[i]+t)%mod;
33                 w=1ll*w*nw%mod;
34             }
35         }
36     }
37     if(f==-1)
38     {
39         int inv=power(n,mod-2);
40         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
41     }
42 }
43 void pinv(int *a,int *b,int n)
44 {
45     if(n==1){b[0]=power(a[0],mod-2);return;}
46     pinv(a,b,n>>1);n<<=1;
47     for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
48     ntt(c,n,1);ntt(b,n,1);
49     for(int i=0;i<n;i++)
50         b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
51     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
52 }
53 void psqrt(int *a,int *b,int *invb,int n)
54 {
55     if(n==1){b[0]=1;return;}
56     psqrt(a,b,invb,n>>1);
57     memset(invb,0,4*n);pinv(b,invb,n);
58     for(int i=0;i<n;i++)c[i]=a[i];
59     for(int i=n;i<n+n;i++)c[i]=0;
60     n<<=1;ntt(c,n,1);ntt(b,n,1);ntt(invb,n,1);
61     for(int i=0;i<n;i++)
62         b[i]=(b[i]+1ll*invb[i]*c[i]%mod)%mod*inv%mod;
63     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
64 }
View Code

〖相关题目

1.【bzoj2179】FFT快速傅立叶

题意:给出两个n位10进制整数x和y,输出x*y。

分析:NTT裸题

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #define LL long long
 5 using namespace std;
 6 const int N=131072+5;
 7 const int mod=998244353;
 8 int n,len,a[N],b[N],c[N];
 9 char s[N];
10 int qmod(int a,int b)
11 {
12     int ans=1;
13     while(b)
14     {
15         if(b&1)ans=1ll*ans*a%mod;
16         a=1ll*a*a%mod;b>>=1;
17     }
18     return ans;
19 }
20 void ntt(int *a,int n,int f)
21 {
22      int k=0;while((1<<k)<n)k++;
23      for(int i=0;i<n;i++)
24      {
25         int t=0;
26         for(int j=0;j<k;j++)
27             if(i&(1<<j))t|=(1<<(k-j-1));
28         if(i<t)swap(a[i],a[t]);
29      }
30      for(int l=2;l<=n;l<<=1)
31      {
32         int m=l>>1,nw=qmod(3,(mod-1)/l);
33         if(f==-1)nw=qmod(nw,mod-2);
34         for(int *p=a;p!=a+n;p+=l)
35         {
36             int w=1;
37             for(int i=0;i<m;i++)
38             {
39                 int t=1ll*p[m+i]*w%mod;
40                 p[m+i]=(p[i]-t+mod)%mod;
41                 p[i]=(p[i]+t)%mod;
42                 w=1ll*w*nw%mod;
43             }
44         }
45      }
46      if(f==-1)
47      {
48         int inv=qmod(n,mod-2);
49         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
50      }
51 }
52 int main()
53 {
54     scanf("%d",&len);
55     n=1;while(n<len*2)n<<=1;
56     scanf("%s",s+1);
57     for(int i=1;i<=len;i++)a[len-i]=s[i]-0;
58     scanf("%s",s+1);
59     for(int i=1;i<=len;i++)b[len-i]=s[i]-0;
60     ntt(a,n,1);ntt(b,n,1);
61 //  for(int i=0;i<n;i++)printf("%d %d\n",a[i],b[i]);
62     for(int i=0;i<n;i++)c[i]=1ll*a[i]*b[i]%mod;
63     ntt(c,n,-1);
64     for(int i=0;i<n;i++)
65         if(c[i]>=10)
66         {
67             c[i+1]+=c[i]/10;c[i]%=10;
68             if(i==n-1)n++;
69         }
70     n--;while(!c[n]&&n)n--;
71     for(int i=n;i>=0;i--)printf("%d",c[i]);
72     return 0;
73 }
View Code

2.【51nod1135】原根

题意:给出1个质数P,找出P最小的原根。

分析:无

技术分享图片
 1 #include<cstdio>
 2 #include<cmath>
 3 #include<cstring>
 4 #include<algorithm>
 5 using namespace std;
 6 int mod,num,mx,tot,b[1005];
 7 int qmod(int a,int b)
 8 {
 9     int ans=1;
10     while(b)
11     {
12         if(b&1)ans=1ll*ans*a%mod;
13         a=1ll*a*a%mod;b>>=1;
14     }
15     return ans;
16 }
17 int main()
18 {
19     scanf("%d",&mod);
20     mx=(int)(sqrt(mod)+0.5);num=mod-1;
21     for(int i=2;i<=mx;i++)
22         if(num%i==0)
23         {
24             b[++tot]=i;
25             while(num%i==0)num/=i;
26         }
27     if(num!=1)b[++tot]=num;
28     for(int i=2;i<=mod;i++)
29     {
30         bool ok=true;
31         for(int j=1;j<=tot;j++)
32             if(qmod(i,(mod-1)/b[j])==1){ok=false;break;}
33         if(ok){printf("%d",i);return 0;}
34     }
35     return 0;
36 }
View Code

3.【bzoj3992】[SDOI2015]序列统计

题意:见原题。

分析:WerKeyTom_FTDの博客

技术分享图片
  1 #include<cstdio>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<cmath>
  5 #include<cstdlib>
  6 using namespace std;
  7 const int N=2e4+5;
  8 const int mod=1004535809;
  9 int n,m,x,num,tmp;
 10 int lg[N],b[N],c[N],ans[N],f[N];
 11 int read()
 12 {
 13     int x=0,f=1;char c=getchar();
 14     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
 15     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
 16     return x*f;
 17 }
 18 int power(int a,int b,int mod)
 19 {
 20     int ans=1;
 21     while(b)
 22     {
 23         if(b&1)ans=1ll*ans*a%mod;
 24         a=1ll*a*a%mod;b>>=1;
 25     }
 26     return ans;
 27 }
 28 int find(int n)
 29 {
 30     int mx=(int)(sqrt(n)+0.5),num=n-1,tot=0;
 31     for(int i=2;i<=mx;i++)
 32         if(num%i==0)
 33         {
 34             b[++tot]=i;
 35             while(num%i==0)num/=i;
 36         }
 37     if(num!=1)b[++tot]=num;
 38     for(int i=2;i<=n;i++)
 39     {
 40         bool ok=true;
 41         for(int j=1;j<=tot;j++)
 42             if(power(i,(n-1)/b[j],n)==1){ok=false;break;}
 43         if(ok)return i;
 44     }
 45     return 0;
 46 }
 47 void pre()
 48 {
 49     int g=find(m),sum=1;
 50     for(int i=0;i<m-1;i++)
 51     {
 52         lg[sum]=i;
 53         sum=1ll*sum*g%m;
 54     }
 55 }
 56 void ntt(int *a,int n,int f)
 57 {
 58     int k=0;while((1<<k)<n)k++;
 59     for(int i=0;i<n;i++)
 60     {
 61         int t=0;
 62         for(int j=0;j<k;j++)
 63             if(i&(1<<j))t|=(1<<(k-j-1));
 64         if(i<t)swap(a[i],a[t]);
 65     }
 66     for(int l=2;l<=n;l<<=1)
 67     {
 68         int m=l>>1,nw=power(3,(mod-1)/l,mod);
 69         if(f==-1)nw=power(nw,mod-2,mod);
 70         for(int *p=a;p!=a+n;p+=l)
 71         {
 72             int w=1;
 73             for(int i=0;i<m;i++)
 74             {
 75                 int t=1ll*p[m+i]*w%mod;
 76                 p[m+i]=(p[i]-t+mod)%mod;
 77                 p[i]=(p[i]+t)%mod;
 78                 w=1ll*w*nw%mod;
 79             }
 80         }
 81     }
 82     if(f==-1)
 83     {
 84         int inv=power(n,mod-2,mod);
 85         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
 86     }
 87 }
 88 void mul(int *a,int *b,int n)
 89 {
 90     for(int i=0;i<n;i++)c[i]=b[i];
 91     ntt(a,n,1);ntt(c,n,1);
 92     for(int i=0;i<n;i++)a[i]=1ll*a[i]*c[i]%mod;
 93     ntt(a,n,-1);
 94     for(int i=m-1;i<n;i++)
 95         a[i%(m-1)]=(a[i%(m-1)]+a[i])%mod,a[i]=0;
 96 }
 97 void work(int *a,int b)
 98 {
 99     int n=1;while(n<m+m)n<<=1;
100     ans[0]=1;
101     while(b)
102     {
103         if(b&1)mul(ans,a,n);
104         mul(a,a,n);b>>=1;
105     }
106 }
107 int main()
108 {
109     n=read();m=read();x=read();num=read();
110     pre();
111     for(int i=1;i<=num;i++)
112     {
113         tmp=read();
114         if(tmp)f[lg[tmp]]=1;
115     }
116     work(f,n);
117     printf("%d",ans[lg[x]]);
118     return 0;
119 }
View Code

4.【bzoj4555】[Tjoi2016&Heoi2016]求和

题意:见原题。

分析:zltttttの博客

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 #include<cstdlib>
 6 using namespace std;
 7 const int N=262144+5;
 8 const int mod=998244353;
 9 int n,m,ans;
10 int fac[N],inv[N],f[N],g[N],c[N];
11 int power(int a,int b)
12 {
13     int ans=1;
14     while(b)
15     {
16         if(b&1)ans=1ll*ans*a%mod;
17         a=1ll*a*a%mod;b>>=1;
18     }
19     return ans;
20 }
21 void ntt(int *a,int n,int f)
22 {
23     int k=0;while((1<<k)<n)k++;
24     for(int i=0;i<n;i++)
25     {
26         int t=0;
27         for(int j=0;j<k;j++)
28             if(i&(1<<j))t|=(1<<(k-j-1));
29         if(i<t)swap(a[i],a[t]);
30     }
31     for(int l=2;l<=n;l<<=1)
32     {
33         int m=l>>1,nw=power(3,(mod-1)/l);
34         if(f==-1)nw=power(nw,mod-2);
35         for(int *p=a;p!=a+n;p+=l)
36         {
37             int w=1;
38             for(int i=0;i<m;i++)
39             {
40                 int t=1ll*p[m+i]*w%mod;
41                 p[m+i]=(p[i]-t+mod)%mod;
42                 p[i]=(p[i]+t)%mod;
43                 w=1ll*w*nw%mod;
44             }
45         }
46     }
47     if(f==-1)
48     {
49         int inv=power(n,mod-2);
50         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
51     }
52 }
53 void pinv(int *a,int *b,int n)
54 {
55     if(n==0){b[0]=power(a[0],mod-2);return;}
56     pinv(a,b,n>>1);n<<=1;
57     for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
58     ntt(c,n,1);ntt(b,n,1);
59     for(int i=0;i<n;i++)b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
60     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
61 }
62 int main()
63 {
64     scanf("%d",&m);
65     n=1;while(n<=m)n<<=1;
66     fac[0]=1;
67     for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
68     inv[n]=power(fac[n],mod-2);
69     for(int i=n;i>=1;i--)inv[i-1]=1ll*inv[i]*i%mod;
70     for(int i=1;i<=n;i++)f[i]=1ll*(mod-inv[i])*2%mod;
71     f[0]=1;pinv(f,g,n);
72     for(int i=0;i<=m;i++)ans=(ans+1ll*g[i]*fac[i]%mod)%mod;
73     printf("%d",ans);
74     return 0;
75 }
View Code

5.【bzoj3456】城市规划

题意:求n个点的无向简单连通图个数。

分析:PoPoQQQの博客

技术分享图片
 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 using namespace std;
 6 const int N=262144+5;
 7 const int mod=1004535809;
 8 int n,mx;
 9 int fac[N],inv[N],c[N];
10 int A[N],B[N],C[N];
11 int power(int a,int b)
12 {
13     int ans=1;
14     while(b)
15     {
16         if(b&1)ans=1ll*ans*a%mod;
17         a=1ll*a*a%mod;b>>=1;
18     }
19     return ans;
20 }
21 void ntt(int *a,int n,int f)
22 {
23     int k=0;while((1<<k)<n)k++;
24     for(int i=0;i<n;i++)
25     {
26         int t=0;
27         for(int j=0;j<k;j++)
28             if(i&(1<<j))t|=(1<<(k-j-1));
29         if(i<t)swap(a[i],a[t]);
30     }
31     for(int l=2;l<=n;l<<=1)
32     {
33         int m=l>>1,nw=power(3,(mod-1)/l);
34         if(f==-1)nw=power(nw,mod-2);
35         for(int *p=a;p!=a+n;p+=l)
36         {
37             int w=1;
38             for(int i=0;i<m;i++)
39             {
40                 int t=1ll*p[m+i]*w%mod;
41                 p[m+i]=(p[i]-t+mod)%mod;
42                 p[i]=(p[i]+t)%mod;
43                 w=1ll*w*nw%mod;
44             }
45         }
46     }
47     if(f==-1)
48     {
49         int inv=power(n,mod-2);
50         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
51     }
52 }
53 void pinv(int *a,int *b,int n)
54 {
55     if(n==1){b[0]=power(a[0],mod-2);return;}
56     pinv(a,b,n>>1);n<<=1;
57     for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
58     ntt(c,n,1);ntt(b,n,1);
59     for(int i=0;i<n;i++)
60         b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
61     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
62 }
63 int main()
64 {
65     scanf("%d",&n);
66     mx=1;while(mx<=n)mx<<=1;
67     fac[0]=1;
68     for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
69     inv[n]=power(fac[n],mod-2);
70     for(int i=n;i>=1;i--)inv[i-1]=1ll*inv[i]*i%mod;
71     for(int i=0;i<=n;i++)
72         C[i]=1ll*power(2,1ll*i*(i-1)/2%(mod-1))*inv[i]%mod;
73     pinv(C,A,mx);
74 //  for(int i=1;i<=n;i++)printf("%d\n",A[i]);
75     for(int i=1;i<=n;i++)
76         B[i]=1ll*power(2,1ll*i*(i-1)/2%(mod-1))*inv[i-1]%mod;
77     mx<<=1;ntt(A,mx,1);ntt(B,mx,1);
78     memset(C,0,sizeof(C));
79     for(int i=0;i<mx;i++)C[i]=1ll*A[i]*B[i]%mod;
80     ntt(C,mx,-1);
81     printf("%lld",1ll*C[n]*fac[n-1]%mod);
82     return 0;
83 }
View Code

6.【bzoj3625】[Codeforces Round #250]小朋友和二叉树

题意:见原题。

分析:wzq_QwQの博客

技术分享图片
 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int N=262144+5;
 8 const int mod=998244353;
 9 int n,num,m,x,inv;
10 int f[N],g[N],invg[N],c[N];
11 int read()
12 {
13     int x=0,f=1;char c=getchar();
14     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
15     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
16     return x*f;
17 }
18 int power(int a,int b)
19 {
20     int ans=1;
21     while(b)
22     {
23         if(b&1)ans=1ll*ans*a%mod;
24         a=1ll*a*a%mod;b>>=1;
25     }
26     return ans;
27 }
28 void ntt(int *a,int n,int f)
29 {
30     int k=0;while((1<<k)<n)k++;
31     for(int i=0;i<n;i++)
32     {
33         int t=0;
34         for(int j=0;j<k;j++)
35             if(i&(1<<j))t|=(1<<(k-j-1));
36         if(i<t)swap(a[i],a[t]);
37     }
38     for(int l=2;l<=n;l<<=1)
39     {
40         int m=l>>1,nw=power(3,(mod-1)/l);
41         if(f==-1)nw=power(nw,mod-2);
42         for(int *p=a;p!=a+n;p+=l)
43         {
44             int w=1;
45             for(int i=0;i<m;i++)
46             {
47                 int t=1ll*p[m+i]*w%mod;
48                 p[m+i]=(p[i]-t+mod)%mod;
49                 p[i]=(p[i]+t)%mod;
50                 w=1ll*w*nw%mod;
51             }
52         }
53     }
54     if(f==-1)
55     {
56         int inv=power(n,mod-2);
57         for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
58     }
59 }
60 void pinv(int *a,int *b,int n)
61 {
62     if(n==1){b[0]=power(a[0],mod-2);return;}
63     pinv(a,b,n>>1);n<<=1;
64     for(int i=0;i<n/2;i++)c[i]=a[i],c[i+n/2]=0;
65     ntt(c,n,1);ntt(b,n,1);
66     for(int i=0;i<n;i++)
67         b[i]=(2*b[i]%mod-1ll*c[i]*b[i]%mod*b[i]%mod+mod)%mod;
68     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
69 }
70 void psqrt(int *a,int *b,int *invb,int n)
71 {
72     if(n==1){b[0]=1;return;}
73     psqrt(a,b,invb,n>>1);
74     memset(invb,0,4*n);pinv(b,invb,n);
75     for(int i=0;i<n;i++)c[i]=a[i];
76     for(int i=n;i<n+n;i++)c[i]=0;
77     n<<=1;ntt(c,n,1);ntt(b,n,1);ntt(invb,n,1);
78     for(int i=0;i<n;i++)
79         b[i]=(b[i]+1ll*invb[i]*c[i]%mod)%mod*inv%mod;
80     ntt(b,n,-1);for(int i=n/2;i<n;i++)b[i]=0;
81 }
82 int main()
83 {
84     num=read();m=read();inv=power(2,mod-2);
85     for(int i=1;i<=num;i++)x=read(),f[x]=1;
86     for(int i=1;i<=m;i++)f[i]=(-f[i]*4%mod+mod)%mod;
87     n=1;while(n<=m)n<<=1;
88     f[0]=1;psqrt(f,g,invg,n);
89     g[0]=(g[0]+1)%mod;
90     memset(f,0,sizeof(f));
91     pinv(g,f,n);
92     for(int i=1;i<=m;i++)printf("%d\n",f[i]*2%mod);
93     return 0;
94 }
View Code

【拉格朗日插值】

〖相关题目

1.【bzoj4559】成绩比较

题意:见原题。

分析:cdsszjjの博客

技术分享图片
 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long
 6 using namespace std;
 7 const int mod=1e9+7;
 8 const int N=105;
 9 int n,m,k;
10 int u[N],r[N],g[N],w[N];
11 int c[N][N],f[N][N];
12 int read()
13 {
14     int x=0,f=1;char c=getchar();
15     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
16     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
17     return x*f;
18 }
19 int power(int a,int b)
20 {
21     int ans=1;
22     while(b)
23     {
24         if(b&1)ans=1ll*ans*a%mod;
25         a=1ll*a*a%mod;b>>=1;
26     }
27     return ans;
28 }
29 int lag(int u,int r)
30 {
31     memset(g,0,sizeof(g));
32     for(int x=1;x<=n+1;x++)
33     {
34         for(int i=1;i<=x;i++)
35             g[x]=(g[x]+1ll*power(x-i,r-1)*power(i,n-r)%mod)%mod;
36         if(u==x)return g[x];
37     }
38     for(int i=1;i<=n+1;i++)
39     {
40         w[i]=1;
41         for(int j=1;j<=n+1;j++)
42             if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
43         w[i]=power(w[i],mod-2);
44     }
45     int sum=0;
46     for(int i=1;i<=n+1;i++)
47     {
48         int tmp=1ll*w[i]*g[i]%mod;
49         for(int j=1;j<=n+1;j++)
50             if(i!=j)tmp=1ll*tmp*(u-j)%mod;
51         sum=(sum+tmp)%mod;
52     }
53     return sum;
54 }
55 int C(int n,int m)
56 {
57     if(n<0||m<0||n<m)return 0;
58     return c[n][m];
59 }
60 int main()
61 {
62     n=read();m=read();k=read();
63     for(int i=1;i<=m;i++)u[i]=read();
64     for(int i=1;i<=m;i++)r[i]=read();
65     for(int i=0;i<=n;i++)c[i][0]=1;
66     for(int i=1;i<=n;i++)
67         for(int j=1;j<=i;j++)
68             c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
69     f[0][0]=1;
70     for(int i=1;i<=m;i++)
71     {
72         int tmp=lag(u[i],r[i]);
73         for(int j=0;j<=n-1;j++)
74         {
75             for(int w=0;w<=j;w++)
76                 f[i][j]=(f[i][j]+1ll*C(n-1-w,j-w)*C(w,r[i]-1-j+w)%mod*f[i-1][w]%mod)%mod;
77             f[i][j]=1ll*f[i][j]*tmp%mod;
78         }
79     }
80     printf("%d",f[m][n-1-k]);
81     return 0;
82 }
View Code

2.【bzoj2655】calc

题意:见原题。

分析:ez_ywwの博客

技术分享图片
 1 #include<cstdio>
 2 #include<algorithm>
 3 #include<cstring>
 4 #include<cmath>
 5 using namespace std;
 6 const int N=505;
 7 int m,n,mod,ans;
 8 int w[N],f[N*2][N];
 9 int power(int a,int b)
10 {
11     int ans=1;
12     while(b)
13     {
14         if(b&1)ans=1ll*ans*a%mod;
15         a=1ll*a*a%mod;b>>=1;
16     }
17     return ans;
18 }
19 int main()
20 {
21     scanf("%d%d%d",&m,&n,&mod);
22     f[0][0]=1;
23     for(int i=1;i<=2*n;i++)
24     {
25         f[i][0]=f[i-1][0];
26         for(int j=1;j<=n;j++)
27             f[i][j]=(1ll*f[i-1][j-1]*i%mod*j%mod+f[i-1][j])%mod;
28     }
29     if(m<=2*n){printf("%d",f[m][n]);return 0;}
30     for(int i=0;i<=2*n;i++)
31     {
32         w[i]=1;
33         for(int j=0;j<=2*n;j++)
34             if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
35         w[i]=power(w[i],mod-2);
36     }
37     for(int i=0;i<=2*n;i++)
38     {
39         int sum=1ll*f[i][n]*w[i]%mod;
40         for(int j=0;j<=2*n;j++)
41             if(i!=j)sum=1ll*sum*(m-j)%mod;
42         ans=(ans+sum)%mod;
43     }
44     printf("%d\n",ans);
45     return 0;
46 }
View Code

3.【bzoj3453】tyvj 1858 XLkxc

题意:见原题。

分析:DZYOの博客

技术分享图片
 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #include<iostream>
 6 using namespace std;
 7 const int N=150;
 8 const int mod=1234567891;
 9 int T,k,a,n,d,m,sum,ans;
10 int g[N],w[N],f[N];
11 int read()
12 {
13     int x=0,f=1;char c=getchar();
14     while(c<0||c>9){if(c==-)f=-1;c=getchar();}
15     while(c>=0&&c<=9){x=x*10+c-0;c=getchar();}
16     return x*f;
17 }
18 int power(int a,int b)
19 {
20     int ans=1;
21     while(b)
22     {
23         if(b&1)ans=1ll*ans*a%mod;
24         a=1ll*a*a%mod;b>>=1;
25     }
26     return ans;
27 }
28 void work()
29 {
30     k=read();a=read();n=read();d=read();
31     for(int i=1;i<=k+3;i++)g[i]=power(i,k);
32     for(int i=2;i<=k+3;i++)g[i]=(1ll*g[i-1]+g[i])%mod;
33     for(int i=2;i<=k+3;i++)g[i]=(1ll*g[i-1]+g[i])%mod;
34     for(int i=1;i<=k+3;i++)
35     {
36         w[i]=1;
37         for(int j=1;j<=k+3;j++)
38             if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
39         w[i]=power(w[i],mod-2);
40     }
41     for(int x=0;x<=k+4;x++)
42     {
43         f[x]=0;m=(a+1ll*x*d%mod)%mod;
44         for(int i=1;i<=k+3;i++)
45         {
46             sum=1ll*g[i]*w[i]%mod;
47             for(int j=1;j<=k+3;j++)
48                 if(i!=j)sum=(1ll*sum*(m-j)%mod+mod)%mod;
49             f[x]=(1ll*f[x]+sum)%mod;
50         }
51     }
52     for(int i=1;i<=k+4;i++)f[i]=(1ll*f[i-1]+f[i])%mod;
53     for(int i=0;i<=k+4;i++)
54     {
55         w[i]=1;
56         for(int j=0;j<=k+4;j++)
57             if(i!=j)w[i]=(1ll*w[i]*(i-j)%mod+mod)%mod;
58         w[i]=power(w[i],mod-2);
59     }
60     ans=0;
61     for(int i=0;i<=k+4;i++)
62     {
63         sum=1ll*f[i]*w[i]%mod;
64         for(int j=0;j<=k+4;j++)
65             if(i!=j)sum=(1ll*sum*(n-j)%mod+mod)%mod;
66         ans=(1ll*ans+sum)%mod;
67     }
68     printf("%d\n",ans);
69 }
70 int main()
71 {
72     T=read();
73     while(T--)work();
74     return 0;
75 }
View Code

                                                                                                                                                                                                                   

【算法总结】多项式相关

标签:better   sig   number   ring   sed   cal   swa   include   fft   

原文地址:https://www.cnblogs.com/zsnuo/p/8521067.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!