一个循环卷积的问题

无心人

麻烦你把那个PDF发上来

ssikkiss

或者你在网上搜索Winograd WFTA也能找到相应的信息
在WFTA中,Winograd 为了计算FFT,把大N点FFT转化成小n的循环卷积问题,
而我们要做乘法,就是一个循环卷积问题
过去我们是用FFT来做的,也就是把循环卷积问题转化为FFT的计算
现在,如果我们对Winograd的那个"小n点快速卷及算法"弄懂了,就不必绕圈子了!

无心人

那个和FNT关系不大
那个算法是把高度复合数N点变换
分解成小n变换，高效的2-12点的变换

ssikkiss

发了3次!发不了!太大了,9MB,

ssikkiss

那个和FNT关系不大
那个算法是把高度复合数N点变换
分解成小n变换，高效的2-12点的变换

你说的这个是针对整个WFTA,而我说的是WFTA中Winograd提出的那个 小n点循环卷及算法

不知怎么发不上来文件,我把它用WINRAR分成了24分,还是发不上来
要么,留个email?

无心人

yaojialin@126.com

ssikkiss

这是用Z变换来计算卷及,不是与TAOCP上的一样么?

ssikkiss

发了,你看看收到没有

无心人

谢谢

ssikkiss

此题网上以有人做出。
我抄写了它的c++代码，编译出错，稍更改了下，可以通过，不过请自己加入main函数

1. 
   
2. const int NMax = 1<<22;
   
3. const int SomeMore = 600;
   
4. const int PolyTempMax = 1<<11;
   
5. #define NUMBERTYPE long
   
6. NUMBERTYPE * Global = new NUMBERTYPE[4*NMax + SomeMore];
   
7. long CountAdd;
   
8. long CountMul;
   
9. // Standard circular convolution
   
10. void CircConv(long N, NUMBERTYPE *X, NUMBERTYPE *H, NUMBERTYPE *Y) {
    
11.     long i, j;
    
12.     for (i = 0; i < N ;i++)
    
13.         Y[i] = 0;
    
14.     for (i = 0; i < N ;i++)
    
15.         for (j = 0; j < N ;j++)
    
16.             if (i + j < N)
    
17.                 Y[i + j] += X[j] * H[i];
    
18.             else
    
19.                 Y[i + j - N] += X[j] * H[i];
    
20. }
    
21. void WriteResult(long N, NUMBERTYPE *X) {
    
22.     long i;
    
23.     long t=0;
    
24.     for (i = 0; i < N; i++) {
    
25.         t+=X[i];
    
26.         cout<<t%10;
    
27.         t/=10;
    
28.     }
    
29.     printf("\n");
    
30. }
    
31. void ModuloPlusMinus(long N, long X) {
    
32.     long i;
    
33.     NUMBERTYPE Temp;
    
34.     for (i = 0; i < N; i++) {
    
35.         Temp = Global[X + i];
    
36.         Global[X + i] = Temp - Global[X + i + N];
    
37.         Global[X + i + N] = Temp + Global[X + i + N];
    
38.     }
    
39.     CountAdd += N<<1;
    
40. }
    
41. void MulAddDiv(long N, long X1, long X2) {
    
42.     long i;
    
43.     NUMBERTYPE Temp;
    
44.     for (i = 0; i < N; i++) {
    
45.         Temp = Global[X1 + i];
    
46.         Global[X1 + i] = (Temp + Global[X2 + i])>>1;
    
47.         Global[X1 + i + N] = (-Temp + Global[X2 + i])>>1;
    
48.     }
    
49.     CountAdd += N<<1;
    
50. }
    
51. void PermutePolySub(long N, long K, long X, NUMBERTYPE * Y, long Z) {
    
52.     long i;
    
53.     long adr;
    
54.     long ex;
    
55.     for (i = 0; i < N; i++) {
    
56.         adr = (i + K +  2 * N * N) % N;
    
57.         ex = ((i + K +  2 * N * N) / N) & 1;
    
58.         if (ex == 0) {
    
59.             Global[Z + adr] = Global[X + i] - Y[i];
    
60.         } else {
    
61.             Global[Z + adr] = -Global[X + i] + Y[i];
    
62.         }
    
63.     }
    
64.     CountAdd +=N;
    
65. }
    
66. 
    
67. void PolyAdd(long N, long X, NUMBERTYPE * Y, long Z) {
    
68.     long i;
    
69.     for (i = 0; i < N; i++) {
    
70.         Global[Z + i] = Global[X + i] + Y[i];
    
71.     }
    
72.     CountAdd +=N;
    
73. }
    
74. 
    
75. void PolySubPermute(long N, long K, long X, NUMBERTYPE * Y, long Z) {
    
76.     long i;
    
77.     long adr;
    
78.     long ex;
    
79. 
    
80.     for (i = 0; i < N; i++) {
    
81.         adr = (i + K + 2 * N * N) % N;
    
82.         ex = ((i + K + 2 * N * N) / N) & 1;
    
83.         if (ex == 0) {
    
84.             Global[Z + adr] = Global[X + adr] - Y[i];
    
85.         } else {
    
86.             Global[Z + adr] = Global[X + adr] + Y[i];
    
87.         }
    
88.     }
    
89.     CountAdd +=N;
    
90. }
    
