关于万恶的多项式

时间：2018-08-12 22:22:42 阅读：179 评论：0 收藏：0 [点我收藏+]

模板：

 1 //Achen
 2 #include<algorithm>
 3 #include<iostream>
 4 #include<cstring>
 5 #include<cstdlib>
 6 #include<vector>
 7 #include<cstdio>
 8 #include<queue>
 9 #include<cmath>
10 #include<set>
11 #include<map>
12 #define pi acos(-1)
13 #define Formylove return 0
14 #define For(i,a,b) for(int i=(a);i<=(b);i++)
15 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
16 const int N=2097155; 
17 typedef long long LL;
18 typedef double db;
19 using namespace std;
20 int n,m,ans[N];
21 
22 template<typename T>void read(T &x)  {
23     char ch=getchar(); x=0; T f=1;
24     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
25     if(ch==‘-‘) f=-1,ch=getchar();
26     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
27 }
28 
29 struct E {
30     db a,b;
31     E(db a=0,db b=0):a(a),b(b){}
32 }a[N],b[N];
33 E operator +(const E&A,const E&B) { return E(A.a+B.a,A.b+B.b); }
34 E operator -(const E&A,const E&B) { return E(A.a-B.a,A.b-B.b); }
35 E operator *(const E&A,const E&B) { return E(A.a*B.a-A.b*B.b,A.a*B.b+A.b*B.a); }
36 E operator /(const E&A,const db&B) { return E(A.a/B,A.b/B); }
37 
38 int rev[N];
39 void FFT(E a[],int n,int l,int f) {
40     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
41     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
42     for(int i=1;i<n;i<<=1) {
43         E wi(cos(pi/i),f*sin(pi/i));
44         for(int j=0,pp=(i<<1);j<n;j+=pp) {
45             E w(1,0);
46             for(int k=0;k<i;k++,w=w*wi) {
47                 E x=a[j+k],y=w*a[j+k+i];
48                 a[j+k]=x+y; a[j+k+i]=x-y;
49             }
50         }
51     }
52     if(f==-1) For(i,0,n) a[i].a=(int)(a[i].a/n+0.5);
53 }
54 
55 void mul(E a[],E b[],int c[],int n,int m) {
56     static E A[N],B[N];
57     int nn,l;
58     for(nn=1,l=0;nn<=n+m;nn<<=1) l++;
59     For(i,0,n) A[i]=a[i]; For(i,n+1,nn) A[i]=E(0,0);
60     For(i,0,m) B[i]=b[i]; For(i,m+1,nn) B[i]=E(0,0);
61     FFT(A,nn,l,1); FFT(B,nn,l,1);
62     For(i,0,nn) A[i]=A[i]*B[i];
63     FFT(A,nn,l,-1);
64     For(i,0,n+m) c[i]=A[i].a;
65 }
66 
67 int main() {
68 #ifdef ANS
69     freopen(".in","r",stdin);
70     freopen(".out","w",stdout);
71 #endif
72     read(n); read(m);
73     For(i,0,n) read(a[i].a);
74     For(i,0,m) read(b[i].a);
75     mul(a,b,ans,n,m);
76     For(i,0,n+m) printf("%d ",ans[i]);
77     Formylove;
78 }

FFT

 1 //Achen
 2 #include<algorithm>
 3 #include<iostream>
 4 #include<cstring>
 5 #include<cstdlib>
 6 #include<vector>
 7 #include<cstdio>
 8 #include<queue>
 9 #include<cmath>
10 #include<set>
11 #include<map>
12 #define Formylove return 0
13 #define For(i,a,b) for(int i=(a);i<=(b);i++)
14 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
15 const int N=2097155,p=998244353,g=3,gi=332748118; 
16 typedef long long LL;
17 typedef double db;
18 using namespace std;
19 int n,m;
20 LL a[N],b[N],ans[N];
21 
22 template<typename T>void read(T &x)  {
23     char ch=getchar(); x=0; T f=1;
24     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
25     if(ch==‘-‘) f=-1,ch=getchar();
26     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
27 }
28 
29 LL ksm(LL a,LL b) {
30     LL rs=1,bs=a;
31     while(b) {
32         if(b&1) rs=rs*bs%p;
33         bs=bs*bs%p;
34         b>>=1;
35     }
36     return rs;
37 }
38 
39 LL mo(LL x) { return x<0?x+p:(x>=p?x-p:x); }
40 
41 int rev[N];
42 void FNT(LL a[],int n,int l,int f) {
43     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
44     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
45     for(int i=1;i<n;i<<=1) {
46         LL wi=ksm((f==1)?g:gi,(p-1)/(i<<1));
47         for(int j=0,pp=(i<<1);j<n;j+=pp) {
48             LL w=1;
49             for(int k=0;k<i;k++,w=w*wi%p) {
50                 LL x=a[j+k],y=w*a[j+k+i]%p;
51                 a[j+k]=mo(x+y); a[j+k+i]=mo(x-y);
52             }
53         }
54     }
55     if(f==-1) {
56         LL inv_n=ksm(n,p-2);
57         For(i,0,n) a[i]=a[i]*inv_n%p;
58     }
59 }
60 
61 void mul(LL a[],LL b[],LL c[],int n,int m) {
62     static LL A[N],B[N];
63     int nn,l;
64     for(nn=1,l=0;nn<=n+m;nn<<=1) l++;
65     For(i,0,n) A[i]=a[i]; For(i,n+1,nn) A[i]=0;
66     For(i,0,m) B[i]=b[i]; For(i,m+1,nn) B[i]=0;
67     FNT(A,nn,l,1); FNT(B,nn,l,1);
68     For(i,0,nn) A[i]=A[i]*B[i]%p;
69     FNT(A,nn,l,-1);
70     For(i,0,n+m) c[i]=A[i];
71 }
72 
73 int main() {
74 #ifdef ANS
75     freopen(".in","r",stdin);
76     freopen(".out","w",stdout);
77 #endif
78     read(n); read(m);
79     For(i,0,n) read(a[i]);
80     For(i,0,m) read(b[i]);
81     mul(a,b,ans,n,m);
82     For(i,0,n+m) printf("%lld ",ans[i]);
83     Formylove;
84 }

