Partial Sums V2
题目连接：
https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1172
Description
给出一个数组A，经过一次处理，生成一个 数组S，数组S中的每个值相当于数组A的累加，比如：A = {1 3 5 6} => S = {1 4 9 15}。如果对生成的数组S再进行一次累加操作，{1 4 9 15} => {1 5 14 29}，现在给出数组A，问进行K次操作后的结果。(输出结果 Mod 10^9 + 7)
Input
第1行，2个数N和K，中间用空格分隔，N表示数组的长度，K表示处理的次数（2 <= n <= 50000, 0 <= k <= 10^9, 0 <= a[i] <= 10^9)
Output
共N行，每行一个数，对应经过K次处理后的结果。输出结果 Mod 10^9 + 7。
Sample Input
4 2 1 3 5 6
Sample Output
1 5 14 29
题意：
题解：
考虑每个位置的贡献，所求即$ans[n]=\sum_{i+j=n}^{}C(k-1+i,i)\cdot a[j]$ ，其中 $C(n,i) $表示从 n 个物品里选出 i 个物品的方案数。
上述卷积形式可以利用FFT或NTT加速运算，难点在于答案模 109+7，设其为 P ，这意味这参与卷积的每个元素值为不大于 (P−1) 的非负整数，卷积后的值为不大于 (n+1)⋅(P−1)2=O(nP2) 的值，这个上界可以达到 5⋅1022 。
一种想法是利用3种模意义下的NTT结果进行合并，得到这个可能达到 5⋅1022 的数，再取模，一共需要9次NTT，涉及高精度。
还有一种想法就是降低卷积前的元素值上界，令 M 为满足 P≤M2 的最小正整数，将多项式 f 拆成 f1+f2⋅M 的形式，那么每个元素的值上界是 O(M)=O($\sqrt{p}$) 的，卷积后的值上界是 O(nP) 的，用long double来做FFT差不多可以接受，而卷积可以写为 f∗g=(f1+f2⋅M)∗(g1+g2⋅M)=(f1∗f2)+(f1∗g2+f2∗g1)⋅M+(f2∗g2)⋅M2 ，一共需要7次FFT，不涉及高精度。
很邪，51nod上c++会wa，c++11会过..
以前只在大牛的博客上看到过第一种做法，但没有实用性，第二种方法很强。
昨天刚好看到，拉出来封个版
代码
1 //#include <bits/stdc++.h>
2 #include <stdio.h>
3 #include <iostream>
4 #include <string.h>
5 #include <math.h>
6 #include <stdlib.h>
7 #include <limits.h>
8 #include <algorithm>
9 #include <queue>
10 #include <vector>
11 #include <set>
12 #include <map>
13 #include <stack>
14 #include <bitset>
15 #include <string>
16 #include <time.h>
17 using namespace std;
18 long double esp=1e-11;
19 //#pragma comment(linker, "/STACK:1024000000,1024000000")
20 #define fi first
21 #define se second
22 #define all(a) (a).begin(),(a).end()
23 #define cle(a) while(!a.empty())a.pop()
24 #define mem(p,c) memset(p,c,sizeof(p))
25 #define mp(A, B) make_pair(A, B)
26 #define pb push_back
27 #define lson l , m , rt << 1
28 #define rson m + 1 , r , rt << 1 | 1
29 typedef long long int LL;
30 const long double PI = acos((long double)-1);
31 const LL INF=0x3f3f3f3fll;
32 const int MOD =1000000007ll;
33 const int maxn=140100;
34 struct complex
35 {
36 long double r,i;
37 complex(long double _r = 0.0,long double _i = 0.0){r = _r; i = _i;}
38 complex operator +(const complex &b){return complex(r+b.r,i+b.i);}
39 complex operator -(const complex &b){return complex(r-b.r,i-b.i);}
40 complex operator *(const complex &b){return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
41 };
42 void change(complex y[],int len)
43 {
44 int i,j,k;
45 for(i = 1, j = len/2;i < len-1; i++)
46 {
47 if(i < j)swap(y[i],y[j]);
48 k = len/2;
49 while( j >= k)
50 {
51 j -= k;
52 k /= 2;
53 }
54 if(j < k) j += k;
55 }
56 }
57 void FFT(complex y[],int len,int on) //len=2^k
58 {
59 change(y,len);
60 for(int h = 2; h <= len; h <<= 1)
61 {
62 complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
63 for(int j = 0;j < len;j+=h)
64 {
65 complex w(1,0);
66 for(int k = j;k < j+h/2;k++)
67 {
68 complex u = y[k];
69 complex t = w*y[k+h/2];
70 y[k] = u+t;
71 y[k+h/2] = u-t;
72 w = w*wn;
73 }
74 }
75 }
76 if(on == -1)
77 for(int i = 0;i < len;i++)
78 y[i].r /= len;
79 }
80 int callen(int len1,int len2){
81 int len=1;
82 while(len < (len1<<1) || len < (len2<<1))len<<=1;
83 return len;
84 }
85 LL fftans[maxn];
86 complex tf[5][maxn];
87 void fft(LL* y1,int len1,LL* y2,int len2,LL mod){
88 LL p=sqrt(mod)+10;
89 int len=callen(len1,len2);
90 for(int i = 0;i < len1;i++)
91 {
92 tf[0][i] = complex(y1[i]%p,0);
93 tf[1][i] = complex(y1[i]/p,0);
94 }
95 for(int i = len1;i < len;i++)
96 {
97 tf[0][i] = complex(0,0);
98 tf[1][i] = complex(0,0);
99 }
100 for(int i = 0;i < len2;i++)
101 {
102 tf[2][i] = complex(y2[i]%p,0);
103 tf[3][i] = complex(y2[i]/p,0);
104 }
105 for(int i = len2;i < len;i++)
106 {
107 tf[2][i] = complex(0,0);
108 tf[3][i] = complex(0,0);
109 }
110 for(int i = 0;i < 4;i++)
111 FFT(tf[i],len,1);
112 memset(fftans,0,(len1+len2-1)*sizeof(LL));
113
114 for(int i = 0;i < len;i++)
115 tf[4][i] = tf[0][i]*tf[2][i];
116 FFT(tf[4],len,-1);
117 for(int i = 0;i < len1+len2-1;i++)
118 fftans[i]=(fftans[i]+(LL)(tf[4][i].r+0.5))%mod;
119
120 for(int i = 0;i < len;i++)
121 tf[4][i] = tf[0][i]*tf[3][i]+tf[1][i]*tf[2][i];
122 FFT(tf[4],len,-1);
123 for(int i = 0;i < len1+len2-1;i++)
124 fftans[i]=(fftans[i]+((LL)(tf[4][i].r+0.5))%mod*p)%mod;
125
126 for(int i = 0;i < len;i++)
127 tf[4][i] = tf[1][i]*tf[3][i];
128 FFT(tf[4],len,-1);
129 for(int i = 0;i < len1+len2-1;i++)
130 fftans[i]=(fftans[i]+((LL)(tf[4][i].r+0.5))%mod*p%mod*p)%mod;
131 }
132 LL inv[50001];
133 LL a[maxn],b[maxn];
134
135 int main()
136 {
137 //freopen("in.txt", "r", stdin);
138 //freopen("inlay.in", "r", stdin);
139 //freopen("out.txt", "w", stdout); //%I64d
140 //vector<int>::iterator iter;
141 //for(int x=1;x<=n;x++)
142 //for(int y=1;y<=n;y++)
143 //scanf("%d",&a);
144 //printf("%d\n",ans);
145 int n;
146 LL k;
147 scanf("%d%lld",&n,&k);
148 inv[1]=1;
149 for(int x=2;x<n;x++)inv[x]=(MOD-MOD/x)*inv[MOD%x]%MOD;
150 for(int x=0;x<n;x++)
151 scanf("%lld",&a[x]);
152 LL tem=1;
153 b[0]=1;
154 for(int x=1;x<n;x++)
155 {
156 tem=tem*(x+k-1)%MOD*inv[x]%MOD;
157 b[x]=tem;
158 }
159 fft(a,n,b,n,MOD);
160 for(int x=0;x<n;x++)
161 printf("%lld\n",(fftans[x]+MOD)%MOD);
162 return 0;
163 }
2016-07-30
时间： 2024-10-06 11:21:02