91. 
    
92. void PolyAddPermute(long N, long K, long X, NUMBERTYPE * Y, long Z) {
    
93.     long i;
    
94.     long adr;
    
95.     long ex;
    
96. 
    
97.     for (i = 0; i < N; i++) {
    
98.         adr = (i + K + 2 * N * N) % N;
    
99.         ex = ((i + K + 2 * N * N) / N) & 1;
    
100.         if (ex == 0) {
     
101.             Global[Z + adr] = Global[X + adr] + Y[i];
     
102.         } else {
     
103.             Global[Z + adr] = Global[X + adr] - Y[i];
     
104.         }
     
105.     }
     
106.     CountAdd +=N;
     
107. }
     
108. 
     
109. void PolyTransNB(long N, long L2, long K, long X) {
     
110.     long i;
     
111.     long M;
     
112.     long m1;
     
113.     long P;
     
114.     long p1;
     
115.     long Q;
     
116.     long q1;
     
117.     long adr;
     
118.     NUMBERTYPE *Temp = new NUMBERTYPE[PolyTempMax];;
     
119.     M = (long)floor(0.5 + log(N) / log(2));
     
120.     P = 1;
     
121.     Q = N / 2;
     
122.     for (m1 = 0; m1 < M; m1++) {
     
123.         for (p1 = 0; p1 < P; p1++)
     
124.             for (q1 = 0; q1 < Q; q1++) {
     
125.                 adr = q1 +  p1 * (Q<<1);
     
126.                 for (i = 0; i < L2; i++)
     
127.                     Temp[i] = Global[X + L2 * (adr +Q) + i];
     
128.                 PermutePolySub(L2, P * q1 * K, X + L2 * adr, Temp, X + L2 * (adr + Q));
     
129.                 PolyAdd(L2, X + L2 * adr, Temp, X + L2 * adr);
     
130.             }
     
131.         P <<=1;
     
132.         Q >>=1;
     
133.     }
     
134.     delete []Temp;
     
135. }
     
136. 
     
137. void PolyTransBN(long N, long L2, long K, long X) {
     
138.     long i;
     
139.     long M;
     
140.     long m1;
     
141.     long P;
     
142.     long p1;
     
143.     long Q;
     
144.     long q1;
     
145.     long adr;
     
146.     NUMBERTYPE *Temp = new NUMBERTYPE[PolyTempMax];;
     
147. 
     
148.     M = (long)floor(0.5 + log(N) /log(2));
     
149.     P = N / 2;
     
150.     Q = 1;
     
151.     for (m1 = 0; m1 < M; m1++) {
     
152.         for (p1 = 0; p1 < P; p1++)
     
153.             for (q1 = 0; q1 < Q; q1++) {
     
154.                 adr = q1 + 2 * p1 * Q;
     
155.                 for (i = 0; i < L2; i++)
     
156.                     Temp[i] = Global[X + L2 * (adr+Q)+i];
     
157.                 PolySubPermute(L2, P * q1 * K,X+L2*adr,Temp,X+L2*(adr+Q));
     
158.                 PolyAddPermute(L2, P * q1 * K,X + L2 * adr, Temp, X + L2 * adr);
     
159.             }
     
160.         P >>=1;
     
161.         Q <<=1;
     
162.     }
     
163.     delete []Temp;
     
164. }
     
165. 
     
166. void NegaConvolution(long N, long X, long H, long Y) {
     
167.     long i, k;
     
168.     long N2, N3, N4;
     
169.     long L1, L2;
     
170.     NUMBERTYPE Tp0, Tp1, Tp2;
     
171. 
     
172.     if (N==2) {
     
173.         Tp0= (Global[X] + Global[X + 1]) * Global[H];
     
174.         Tp1=(Global[H] + Global[H + 1]) * Global[X + 1];
     
175.         Tp2=(Global[H] - Global[H + 1]) * Global[X];
     
176. 
     
177.         CountMul +=3;
     
178.         Global[X] = Tp0 - Tp1;
     
179.         Global[X+1] = Tp0 - Tp2;
     
180. 
     
181.         CountAdd+=5;
     
182.     } else {
     
183.         N2=N<<1;
     
184.         N3=N2+N;
     
185.         N4=N2<<1;
     
186. 
     
187.         L2 = (long)floor(0.5 + sqrt(N));
     
188.         if (L2 * L2 != N)
     
189.             L2 = (long)floor(0.5 + sqrt(N2));
     
190.         L1 = N / L2;
     
191. 
     
192.         for (k=0;k<L1;k++)
     
193.             for (i=0;i<L2;i++) {
     
194.                 Global[Y+L2*k+i]=Global[X + L1 * i + k];
     
195.                 Global[Y+N2+L2*k+i]=Global[H + L1 * i + k];
     
196.                 Global[Y+N+L2*k+i]=0;
     
197.                 Global[Y+N3+L2*k+i]=0;
     
198.             }
     
199.         PolyTransNB(2 * L1, L2, L2 / L1, Y);
     
200.         PolyTransNB(2 * L1, L2, L2 / L1, Y + N2);
     
201.         for (k = 0; k < 2 * L1; k++)
     
202.             NegaConvolution(L2, Y + L2 * k, Y + N2 + L2 * k, Y + N4);
     
203.         PolyTransBN(2 * L1, L2, -(L2 / L1), Y);
     
204.         for (k = 0; k < L1; k++) {
     
205.             for (i = 1; i < L2; i++) {
     
206.                 Global[X + L1 * i + k] = (Global[Y + L2 * k + i] +
     
207.                                           Global[Y + N + L2 * k + (i - 1)]) / (2 * L1);
     
208.             }
     
209.             Global[X + k] = (Global[Y + L2 * k] - Global[Y + N + L2 * k + L2 - 1]) / (2 * L1);
     
210.             CountAdd+=L2+1;
     
211.         }
     
212.     }
     
213. }
     