FNT

 1 //易错:get_inv中传入的b必须是空的,空间开到 get_inv中的n的2倍 
 2 //Achen
 3 #include<algorithm>
 4 #include<iostream>
 5 #include<cstring>
 6 #include<cstdlib>
 7 #include<vector>
 8 #include<cstdio>
 9 #include<queue>
10 #include<cmath>
11 #include<set>
12 #include<map>
13 #define Formylove return 0
14 #define For(i,a,b) for(int i=(a);i<=(b);i++)
15 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
16 const int N=262155,p=998244353,g=3,gi=332748118; 
17 typedef long long LL;
18 typedef double db;
19 using namespace std;
20 int n;
21 LL a[N],ans[N];
22 
23 template<typename T>void read(T &x)  {
24     char ch=getchar(); x=0; T f=1;
25     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
26     if(ch==‘-‘) f=-1,ch=getchar();
27     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
28 }
29  
30 LL ksm(LL a,LL b) {
31     LL rs=1,bs=a;
32     while(b) {
33         if(b&1) rs=rs*bs%p;
34         bs=bs*bs%p;
35         b>>=1;
36     }
37     return rs;
38 }
39 
40 LL mo(LL x) { return x<0?x+p:(x>=p?x-p:x); }
41 
42 int rev[N];
43 void FNT(LL a[],int n,int l,int f) {
44     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
45     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
46     for(int i=1;i<n;i<<=1) {
47         LL wi=ksm((f==1)?g:gi,(p-1)/(i<<1));
48         for(int j=0,pp=(i<<1);j<n;j+=pp) {
49             LL w=1;
50             for(int k=0;k<i;k++,w=w*wi%p) {
51                 LL x=a[j+k],y=w*a[j+k+i]%p;
52                 a[j+k]=mo(x+y); a[j+k+i]=mo(x-y);
53             }
54         }
55     }
56     if(f==-1) {
57         LL inv_n=ksm(n,p-2);
58         For(i,0,n) a[i]=a[i]*inv_n%p;
59     }
60 }
61 
62 LL tp[N];
63 void get_inv(LL a[],LL b[],int n,int l) {
64     if(n==0) {
65         b[0]=ksm(a[0],p-2);
66         return ;
67     }
68     get_inv(a,b,n>>1,l-1);
69     For(i,0,n-1) tp[i]=a[i],tp[i+n]=0;
70     FNT(tp,n<<1,l+1,1);
71     FNT(b,n<<1,l+1,1);
72     For(i,0,(n<<1)) tp[i]=mo(2LL-tp[i]*b[i]%p)*b[i]%p;
73     FNT(tp,n<<1,l+1,-1);
74     For(i,0,n-1) b[i]=tp[i],b[i+n]=0;
75 }
76 
77 int main() {
78 #ifdef ANS
79     freopen(".in","r",stdin);
80     freopen(".out","w",stdout);
81 #endif
82     read(n);
83     For(i,0,n-1) read(a[i]);
84     int nn,l;
85     for(nn=1,l=0;nn<=n;nn<<=1) l++;
86     get_inv(a,ans,nn,l);
87     For(i,0,n-1) printf("%lld ",ans[i]);
88     Formylove;
89 }

多项式求逆

 1 //Achen
 2 #include<algorithm>
 3 #include<iostream>
 4 #include<cstring>
 5 #include<cstdlib>
 6 #include<vector>
 7 #include<cstdio>
 8 #include<queue>
 9 #include<cmath>
