Oh, I see your point, if you are convolving length 4 and length 11, you
only need to go up to length 16, not 32. Yeah, that's right. So that
would be
extend3 =: >:@>.&.(2&^.)@+&# {."_1 ,:
Henry Rich
On 12/28/2010 11:37 AM, Henry Rich wrote:
> 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
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
