Note that f=....@+ is its own inverse. (F^:_1 is equal to F). This makes life easier because it removes the distinction between fourier transform and inverse fourier transform.
--- Den tirs 28/12/10 skrev Henry Rich <[email protected]>: > Fra: Henry Rich <[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 18.42 > OK, all these suggestions are now in > > http://www.jsoftware.com/jwiki/Essays/FFT > > Henry Rich > > On 12/28/2010 12:07 PM, Roger Hui wrote: > > 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 > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