10 #include<set>
11 #include<map>
12 #define Formylove return 0
13 #define For(i,a,b) for(int i=(a);i<=(b);i++)
14 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
15 const int N=262155,p=998244353,g=3,gi=332748118; 
16 typedef long long LL;
17 typedef double db;
18 using namespace std;
19 int n;
20 LL a[N],b[N],tp[N];
21 
22 template<typename T>void read(T &x)  {
23     char ch=getchar(); x=0; T f=1;
24     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
25     if(ch==‘-‘) f=-1,ch=getchar();
26     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
27 }
28 
29 LL ksm(LL a,LL b) {
30     LL rs=1,bs=a;
31     while(b) {
32         if(b&1) rs=rs*bs%p;
33         bs=bs*bs%p;
34         b>>=1; 
35     }    
36     return rs;
37 }
38 
39 LL mo(LL x) { return x<0?x+p:(x>=p?x-p:x); }
40 int rev[N];
41 void FNT(LL a[],int n,int l,int f) {
42     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
43     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
44     for(int i=1;i<n;i<<=1) {
45         LL wi=ksm((f==1?g:gi),(p-1)/(i<<1));
46         for(int j=0,pp=(i<<1);j<n;j+=pp) {
47             LL w=1;
48             for(int k=0;k<i;k++,w=w*wi%p) {
49                 LL x=a[j+k],y=w*a[j+k+i]%p;
50                 a[j+k]=mo(x+y); a[j+k+i]=mo(x-y);
51             }
52         } 
53     }
54     if(f==-1) {
55         LL inv_n=ksm(n,p-2);
56         For(i,0,n) a[i]=a[i]*inv_n%p;
57     }
58 }
59 
60 void get_inv(LL a[],LL b[],LL tp[],int n,int l) {
61     if(n==1) { b[0]=ksm(a[0],p-2); return; }
62     get_inv(a,b,tp,(n>>1),l-1);
63     For(i,0,n-1) tp[i]=a[i],tp[i+n]=0;
64     FNT(tp,(n<<1),l+1,1);
65     FNT(b,(n<<1),l+1,1);
66     For(i,0,(n<<1)) tp[i]=mo(2LL-tp[i]*b[i]%p)*b[i]%p;
67     FNT(tp,(n<<1),l+1,-1);
68     For(i,0,n-1) b[i]=tp[i],b[i+n]=0;
69 }
70 
71 LL inv[N];
72 void get_ln(LL a[],LL b[],int n,int l) {
73     static LL a_dao[N],a_inv[N];
74     For(i,0,(n<<1)) a_inv[i]=a_dao[i]=0;
75     For(i,0,n-1) a_dao[i]=a[i+1]*(i+1)%p;
76     get_inv(a,a_inv,tp,n,l);
77     FNT(a_inv,(n<<1),l+1,1);
78     FNT(a_dao,(n<<1),l+1,1);
79     For(i,0,(n<<1)) a_inv[i]=a_inv[i]*a_dao[i]%p;
80     FNT(a_inv,(n<<1),l+1,-1);
81     b[0]=0;
82     For(i,1,n-1) b[i]=a_inv[i-1]*inv[i]%p;
83 }
84 
85 int main() {
86 #ifdef ANS
87     freopen(".in","r",stdin);
88     freopen(".out","w",stdout);
89 #endif
90     read(n);
91     For(i,0,n-1) read(a[i]);
92     int nn,l;
93     for(nn=1,l=0;nn<=n;nn<<=1) l++;
94     inv[0]=inv[1]=1;
95     For(i,2,n) inv[i]=mo(p-p/i*inv[p%i]%p);
96     get_ln(a,b,nn,l);
97     For(i,0,n-1) printf("%lld ",b[i]);
98     Formylove;
99 }

