Note that ^@(0j2p1&%) is equal to _1^(2&%) so roots can be somewhat simplified.
--- Den tirs 28/12/10 skrev Roger Hui <[email protected]>: > Fra: Roger Hui <[email protected]> > Emne: Re: [Jprogramming] Convolution using FFT WAS: Proposal for special code > voor dyadic f/.@:(g/) > Til: "Programming forum" <[email protected]> > Dato: tirsdag 28. december 2010 08.24 > > 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)@|.) > > The following is somewhat more efficient: > > carry1=: }. @ (((0,~[)+0,])/) @ |: @ (0 10&#:) > carry =: carry1^: _ @ |. > > > > ----- Original Message ----- > From: Henry Rich <[email protected]> > Date: Monday, December 27, 2010 7:32 > Subject: [Jprogramming] Convolution using FFT WAS: Proposal > for special code voor dyadic f/.@:(g/) > To: Programming forum <[email protected]> > > > 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%20f/.%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 > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
