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<rhui...@shaw.ca>:
>
>> Fra: Roger Hui<rhui...@shaw.ca>
>> Emne: Re: [Jprogramming] Convolution using FFT WAS: Proposal for special
>> code voor dyadic f/.@:(g/)
>> Til: "Programming forum"<programming@jsoftware.com>
>> 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<henryhr...@nc.rr.com>
>> Date: Monday, December 27, 2010 7:32
>> Subject: [Jprogramming] Convolution using FFT WAS: Proposal
>> for special code voor dyadic f/.@:(g/)
>> To: Programming forum<programming@jsoftware.com>
>>
>>> 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: programming-boun...@jsoftware.com
>>>> [mailto:programming-boun...@jsoftware.com]
>> 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
