HDU5730 FFT+CDQ分治

题意：dp[n] = ∑ ( dp[n-i]*a[i] )+a[n], ( 1 <= i < n)

cdq分治。

计算出dp[l ~ mid]后，dp[l ~ mid]与a[1 ~ r-l]做卷积运算。

 1 #include <bits/stdc++.h>
 2 using namespace std;
 3 const int N = 4e5+5;
 4 const int mod = 313;
 5 const double pi = acos(-1.0);
 6 struct comp{
 7     double r,i;comp(double _r=0,double _i=0){r=_r;i=_i;}
 8     comp operator+(const comp x){return comp(r+x.r,i+x.i);}
 9     comp operator-(const comp x){return comp(r-x.r,i-x.i);}
10     comp operator*(const comp x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
11 }x[N], y[N];
12 void FFT(comp a[],int n,int t){
13     for(int i=1,j=0;i<n-1;i++){
14         for(int s=n;j^=s>>=1,~j&s;);
15         if(i<j)swap(a[i],a[j]);
16     }
17     for(int d=0;(1<<d)<n;d++){
18         int m=1<<d,m2=m<<1;
19         double o=pi/m*t;comp _w(cos(o),sin(o));
20         for(int i=0;i<n;i+=m2){
21             comp w(1,0);
22             for(int j=0;j<m;j++){
23                 comp &A=a[i+j+m],&B=a[i+j],t=w*A;
24                 A=B-t;B=B+t;w=w*_w;
25             }
26         }
27     }
28     if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;
29 }
30
31 int a[N], dp[N], n;
32 void cdq(int l,int r){
33     if(l == r){
34         dp[l] = (dp[l]+a[l])%mod;
35         return;
36     }
37     int mid = l+r >> 1;
38     cdq(l, mid);
39     int len = 1;
40     while(len <= (r-l+1)) len <<= 1;
41     for(int i = 0; i < len; i++)
42         x[i] = y[i] = comp(0, 0);
43     for(int i = l; i <= mid; i++)
44         x[i-l] = comp(dp[i], 0);
45     for(int i = 1; l+i <= r; i++)
46         y[i-1] = comp(a[i], 0);
47     FFT(x, len, 1); FFT(y, len, 1);
48     for(int i = 0; i < len; i++)
49         x[i] = x[i]*y[i];
50     FFT(x, len, -1);
51
52     for(int i = mid+1; i <= r; i++){
53         dp[i] += x[i-l-1].r+0.5;
54         dp[i] %= mod;
55     }
56     cdq(mid+1,r);
57 }
58
59 int main(){
60     while(scanf("%d", &n), n){
61         for(int i = 1; i <= n; i++){
62             scanf("%d", a+i);
63             a[i] %= mod;
64             dp[i] = 0;
65         }
66         cdq(1, n);
67         printf("%d\n", dp[n]);
68     }
69     return 0;
70 }

补：

因为做多项式逆运算在分治FFT之前，故做此题时首先有了一个多项式求逆的方法。

观察 dp[n] = ∑ ( dp[n-i]*a[i] )+a[n], ( 1 <= i < n)

令a[0] = -1, 则 ∑ ( dp[n-i]*a[i] )+a[n] = 0, ( 0 <= i < n) 对任意n > 0成立

令dp[0] = 1

则dp序列与a序列卷积 = dp[0]*a[0]  +  ∑ ( dp[i-j]*a[j] ) xⁱ = dp[0]*a[0] = -1,

已知a序列，则dp序列为a的逆多项式序列*(-1)。

因为模313，最终各个系数 <= n(m-1)*(m-1) = 9734400000,  故选了一个大素数206158430209

写了一个多项式求逆的代码，不过不知什么原因WA了。

 1 #include<cstdio>
 2 #include<iostream>
 3 typedef long long ll;
 4 using namespace std;
 5 const int N = 262144, K = 17;
 6 ll n, m, i, k;
 7 ll a[N+10], b[N+10], tmp[N+10], tmp2[N+10];
 8 ll P = 206158430209LL, G = 22, g[K+1], ng[K+10], inv[N+10], inv2;
 9 //ll P = 998244353, G = 3, g[K+1], ng[K+10], inv[N+10], inv2;
10 ////////////
11 ll mul(ll a, ll b){
12     ll ans = 0;
13     while(b){
14         if(b&1) ans += a;
15         if(ans >= P) ans -= P;
16         a = a+a;
17         if(a >= P) a -= P;
18         b >>= 1;
19     }
20     return ans;
21 }
22 ll pow(ll a, ll n){
23     ll ans = 1;
24     while(n){
25         if(n&1)
26             ans = mul(ans, a);
27         a = mul(a, a);
28         n >>= 1;
29     }
30     return ans;
31 }
32 ////////////
33 void NTT(ll*a,int n,int t){
34     for(int i=1,j=0;i<n-1;i++){
35         for(int s=n;j^=s>>=1,~j&s;);
36         if(i<j)swap(a[i], a[j]);
37     }
38     for(int d=0;(1<<d)<n;d++){
39         int m=1<<d,m2=m<<1;
40         ll _w=t==1?g[d]:ng[d];///////////////////////////////////
41
42         for(int i=0;i<n;i+=m2)for(ll w=1,j=0;j<m;j++){
43             ll&A=a[i+j+m],&B=a[i+j],t=mul(w, A);////////////(ll)w*A%P;
44             A=B-t;if(A<0)A+=P;
45             B=B+t;if(B>=P)B-=P;
46             w=mul(w, _w);/////////(ll)w*_w%P;/////////////////////
47         }
48     }
49     //if(t==-1)for(int i=0,j=inv[n];i<n;i++)a[i]=mul(a[i], j);///////(ll)a[i]*j%P;
50     if(t==-1){
51         ll j = pow(n, P-2);
52         for(int i=0;i<n;i++)a[i]=mul(a[i], j);
53     }
54 }
55 //给定a,求a的逆元b
56 void getinv(ll*a,ll*b,int n){
57     if(n==1){b[0]=pow(a[0],P-2);return;}
58     getinv(a,b,n>>1);
59     int k=n<<1,i;
60     for(i=0;i<n;i++)tmp[i]=a[i];
61     for(i=n;i<k;i++)tmp[i]=b[i]=0;
62     NTT(tmp,k,1),NTT(b,k,1);
63
64     for(i=0;i<k;i++){
65         ll xx = b[i];
66         b[i] = mul(xx, (P+2-mul(xx, tmp[i]))%P);
67     //b[i]=(ll)b[i]*(2-(ll)tmp[i]*b[i]%P)%P;
68     //if(b[i]<0)b[i]+=P;
69     }
70     NTT(b,k,-1);
71     for(i=n;i<k;i++)b[i]=0;
72 }
73 int main(){
74     for(g[K]=pow(G,(P-1)/N),ng[K]=pow(g[K],P-2),i=K-1;~i;i--){
75         //g[i]=(ll)g[i+1]*g[i+1]%P,ng[i]=(ll)ng[i+1]*ng[i+1]%P;
76         g[i] = mul(g[i+1], g[i+1]);
77         ng[i] = mul(ng[i+1], ng[i+1]);
78     }
79     //for(inv[1]=1,i=2;i<=N;i++)inv[i]=(ll)(P-inv[P%i])*(P/i)%P;inv2=inv[2];
80     while(scanf("%lld", &n), n){
81         int len = 1;
82         while(len <= n) len <<= 1;
83         //len <= 131072
84         a[0] = (P-1)%313;
85         for(i = 1; i <= n; i++){
86             scanf("%lld", a+i);
87             a[i] %= 313;
88         }
89         for(i = n+1; i < len; i++) a[i] = 0;
90
91         getinv(a, b, len);
92         b[n] = b[n]%313*(P-1LL)%313;
93 //        for(i = 0; i < len; i++)
94 //            b[i] = mul(b[i], P-1);
95 //              b[i] = b[i]*(P-1LL)%P;
96         printf("%lld\n", b[n]);
97     }
98     return 0;
99 }