多项式对数函数(ln)

  1 //Achen
  2 #include<algorithm>
  3 #include<iostream>
  4 #include<cstring>
  5 #include<cstdlib>
  6 #include<vector>
  7 #include<cstdio>
  8 #include<queue>
  9 #include<cmath>
 10 #include<set>
 11 #include<map>
 12 #define Formylove return 0
 13 #define For(i,a,b) for(int i=(a);i<=(b);i++)
 14 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
 15 const int N=262155,p=998244353,g=3,gi=332748118; 
 16 typedef long long LL;
 17 typedef double db;
 18 using namespace std;
 19 int n;
 20 LL a[N],b[N];
 21 
 22 template<typename T>void read(T &x)  {
 23     char ch=getchar(); x=0; T f=1;
 24     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
 25     if(ch==‘-‘) f=-1,ch=getchar();
 26     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
 27 }
 28 
 29 LL ksm(LL a,LL b) {
 30     LL rs=1,bs=a;
 31     while(b) {
 32         if(b&1) rs=rs*bs%p;
 33         bs=bs*bs%p;
 34         b>>=1; 
 35     }    
 36     return rs;
 37 }
 38 
 39 LL mo(LL x) { return x<0?x+p:(x>=p?x-p:x); }
 40 int rev[N];
 41 void FNT(LL a[],int n,int l,int f) {
 42     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
 43     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
 44     for(int i=1;i<n;i<<=1) {
 45         LL wi=ksm((f==1?g:gi),(p-1)/(i<<1));
 46         for(int j=0,pp=(i<<1);j<n;j+=pp) {
 47             LL w=1;
 48             for(int k=0;k<i;k++,w=w*wi%p) {
 49                 LL x=a[j+k],y=w*a[j+k+i]%p;
 50                 a[j+k]=mo(x+y); a[j+k+i]=mo(x-y);
 51             }
 52         } 
 53     }
 54     if(f==-1) {
 55         LL inv_n=ksm(n,p-2);
 56         For(i,0,n) a[i]=a[i]*inv_n%p;
 57     }
 58 }
 59 
 60 LL tp[N];
 61 void get_inv(LL a[],LL b[],LL tp[],int n,int l) {
 62     if(n==1) { b[0]=ksm(a[0],p-2); return; }
 63     get_inv(a,b,tp,(n>>1),l-1);
 64     For(i,0,n-1) tp[i]=a[i],tp[i+n]=0;
 65     FNT(tp,(n<<1),l+1,1);
 66     FNT(b,(n<<1),l+1,1);
 67     For(i,0,(n<<1)) tp[i]=mo(2LL-tp[i]*b[i]%p)*b[i]%p;
 68     FNT(tp,(n<<1),l+1,-1);
 69     For(i,0,n-1) b[i]=tp[i],b[i+n]=0;
 70 }
 71 
 72 LL inv[N];
 73 void get_ln(LL a[],LL b[],int n,int l) {
 74     static LL a_dao[N],a_inv[N];
 75     For(i,0,(n<<1)) a_inv[i]=a_dao[i]=0;
 76     For(i,0,n-1) a_dao[i]=a[i+1]*(i+1)%p; 
 77     a_dao[n-1]=0;
 78     get_inv(a,a_inv,tp,n,l);
 79     FNT(a_inv,(n<<1),l+1,1);
 80     FNT(a_dao,(n<<1),l+1,1);
 81     For(i,0,(n<<1)) a_inv[i]=a_inv[i]*a_dao[i]%p;
 82     FNT(a_inv,(n<<1),l+1,-1);
 83     b[0]=0;
 84     For(i,1,n-1) b[i]=a_inv[i-1]*inv[i]%p;
 85 }
 86 
 87 LL ln_b[N];
 88 void get_exp(LL a[],LL b[],int n,int l) {
 89     if(n==1) { b[0]=1; return; }
 90     get_exp(a,b,(n>>1),l-1);
 91     For(i,0,(n<<1)) ln_b[i]=0;
 92     get_ln(b,ln_b,n,l); ////
 93     For(i,0,n-1) ln_b[i]=((LL)(i==0)-ln_b[i]+a[i]+p)%p;
 94     FNT(b,(n<<1),l+1,1);
 95     FNT(ln_b,(n<<1),l+1,1);
 96     For(i,0,(n<<1)) ln_b[i]=ln_b[i]*b[i]%p;
 97     FNT(ln_b,(n<<1),l+1,-1);
 98     For(i,0,n-1) b[i]=ln_b[i],b[i+n]=0; 
 99 }
100 
101 int main() {
102 #ifdef ANS
103     freopen(".in","r",stdin);
104     freopen(".out","w",stdout);
105 #endif
106     read(n);
107     For(i,0,n-1) read(a[i]);
108     int nn,l;
109     for(nn=1,l=0;nn<=n;nn<<=1) l++;
110     inv[0]=inv[1]=1;
111     For(i,2,n) inv[i]=mo(p-p/i*inv[p%i]%p);
112     get_exp(a,b,nn,l);
113     For(i,0,n-1) printf("%lld ",b[i]);
114     Formylove;
115 }

多项式指数函数(exp)

 1 LL h_sqr[N],a_t[N],b_t[N];
 2 void get_sqr(LL a[],LL b[],int n,int l) {
 3     if(n==1) {
 4         b[0]=sqrt(a[0]);
 5         return;
 6     }
 7     get_sqr(a,b,n>>1,l-1);
 8     For(i,0,n<<1) b_t[i]=0;
 9     For(i,0,n-1) a_t[i]=a[i],a_t[i+n]=0;
10     get_inv(b,b_t,n,l);
11     FFT(b,n<<1,l+1,1);
12     FFT(a_t,n<<1,l+1,1);
13     FFT(b_t,n<<1,l+1,1);
14     For(i,0,n<<1) b[i]=(b[i]+a_t[i]*b_t[i]%p)%p*inv_2%p;
15     FFT(b,n<<1,l+1,-1);
16     For(i,0,n-1) b[n+i]=0;
17 }

