I including the following at least as a comment: roots1=: (_1^2%]) ^ i...@-:
As well: roots 4 |domain error: roots | roots 4 ----- Original Message ----- From: Henry Rich <[email protected]> Date: Tuesday, December 28, 2010 7:22 Subject: Re: [Jprogramming] Convolution using FFT WAS: Proposal for special code voor dyadic f/.@:(g/) To: Programming forum <[email protected]> > Actually roots as given produces some numerical > inaccuracies. A better > version is > > roots =: +...@]^:(0>[) [: (, j.) [: (, (j.~%:0.5) , |."1&.:+...@}.) [: ^ > i.@(%&8) * 0j2p1 % ] > > which calculates the exponentials over a limited range and then > uses > symmetry. This reduces the errors. > > I have put this and Roger's improved carry-propagator into the Wiki. > > Henry Rich > > On 12/28/2010 6:31 AM, Bo Jacoby wrote: > > 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
