数学研发网设为首页收藏本站
切换到窄版

 找回密码
 欢迎注册
数学研发论坛»论坛 › 【算法研发】 › 算法交流 › 一个循环卷积的问题
1234下一页
返回列表 发新帖
楼主: 无心人

[讨论] 一个循环卷积的问题

[复制链接]
无心人
 楼主| 发表于 2008-5-21 20:53:36 | 显示全部楼层
麻烦你把那个PDF发上来
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

ssikkiss
发表于 2008-5-21 20:54:58 | 显示全部楼层
或者你在网上搜索Winograd WFTA也能找到相应的信息 在WFTA中,Winograd 为了计算FFT,把大N点FFT转化成小n的循环卷积问题, 而我们要做乘法,就是一个循环卷积问题 过去我们是用FFT来做的,也就是把循环卷积问题转化为FFT的计算 现在,如果我们对Winograd的那个"小n点快速卷及算法"弄懂了,就不必绕圈子了!
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2008-5-21 21:12:20 | 显示全部楼层
那个和FNT关系不大 那个算法是把高度复合数N点变换 分解成小n变换，高效的2-12点的变换
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

ssikkiss
发表于 2008-5-21 21:13:27 | 显示全部楼层
发了3次!发不了!太大了,9MB,
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

ssikkiss
发表于 2008-5-21 21:18:03 | 显示全部楼层
那个和FNT关系不大 那个算法是把高度复合数N点变换 分解成小n变换，高效的2-12点的变换
你说的这个是针对整个WFTA,而我说的是WFTA中Winograd提出的那个 小n点循环卷及算法 不知怎么发不上来文件,我把它用WINRAR分成了24分,还是发不上来 要么,留个email?
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2008-5-21 21:19:44 | 显示全部楼层
yaojialin@126.com
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

ssikkiss
发表于 2008-5-21 21:20:16 | 显示全部楼层
这是用Z变换来计算卷及,不是与TAOCP上的一样么?
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

ssikkiss
发表于 2008-5-21 21:30:51 | 显示全部楼层
发了,你看看收到没有
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

无心人
 楼主| 发表于 2008-5-21 21:33:47 | 显示全部楼层
谢谢
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

ssikkiss
发表于 2010-7-2 21:47:47 | 显示全部楼层
本帖最后由 ssikkiss 于 2010-7-2 22:05 编辑 此题网上以有人做出。 我抄写了它的c++代码，编译出错，稍更改了下，可以通过，不过请自己加入main函数
  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. }
  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. }
  75. void PolySubPermute(long N, long K, long X, NUMBERTYPE * Y, long Z) {
  76. long i;
  77. long adr;
  78. long ex;
  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. }
  92. void PolyAddPermute(long N, long K, long X, NUMBERTYPE * Y, long Z) {
  93. long i;
  94. long adr;
  95. long ex;
  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. }
  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. }
  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];;
  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. }
  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;
  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];
  177. CountMul +=3;
  178. Global[X] = Tp0 - Tp1;
  179. Global[X+1] = Tp0 - Tp2;
  181. CountAdd+=5;
  182. } else {
  183. N2=N<<1;
  184. N3=N2+N;
  185. N4=N2<<1;
  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;
  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. }
  216. void CircConvolution(long N, long X, long H) {
  217. long M;
  218. long N2;
  219. long X1, X2;
  220. NUMBERTYPE Tp0, Tp1, Tp2;
  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);
  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];
  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. }
  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];
  259. long N;
  260. long X,H;
  261. long i;
  263. N=2;
  264. N<<=e;
  265. for (i=0;i<N/2;i++) {
  266. sequenceOne[i]=9;
  267. sequenceTwo[i]=9;
  268. }
  270. for (i=N/2;i<N;i++) {
  271. sequenceOne[i]=0;
  272. sequenceTwo[i]=0;
  273. }
  275. cout<<N<<endl;
  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");
  297. //CircConv(N, sequenceOne, sequenceTwo, resultSlow);
  298. //WriteResult(N, resultSlow);
  299. cout<<"Additions:"<<CountAdd<<endl;
  300. cout<<"Multiplications:"<<CountMul<<endl;
  302. char* ret=(char*)malloc(sizeof(char)*(N+2));
  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;
  312. return ret;
  313. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

1234下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-11-22 00:38 , Processed in 0.030739 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表