多项式开方(牛顿迭代)

  1 //Achen
  2 #include<algorithm>
  3 #include<iostream>
  4 #include<cstring>
  5 #include<cstdlib>
  6 #include<vector>
  7 #include<cstdio>
  8 #include<queue>
  9 #include<cmath>
 10 #include<set>
 11 #include<map>
 12 #define Formylove return 0
 13 #define For(i,a,b) for(int i=(a);i<=(b);i++)
 14 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
 15 const int N=262155,p=998244353,g=3,gi=332748118; 
 16 typedef long long LL;
 17 typedef double db;
 18 using namespace std;
 19 int n,m;
 20 LL a[N],b[N],D[N],R[N];
 21 
 22 template<typename T>void read(T &x)  {
 23     char ch=getchar(); x=0; T f=1;
 24     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
 25     if(ch==‘-‘) f=-1,ch=getchar();
 26     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
 27 }
 28 
 29 LL ksm(LL a,LL b) {
 30     LL rs=1,bs=a;
 31     while(b) {
 32         if(b&1) rs=rs*bs%p;
 33         bs=bs*bs%p;
 34         b>>=1; 
 35     }    
 36     return rs;
 37 }
 38 
 39 LL mo(LL x) { return x<0?x+p:(x>=p?x-p:x); }
 40 int rev[N];
 41 void FNT(LL a[],int n,int l,int f) {
 42     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
 43     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
 44     for(int i=1;i<n;i<<=1) {
 45         LL wi=ksm((f==1?g:gi),(p-1)/(i<<1));
 46         for(int j=0,pp=(i<<1);j<n;j+=pp) {
 47             LL w=1;
 48             for(int k=0;k<i;k++,w=w*wi%p) {
 49                 LL x=a[j+k],y=w*a[j+k+i]%p;
 50                 a[j+k]=mo(x+y); a[j+k+i]=mo(x-y);
 51             }
 52         } 
 53     }
 54     if(f==-1) {
 55         LL inv_n=ksm(n,p-2);
 56         For(i,0,n) a[i]=a[i]*inv_n%p;
 57     }
 58 }
 59 
 60 LL tp[N];
 61 void get_inv(LL a[],LL b[],int n,int l) {
 62     if(n==1) { 
 63         b[0]=ksm(a[0],p-2); 
 64         return; 
 65     }
 66     get_inv(a,b,(n>>1),l-1);
 67     For(i,0,n-1) tp[i]=a[i],tp[i+n]=0;
 68     FNT(tp,(n<<1),l+1,1);
 69     FNT(b,(n<<1),l+1,1);
 70     For(i,0,(n<<1)) tp[i]=mo(2LL-tp[i]*b[i]%p)*b[i]%p;
 71     FNT(tp,(n<<1),l+1,-1);
 72     For(i,0,n-1) b[i]=tp[i],b[i+n]=0;
 73 }
 74 
 75 void get_R(LL a[],LL a_R[],int n) { For(i,0,n) a_R[i]=a[n-i]; }
 76 
 77 void mul(LL a[],LL b[],LL c[],int n,int m) {
 78     static LL A_l[N],B_l[N];
 79     int nn,l;
 80     for(nn=1,l=0;nn<=n+m;nn<<=1) l++;
 81     For(i,0,n) A_l[i]=a[i]; For(i,n+1,nn) A_l[i]=0;
 82     For(i,0,m) B_l[i]=b[i]; For(i,m+1,nn) B_l[i]=0;
 83     FNT(A_l,nn,l,1);
 84     FNT(B_l,nn,l,1);
 85     For(i,0,nn) (A_l[i]*=B_l[i])%=p;
 86     FNT(A_l,nn,l,-1);
 87     For(i,0,n+m) c[i]=A_l[i];
 88 } 
 89 
 90 void div(LL a[],int n,LL b[],int m,LL D[],LL R[]) {
 91     static LL A_R[N],B_R[N],D_R[N];
 92     int nn,l;
 93     for(nn=1,l=0;nn<=n-m+1;nn<<=1) l++;
 94     For(i,0,(nn<<1)) D_R[i]=0;
 95     get_R(a,A_R,n); 
 96     get_R(b,B_R,m);
 97     get_inv(B_R,D_R,nn,l);
 98     mul(A_R,D_R,D_R,n,n-m);
 99     get_R(D_R,D,n-m);
100     for(nn=1,l=0;nn<=(n<<1);nn<<=1) l++;
101     mul(b,D,R,m,n-m);
102     For(i,0,nn) R[i]=mo(a[i]-R[i]);
103 }
104 
105 //#define ANS
106 int main() {
107 #ifdef ANS
108     freopen("1.in","r",stdin);
109     //freopen(".out","w",stdout);
110 #endif
111     read(n); read(m);
112     For(i,0,n) read(a[i]);
113     For(i,0,m) read(b[i]);
114     div(a,n,b,m,D,R);
115     For(i,0,n-m) printf("%lld ",D[i]); puts("");
116     For(i,0,m-1) printf("%lld ",R[i]); puts("");
117     Formylove;
118 }

多项式除法/取模

一个通往之前写的题的传送门

1.bzoj3527: [Zjoi2014]力

$F_j= \sum_{i<j} q_i/(i-j)^2-\sum_{i>j}q_i/(i-j)^2$

设$E1_j=\sum_{i<j} q_i/(j-i)^2$

$E2_j=\sum_{i>j}q_i/(i-j)^2$

$B={1/i^2},B_0=0$

$F=E1-E2$

$E1=q*B$

$E2_j= \sum_{i=j}^{n-1}q_i/(i-j)^2$

设$qR_{i}=q_{n-i-1}$

$E2_j= \sum_{i=j}^{n-1}qR_{n-i-1}/(i-j)^2$

$=\sum_{j+k=n-j-1}qR_j*B_k$

$=qR*B(n-j-1)$

 1 // luogu-judger-enable-o2
 2 //Achen
 3 #include<algorithm>
 4 #include<iostream>
 5 #include<cstring>
 6 #include<cstdlib>
 7 #include<vector>
 8 #include<cstdio>
 9 #include<queue>
