I meant
extend =: >.&.(2&^.)@<:@+&# {."_1 ,:
iconvolve =: roundimag@((- _1&|.@|.@:+)@:*:@:-:&.fft)@(j./)@extend
Henry Rich
On 12/28/2010 11:48 AM, Henry Rich wrote:
> 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
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
