> One note--cube is better as the ($~ q:@#) that we came up with. This is very cute.
----- Original Message ----- From: Marshall Lochbaum <[email protected]> Date: Monday, December 27, 2010 8:00 Subject: Re: [Jprogramming] Convolution using FFT WAS: Proposal for special code voor dyadic f/.@:(g/) To: 'Programming forum' <[email protected]> > 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. > > > > ) ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