214. 
     
215. 
     
216. void CircConvolution(long N, long X, long H) {
     
217.     long M;
     
218.     long N2;
     
219.     long X1, X2;
     
220.     NUMBERTYPE Tp0, Tp1, Tp2;
     
221. 
     
222.     N2 = 2 * N;
     
223.     M = N / 2;
     
224.     do {
     
225.         ModuloPlusMinus(M, X);
     
226.         ModuloPlusMinus(M, H);
     
227.         NegaConvolution(M, X, H, N2);
     
228.         X += M;
     
229.         H += M;
     
230.         M /= 2;
     
231.     } while (M != 1);
     
232. 
     
233.     Tp0=(Global[X]+Global[X+1])*Global[H];
     
234.     Tp2=(Global[H]-Global[H+1]);
     
235.     Tp1=Tp2*Global[X+1];
     
236.     Tp2=Tp2*Global[X];
     
237. 
     
238.     CountAdd +=5;
     
239.     CountMul +=3;
     
240.     Global[X] = Tp0 - Tp1;
     
241.     Global[X+1] = Tp0 - Tp2;
     
242.     M = 2;
     
243.     X1 = X - M;
     
244.     X2 = X;
     
245.     do {
     
246.         MulAddDiv(M, X1, X2);
     
247.         X2 -= M;
     
248.         M <<=1;
     
249.         X1 -= M;
     
250.     } while (M != N);
     
251. }
     
252. 
     
253. char* testmul(int e) {
     
254.     NUMBERTYPE* sequenceOne=new NUMBERTYPE[NMax];
     
255.     NUMBERTYPE* sequenceTwo=new NUMBERTYPE[NMax];
     
256.     NUMBERTYPE* resultFast=new NUMBERTYPE[NMax];
     
257.     NUMBERTYPE* resultSlow=new NUMBERTYPE[NMax];
     
258. 
     
259.     long N;
     
260.     long X,H;
     
261.     long i;
     
262. 
     
263.     N=2;
     
264.     N<<=e;
     
265.     for (i=0;i<N/2;i++) {
     
266.         sequenceOne[i]=9;
     
267.         sequenceTwo[i]=9;
     
268.     }
     
269. 
     
270.     for (i=N/2;i<N;i++) {
     
271.         sequenceOne[i]=0;
     
272.         sequenceTwo[i]=0;
     
273.     }
     
274. 
     
275.     cout<<N<<endl;
     
276. 
     
277. 
     
278. 
     
279. 
     
280.     X=0;
     
281.     H=N;
     
282.     Timer timer;
     
283.     timer.start();
     
284.     for (i = 0; i < N; i++) {
     
285.         Global[X + i] = sequenceOne[i];
     
286.         Global[H + i] = sequenceTwo[i];
     
287.     }
     
288.     CountAdd = 0;
     
289.     CountMul = 0;
     
290.     CircConvolution(N, X, H);
     
291.     timer.end();
     
292.     cout<<"test mul:耗时"<<timer.getSecond()<<"秒"<<endl;
     
293. //WriteResult(N, Global); // First part of Global contains the convolution
     
294.     printf("\n");
     
295. 
     
296. 
     
297. //CircConv(N, sequenceOne, sequenceTwo, resultSlow);
     
298. //WriteResult(N, resultSlow);
     
299.     cout<<"Additions:"<<CountAdd<<endl;
     
300.     cout<<"Multiplications:"<<CountMul<<endl;
     
301. 
     
302.         char* ret=(char*)malloc(sizeof(char)*(N+2));
     
303. 
     
304.         long t=0;
     
305.     for (i = 0; i < N; i++) {
     
306.         t+=Global[i];
     
307.                 ret[i]=t%10+'0';
     
308.         t/=10;
     
309.     }
     
310.         ret[i++]=0;
     
311. 
     
312.     return ret;
     
313. }
     
314. 
     
315. 
     
316. 
     
317. 
     
318. 
     
319. 
     
320. 
     
321. 
     
322. 
     
323. 
     
324. 
     
325. 
     
326. 
     
327. 
     
328. 
     
329. 
     
330. 
     
331. 
     
332. 
     
333. 
     

复制代码