10 #include<cmath>
11 #include<set>
12 #include<map>
13 #define pi acos(-1)
14 #define Formylove return 0
15 #define For(i,a,b) for(int i=(a);i<=(b);i++)
16 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
17 const int N=2100007;
18 typedef long long LL;
19 typedef double db;
20 using namespace std;
21 int n,nn,m,l,rev[N];
22 
23 template<typename T>void read(T &x)  {
24     char ch=getchar(); x=0; T f=1;
25     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
26     if(ch==‘-‘) f=-1,ch=getchar();
27     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
28 }
29 
30 struct E {
31     db a,b;
32     E(db a=0,db b=0):a(a),b(b){}
33 }p[N],p_R[N],b[N];
34 E operator + (const E&A,const E&B) { return E(A.a+B.a,A.b+B.b); }
35 E operator - (const E&A,const E&B) { return E(A.a-B.a,A.b-B.b); } 
36 E operator * (const E&A,const E&B) { return E(A.a*B.a-A.b*B.b,A.a*B.b+A.b*B.a); } 
37 E operator / (const E&A,const int&B) { return E(A.a/B,A.b/B); } 
38 
39 void FFT(E a[],int n,int f) {
40     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
41     for(int i=1;i<n;i<<=1) {
42         E wi(cos(pi/i),f*sin(pi/i));
43         for(int j=0,pp=(i<<1);j<n;j+=pp) {
44             E w(1,0);
45             for(int k=0;k<i;k++,w=w*wi) {
46                 E x=a[j+k],y=w*a[j+k+i];
47                 a[j+k]=x+y; a[j+k+i]=x-y;
48             }
49         }
50     }
51     if(f==-1) For(i,0,n) a[i]=a[i]/n;
52 } 
53 
54 int main() {
55 #ifdef ANS
56     freopen(".in","r",stdin);
57     freopen(".out","w",stdout);
58 #endif
59     read(n);
60     For(i,0,n-1) { scanf("%lf",&p[i].a); p_R[n-i-1]=p[i]; }
61     For(i,1,n-1) b[i].a=1.0/i/i;
62     for(nn=1;nn<=n+n+1;nn<<=1) l++;
63     For(i,0,nn) rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
64     FFT(p,nn,1); FFT(p_R,nn,1); FFT(b,nn,1);
65     For(i,0,nn) p[i]=p[i]*b[i],p_R[i]=p_R[i]*b[i];
66     FFT(p,nn,-1); FFT(p_R,nn,-1);
67     For(i,0,n-1) printf("%lf\n",p[i].a-p_R[n-i-1].a);
68     Formylove;
69 }

View Code

2.bzoj3160: 万径人踪灭

先把序列中间加上间隔符。

把a的位置设为1，其余为0，自己卷自己得到第2*i项就是关于i对称的位置的数量。同理处理b。然后每个位置计算答案累加。

不允许连续，再跑一遍马拉车减去答案即可。

  1 //Achen
  2 #include<algorithm>
  3 #include<iostream>
  4 #include<cstring>
  5 #include<cstdlib>
  6 #include<vector>
  7 #include<cstdio>
  8 #include<queue>
  9 #include<cmath>
 10 #include<set>
 11 #include<map>
 12 #define pi acos(-1)
 13 #define Formylove return 0
 14 #define For(i,a,b) for(int i=(a);i<=(b);i++)
 15 #define Rep(i,a,b) for(int i=(a);i>=(b);i--)
 16 const int N=1e6+7,mod=1e9+7;
 17 typedef long long LL;
 18 typedef double db;
 19 using namespace std;
 20 int len,n,l,cnt,rev[N];
 21 char s[N],ss[N*2];
 22 LL ans;
 23 
 24 template<typename T>void read(T &x)  {
 25     char ch=getchar(); x=0; T f=1;
 26     while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
 27     if(ch==‘-‘) f=-1,ch=getchar();
 28     for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
 29 }
 30 
 31 struct E {
 32     db a,b;
 33     E(db a=0,db b=0):a(a),b(b){}
 34 }a[N],b[N];
 35 E operator + (const E&A,const E&B) { return E(A.a+B.a,A.b+B.b); }
 36 E operator - (const E&A,const E&B) { return E(A.a-B.a,A.b-B.b); } 
 37 E operator * (const E&A,const E&B) { return E(A.a*B.a-A.b*B.b,A.a*B.b+A.b*B.a); } 
 38 E operator / (const E&A,const int&B) { return E(A.a/B,A.b/B); } 
 39 
 40 void FFT(E a[],int n,int f) {
 41     For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
 42     for(int i=1;i<n;i<<=1) {
 43         E wi(cos(pi/i),f*sin(pi/i));
 44         for(int j=0,pp=(i<<1);j<n;j+=pp) {
 45             E w(1,0);
 46             for(int k=0;k<i;k++,w=w*wi) {
 47                 E x=a[j+k],y=w*a[j+k+i];
 48                 a[j+k]=x+y; a[j+k+i]=x-y;
 49             }
 50         }
 51     }
 52     if(f==-1) For(i,0,n) a[i]=a[i]/n;
 53 } 
 54 
 55 LL ksm(LL a,LL b) {
 56     LL rs=1,bs=a;
 57     while(b) {
 58         if(b&1) rs=rs*bs%mod;
 59         bs=bs*bs%mod;
 60         b>>=1;
 61     }
 62     return rs;
 63 }
 64 
 65 int rad[N<<1];
 66 void manacher() {
 67     LL tot=len;
 68     for(int i=1,k=0,j;i<cnt;) {
 69         while(ss[i-k-1]==ss[i+k+1]) k++;
 70         rad[i]=k;
 71         for(j=1;j<=k&&rad[i-j]!=rad[i]-j;j++)
 72             rad[i+j]=min(rad[i]-j,rad[i-j]);
 73         i+=j;
 74         k=max(0,k-j);
 75     }
 76     For(i,0,cnt-1) tot=(tot+rad[i]/2)%mod;
 77     ans=(ans-tot+mod)%mod;
 78 }
 79 
 80 int main() {
 81 #ifdef ANS
 82     freopen(".in","r",stdin);
 83     freopen(".out","w",stdout);
 84 #endif
 85     scanf("%s",s);
 86     len=strlen(s);
 87     cnt=0;
 88     ss[cnt++]=‘&‘;
 89     ss[cnt++]=‘#‘;
 90     For(i,0,len-1) {
 91         ss[cnt++]=s[i];
 92         ss[cnt++]=‘#‘;
 93     } 
 94     ss[cnt++]=‘*‘;
 95     For(i,0,cnt-1) 
 96         if(ss[i]==‘a‘) a[i].a=1;
 97         else if(ss[i]==‘b‘) b[i].a=1;
 98     for(n=1;n<=cnt*2;n<<=1) l++;
 99     For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
100     FFT(a,n,1); FFT(b,n,1);
101     For(i,0,n) a[i]=a[i]*a[i],b[i]=b[i]*b[i];
102     FFT(a,n,-1); FFT(b,n,-1);
103     For(i,0,cnt-1) {
104         LL c; 
105         if(ss[i]==‘a‘) c=((int)(a[2*i]+0.5).a+1)/2+((int)(b[2*i].a+0.5))/2;
106         else if(ss[i]==‘b‘) c=((int)(a[2*i].a+0.5))/2+((int)(b[2*i].a+0.5)+1)/2;
107         else c=((int)(a[2*i].a+0.5))/2+((int)(b[2*i].a+0.5))/2;
108         if(c>0) ans=(ans+ksm(2,c)-1+mod)%mod;
109     }
110     manacher();
111     printf("%lld\n",ans);
112     Formylove;
113 }
114 /*
115 abaabaa
116 
117 aaabbbaaa
118 
119 aaaaaaaa
120 */

