No; the fft wasn't new to me. I will have to look at your code and see what I can get from it--my approach hasn't worked.
One note--cube is better as the ($~ q:@#) that we came up with. What I was trying to do was to use the dft verb that I made that just makes a full Fourier matrix and multiplies in something like this: fft=. [: ([: step $:"_1)`dft @.(2>#...@$) |: where step takes n rows (2 in the ideal case) of Fourier series gotten from smaller transforms and combines them into a larger series. Marshall -----Original Message----- From: [email protected] [mailto:[email protected]] On Behalf Of Henry Rich Sent: Monday, December 27, 2010 10:32 AM To: Programming forum Subject: [Jprogramming] Convolution using FFT WAS: Proposal for special code voor dyadic f/.@:(g/) This is a good excuse for doing something I have long wanted to do, exploring multiplication with FFT. (Wasn't this new to you, Marshall, when I mentioned it in passing 3 weeks ago? You must have been doing research again :) ). NB. Start with the standard J FFT library: cube =.($~ #:&.<:@#) :. , roots=.^@(0j2p1&%)@* ^ i...@-:@] floop=.4 : 'for_r. i.#$x do. (y=.{."1 y) ] x=.(+/x) ,&,:"r (-/x)*y end.' fft=:(] floop&.cube 1&ro...@#) f. :. ifft ifft=:(# %~ ] floop&.cube _1&ro...@#) f. :. fft NB. Add a verb to extend each polynomial with zeros extend =: {.~ >:@>.&.(2&^.)@# NB. Define convolution for reference: convolve =. +//.@:(*/) NB. convolution using FFT would be: iconvolve =: *&.fft&extend 1 2 3 4 convolve 2 3 4 5 2 7 16 30 34 31 20 1 2 3 4 iconvolve 2 3 4 5 2j9.32587e_15 7j6.21725e_15 16j4.44089e_15 30j_1.11022e_15 34j_8.43769e_15 31j_6.21725e_15 20j_5.32907e_15 0j1.11022e_15 The real part looks right, but there is some numerical inaccuracy. That might make this approach unsuitable in general, but if we are dealing with integers, we can take the real part and round off: roundreal =: [: <. 0.5 + 9&o. iconvolve =: roundreal@(*&.fft&extend) 1 2 3 4 iconvolve 2 3 4 5 2 7 16 30 34 31 20 0 It works! There are those extra zeros to contend with. Try it on longer operands: 'a b' =. 2 1024 ?...@$ 1000 a (iconvolve -:/@,: convolve) b 1 It matches! How's the speed? 'a b' =. 2 8192 ?...@$ 1000 ts 'a convolve b' 2.61487 5.37003e8 ts 'a iconvolve b' 0.183339 4.72448e6 Pretty good. What about for bigger operands? 'a b' =. 2 16384 ?...@$ 1000 ts 'a convolve b' |limit error: convolve | a convolve b ts 'a iconvolve b' 0.405789 9.44307e6 There is one little thing: as coded, the left and right operands must extend to the same size: 0 1 iconvolve 1 2 3 4 |length error: iconvolve | 0 1 iconvolve 1 2 3 4 But we can fix that by forcing them to the same length. While we're at it, we can use the conjugate symmetry of the FFT to take the two forward FFTs at the same time: roundimag =: [: <. 0.5 + 11&o. iconvolve =: roundimag@((- _1&|.@|.@:+)@:*:@:-:&.fft)@extend@(j./)@,: 'a b' =. 2 1e6 ?...@$ 10 ts 'a iconvolve b' 26.9542 6.71095e8 That's multiplication of two million-digit polynomials in 30 seconds. I think if you wanted to treat the polynomials as numbers you could get the digits by a carry-collection pass like polytodigits =: {:"1@((0 10 #: (+ {.))/\.&.(,&0)@|.) polytodigits 3 2 4 iconvolve 8 3 4 6 2 0 0 0 0 0 0 0 0 1 1 1 8 3 2 7 4 423*26438 11183274 Henry Rich On 12/26/2010 2:13 PM, Marshall Lochbaum wrote: > I would say this is not really worth the time to make into special > code rather than just using the obta conjunction. There just don't > seem to be uses beyond polynomial multiplication. > > Also, it seems like it would be quicker to do a convolution using > Fourier > series: pad each polynomial with zeros, convert to a Fourier series, > add, and convert back. > > Marshall > > -----Original Message----- > From: [email protected] > [mailto:[email protected]] On Behalf Of R.E. Boss > Sent: Sunday, December 26, 2010 1:09 PM > To: 'Programming forum' > Subject: [Jprogramming] Poposal for special code voor dyadic f/.@:(g/) > > In > <http://www.jsoftware.com/jwiki/RE%20Boss/J-blog/Special%20code%20for% > 20f/.% > 40:g> > http://www.jsoftware.com/jwiki/RE%20Boss/J-blog/Special%20code%20for%2 > 0f/.%4 0:g I define the conjunction 'obta' - oblique at table. > > > > f obta g is equivalent to f/.@:(g/) but is much leaner and a bit faster. > > See the wiki-page for the detailed figures. > > > > > > obta=: 2 : 0 > > assert. 'dyad only' > > : > > assert.>{.(x *.&(1...@$) y) ; 'vectors only' > > if. x -: y do. (u@:(v|.)\y) , }.u@:(v|.)\.y > > else. > > 's t'=: x ,&# y > > z=. $0 > > if. s=t do. y=.|.y > > if. x-:y do. z=.(u@:(v@:(,:|.))\y) , }.u@:(v@:(,:|.))\.y > > else. NB. y=.|.y > > z=. i.0 1 > > for_j. i.&.<: s do. z=.z, ((j{.x) u@:v (-j){.y) , ((-j){.x) u@:v > j{.y end. > > end. > > elseif. s<t do. y=.|.y > > for_j. i.&.<:s do. z=.z, (j{.x) u@:v (-j){.y end. > > z=.z, |.s x&(u@:v)\y > > for_j. |.i.&.<: s do. z=.z, ((-j){.x) u@:v j{.y end. > > elseif. s>t do. y=.|.y > > for_j. i.&.<:t do. z=.z, (j{.x) u@:v (-j){.y end. > > z=.z, t (u@:v)&y\x > > for_j. |.i.&.<: t do. z=.z, ((-j){.x) u@:v j{.y end. > > end. > > end. > > ) > > > > > > R.E. Boss > > > > > > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
