This looks to be correct. Take, for example, a calculation of sin(pi/32) and sin(31*pi/32).
a=.1 o. 1p1*1 31%32 5!:5<'a' 0.098017140329560604 0.098017140329560826 |@-/a 2.22045e_16 Mathematica gives the exact value of the sin to be 0.098017140329560602 so the first value is quite close and the second is farther off. Can anyone post timings on polynomial multiplication using fftw? It would seem that a package like this would be the best option for implementing convolution in extended integers. Marshall -----Original Message----- From: [email protected] [mailto:[email protected]] On Behalf Of Henry Rich Sent: Thursday, December 30, 2010 1:07 PM To: Programming forum Subject: Re: [Jprogramming] Convolution using FFT It seems to me that if you compute a =. _1 ^ 2 % N, and then take that value and raise it to a high power, you are going to propagate and multiply any original roundoff error of computing a. Whereas, if you compute _1 ^ (n % N), you are going to get a more accurate result since n % N is exact (N is a power of 2). I admit that numerical analysis is not my best subject. I reasoned that by computing 1/8 or a revolution and using symmetry to produce two reflection of it, it would get a slightly more accurate result for the part it calculated, and then using symmetry would absolutely remove any bias in the coefficients. The results seem to bear this out, in that the error terms after fft/ifft are smaller. I don't think this is an error in the J implementation, just an artifact of floating-point arithmetic. Henry Rich On 12/30/2010 12:35 PM, Bo Jacoby wrote: > http://www.jsoftware.com/jwiki/Essays/FFT says that a long version of the verb roots has better accuracy for N>=8 . > > Has _1&^ for computing complex-valued fractional powers of _1 really low accuracy? > > If that is the case it should be fixed in the implementation of ^ rather than in a work-around, I think. > > Happy new year! > > > --- Den tors 30/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: torsdag 30. december 2010 18.05 I have incorporated these >> further comments into the wiki page. >> >> The J system currently does not yet use FFT to multiply big ints. It >> should, and your code and comments will be helpful in that speed-up. >> >> >> >> ----- Original Message ----- >> From: Henry Rich<[email protected]> >> Date: Thursday, December 30, 2010 7:32 >> Subject: Re: [Jprogramming] Convolution using FFT WAS: >> Proposal for special code voor dyadic f/.@:(g/) >> To: Programming forum<[email protected]> >> >>> The code I posted on the Wiki is for understanding, >> not for >>> performance >>> measurement. Several improvements are possible: >>> >>> 0. The code recomputes roots of _1 (4 times!) every >> time it >>> runs; a >>> serious implementation would reuse the roots 1. The final pass, >>> which produces an imaginary-only >> result, >>> would be >>> replaced by a step to halve the size of the final FFT, >> which >>> would save >>> 25% overall >>> 2. You would just use the fftw library anyway, which >> is >>> optimized >>> already and is a few times faster >>> >>> The important point, to me, is that this in an O(*^.) >> >>> algorithm. It can >>> feasibly multiply polynomials of a million terms, or >> numbers >>> with >>> millions of digits. >>> >>> I wanted to see whether it would work as well as the >> literature >>> suggests, and I think the answer is Yes. I also >> wondered >>> whether Roger >>> should use something like that for multiplication of >> extended >>> integers, >>> and I think the answer is Yes to that too, though for >> all I know >>> he's >>> already doing it. >>> >>> I haven't looked at your obta. >>> >>> I question your use of tm*sz as a figure of merit. >> If I am >>> willing to >>> trade space for time I am happy to take 3x the space >> in 0.5x the >>> time. >>> That is especially true here, where any operands that >> are big >>> enough to >>> cause a memory problem would take days for the >> computation. >>> >>> Henry Rich >>> >>> On 12/30/2010 1:58 AM, R.E. Boss wrote: >>>>> -----Oorspronkelijk bericht----- >>>>> Van: [email protected] >> [mailto:programming- >>>>> [email protected]] >> Namens Henry Rich >>>>> Verzonden: dinsdag 28 december 2010 18:42 >>>>> Aan: Programming forum >>>>> Onderwerp: Re: [Jprogramming] Convolution >> using FFT WAS: >>> Proposal for >>>>> special code voor dyadic f/.@:(g/) >>>>> >>>>> OK, all these suggestions are now in >>>>> >>>>> http://www.jsoftware.com/jwiki/Essays/FFT >>>>> >>>>> Henry Rich >>>>> >>>> >>>> Some observations. >>>> >>>> 'X Y'=. 8192 8193<@(?...@$)"(0) 1000 >>>> rank'iconvolve~ X';'iconvolve~ Y' >>>> +----+-----+----+----+ >>>> |rank|tm*sz|time|size| >>>> +----+-----+----+----+ >>>> | 0 |1.00 |1.00|1.00| >>>> +----+-----+----+----+ >>>> | 1 |5.22 |2.61|2.00| >>>> +----+-----+----+----+ >>>> NB. jump in performance at >> 2^n >>>> >>>> 'X Y' =. 2 400 ?...@$ 1000x >>>> datatype&> X;Y;(X +/ obta * Y);X >> >>> iconvolve Y >>>> extended >>>> extended >>>> extended >>>> integer >>>> NB. no extended output >>>> >>>> 'X Y' =. 2 8192 ?...@$ 1000 >>>> rank'X +/ obta * Y';'X iconvolve Y' >>>> +----+-----+----+----+ >>>> |rank|tm*sz|time|size| >>>> +----+-----+----+----+ >>>> | 0 |1.00 |6.49|1.00| >>>> +----+-----+----+----+ >>>> | 1 |1.36 |1.00|8.83| >>>> +----+-----+----+----+ >>>> NB. fast, not lean >>>> >>>> 'X Y' =. 2 8193 ?...@$ 1000 >>>> rank'X +/ obta * Y';'X iconvolve Y' >>>> +----+-----+----+-----+ >>>> |rank|tm*sz|time|size | >>>> +----+-----+----+-----+ >>>> | 0 |1.00 |2.69| 1.00| >>>> +----+-----+----+-----+ >>>> | 1 |6.56 |1.00|17.65| >>>> +----+-----+----+-----+ >>>> NB. overall less efficient >>>> >>>> 'X Y' =. 8000 9000<@(?...@$)"(0) 1000 >>>> rank'X +/ obta * Y';'X iconvolve Y' >>>> +----+-----+----+-----+ >>>> |rank|tm*sz|time|size | >>>> +----+-----+----+-----+ >>>> | 0 | 1.00|2.58| 1.00| >>>> +----+-----+----+-----+ >>>> | 1 |13.26|1.00|34.15| >>>> +----+-----+----+-----+ >>>> NB. especially for arbitrary >> input >>>> >>>> 'X Y' =. 2 16384 ?...@$ 1000 >>>> rank'X +/ obta * Y';'X iconvolve Y' >>>> +----+-----+----+----+ >>>> |rank|tm*sz|time|size| >>>> +----+-----+----+----+ >>>> | 1 |1.41 |9.40|1.00| >>>> +----+-----+----+----+ >>>> | 0 |1.00 |1.00|6.65| >>>> +----+-----+----+----+ >>>> NB. finally more efficient >>>> >>>> 'X Y' =. 16000 17000<@(?...@$)"(0) 1000 >>>> rank'X +/ obta * Y';'X iconvolve Y' >>>> +----+-----+----+-----+ >>>> |rank|tm*sz|time|size | >>>> +----+-----+----+-----+ >>>> | 0 |1.00 |4.13| 1.00| >>>> +----+-----+----+-----+ >>>> | 1 |5.14 |1.00|21.21| >>>> +----+-----+----+-----+ >>>> NB. However ... >>>> >>>> If f/.@:(g/) would be replaced by special code >> like obta, >>> figures would be >>>> even (slightly?) more in favor of obta. >> --------------------------------------------------------------------- >> - 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