View Code

3.bzoj4259: 残缺的字符串

把*值设为0，定义比较函数f(i,j)=(a[i]-b[j])^2*a[i]*b[j]。

把b取反化成卷积，拆开后FFT即可。

把(len1-1)/2写成len1/2交换的时候多换了一位，调了一晚上。。。然后又数组开大了一会T一会MLE。。。

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<queue>
#include<cmath>
#include<set>
#include<map>
#define pi acos(-1)
#define Formylove return 0
#define For(i,a,b) for(int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(int i=(a);i>=(b);i--)
const int N=1048576;
typedef long long LL;
typedef double db;
using namespace std;
int len1,len2,n,l,rev[N],q[N],a[N],b[N];
char s1[N],s2[N];

template<typename T>void read(T &x)  {
    char ch=getchar(); x=0; T f=1;
    while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
    if(ch==‘-‘) f=-1,ch=getchar();
    for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
}

struct E {
    db a,b;
    E(db a=0,db b=0):a(a),b(b){}
}A[N],B[N],t[N];
E operator + (const E&A,const E&B) { return E(A.a+B.a,A.b+B.b); }
E operator - (const E&A,const E&B) { return E(A.a-B.a,A.b-B.b); } 
E operator * (const E&A,const E&B) { return E(A.a*B.a-A.b*B.b,A.a*B.b+A.b*B.a); } 
E operator / (const E&A,const int&B) { return E(A.a/B,A.b/B); } 

void FFT(E a[],int n,int f) {
    For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1) {
        E wi(cos(pi/i),f*sin(pi/i));
        for(int j=0,pp=(i<<1);j<n;j+=pp) {
            E w(1,0);
            for(int k=0;k<i;k++,w=w*wi) {
                E x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y; a[j+k+i]=x-y;
            }
        }
    }
    if(f==-1) For(i,0,n) a[i]=a[i]/n;
} 

int main() {
    read(len1); read(len2);
    scanf("%s%s",s1,s2);
    For(i,0,(len1-1)/2) swap(s1[i],s1[len1-i-1]);
    For(i,0,len1-1) if(s1[i]!=‘*‘) b[i]=s1[i]-‘a‘+1;
    For(i,0,len2-1) if(s2[i]!=‘*‘) a[i]=s2[i]-‘a‘+1;
    for(n=1;n<=len1+len2;n<<=1) l++;
    For(i,1,n) rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));
    
    For(i,0,len1-1) B[i].a=b[i]*b[i]*b[i];
    For(i,0,len2-1) A[i].a=a[i];
    FFT(A,n,1),FFT(B,n,1); 
    For(i,0,n) t[i]=A[i]*B[i];
    
    For(i,0,len1-1) B[i]=E(b[i],0); For(i,len1,n) B[i]=E(0,0);
    For(i,0,len2-1) A[i]=E(a[i]*a[i]*a[i],0); For(i,len2,n) A[i]=E(0,0);
    FFT(A,n,1),FFT(B,n,1); For(i,0,n) t[i]=t[i]+A[i]*B[i];
    
    For(i,0,len1-1) B[i]=E(b[i]*b[i],0); For(i,len1,n) B[i]=E(0,0);
    For(i,0,len2-1) A[i]=E(a[i]*a[i],0); For(i,len2,n) A[i]=E(0,0);
    FFT(A,n,1),FFT(B,n,1); For(i,0,n) t[i]=t[i]-A[i]*B[i]*E(2,0);
    
    /*For(i,0,n) { 
        t[i]=A[2][i]*B[0][i]+A[0][i]*B[2][i]-A[1][i]*B[1][i]-A[1][i]*B[1][i];
    }*/
    FFT(t,n,-1);
    int ans=0;
    For(i,len1-1,len2-1) 
        if((int)(t[i].a+0.5)==0) ans++,q[++q[0]]=i-len1+2;
    printf("%d\n",ans);
    if(q[0]>1) For(i,1,q[0]-1) printf("%d ",q[i]);
    if(q[0]) printf("%d",q[ans]);
    Formylove;
}
/*
4 30
d*ad
*cbacbeeae*b**debabcecdb*aaacb
*/

