I was thinking ahead to the final solution, which solves the padding
problem by laminating the two operands and then applying extend just
once, to the combined complex vector. Then extend is used only monadically.
And anyway, you have to make the operand a power-of-2 length for FFT to
work.
Henry Rich
On 12/28/2010 11:21 AM, Marshall Lochbaum wrote:
> Looking back at extend, I think it is better to make it a dyad that pads
> both numbers:
>
> extend2=: +&# {."_1 ,:
>
> iconvolve =: *&.fft/ @:extend2
>
> This removes the need for inexact padding and gives you exactly the right
> length result.
>
> Marshall
>
>
> -----Original Message-----
> From: [email protected]
> [mailto:[email protected]] On Behalf Of Henry Rich
> Sent: Tuesday, December 28, 2010 10:22 AM
> To: Programming forum
> Subject: Re: [Jprogramming] Convolution using FFT WAS: Proposal for special
> code voor dyadic f/.@:(g/)
>
> 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
>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
