标签：
题意：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] ) xi = 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 }
标签：
原文地址：http://www.cnblogs.com/dirge/p/5933103.html
