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