View Code

4.hdu4609 3-idiots

FFT求出和为i的木棍对的数量，枚举三角形最长的一边，把和比它大的木棍对的数量减去不合法的（两根都比它小为合法），再加上等边和等腰即可。

因为清0和开LL又wa了好久。。。

//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<queue>
#include<cmath>
#include<set>
#include<map>
#define Formylove return 0
#define For(i,a,b) for(int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(int i=(a);i>=(b);i--)
typedef long long LL;
typedef double db;
const int N=262149;
const db pi=acos(-1);
using namespace std;
int T,n,m,l,rev[N],L[N];
LL ans,suml[N],cnt[N];

template<typename T>void read(T &x)  {
    char ch=getchar(); x=0; T f=1;
    while(ch!=‘-‘&&(ch<‘0‘||ch>‘9‘)) ch=getchar();
    if(ch==‘-‘) f=-1,ch=getchar();
    for(;ch>=‘0‘&&ch<=‘9‘;ch=getchar()) x=x*10+ch-‘0‘; x*=f;
}

struct E {
    db a,b;
    E(db a=0,db b=0):a(a),b(b){}
}a[N],b[N];
E operator +(const E&A,const E&B) { return E(A.a+B.a,A.b+B.b); }
E operator -(const E&A,const E&B) { return E(A.a-B.a,A.b-B.b); }
E operator *(const E&A,const E&B) { return E(A.a*B.a-A.b*B.b,A.a*B.b+A.b*B.a); }
E operator /(const E&A,const db&B) { return E(A.a/B,A.b/B); }

void FFT(E a[],int n,int f) {
    For(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1) {
        E wi(cos(pi/i),f*sin(pi/i));
        for(int j=0,pp=(i<<1);j<n;j+=pp) {
            E w(1,0);
            for(int k=0;k<i;k++,w=w*wi) {
                E x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y; a[j+k+i]=x-y;
            } 
        }
    }
    if(f==-1) For(i,0,n) a[i]=a[i]/n; 
}

LL C_2(int n) { return n<2?0:(LL)n*(n-1)/2;}
LL C_3(int n) { return n<3?0:(LL)n*(n-1)*(n-2)/6; }

//#define ANS
int main() {
#ifdef ANS
    freopen("std.in","r",stdin);
    //freopen(".out","w",stdout);
#endif
    read(T);
    while(T--) {
        read(m); int mxlen=0;
        For(i,1,m) { read(L[i]); mxlen=max(mxlen,L[i]); }
        for(n=1,l=0;n<=mxlen+mxlen;n<<=1) l++;
        For(i,1,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1)); 
        For(i,0,mxlen*2) cnt[i]=0;
        For(i,0,n) a[i]=E(0,0);
        For(i,1,m) a[L[i]].a++,cnt[L[i]]++;
        FFT(a,n,1);
        For(i,0,n) a[i]=a[i]*a[i];
        FFT(a,n,-1);
        For(i,0,n) a[i].a=(LL)(a[i].a+0.5);
        For(i,1,mxlen) a[i*2].a-=cnt[i];
        For(i,0,n) a[i].a/=2;
        For(i,1,mxlen) suml[i]=suml[i-1]+cnt[i];
        LL sum=0; ans=0;
        Rep(i,mxlen*2,1) {
            if(cnt[i]) {
                LL tp=suml[mxlen]-suml[i];
                ans+=(sum-tp*suml[i-1]-C_2(tp)-C_2(cnt[i])-(LL)cnt[i]*(m-cnt[i]))*cnt[i]+C_3(cnt[i])+C_2(cnt[i])*suml[i-1];
            }
            sum+=a[i].a;
        }
        db pt=(db)ans/C_3(m);
        printf("%.7lf\n",pt); 
    }
    Formylove;
}
/*
1
7
2 3 2 2 3 3 4
*/

View Code

关于万恶的多项式

标签：就是 -o 取反 ++ 传送门 aac 除法 sed include

原文地址：https://www.cnblogs.com/Achenchen/p/9456439.html

踩

(0)

评论一句话评论（0）

分享档案

更多>

2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)

周排行