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.
>
> )
>
>
>
>
>
> R.E. Boss
>
>
>
>
>
>
>
> ----------------------------------------------------------------------
> 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
