> 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
