(_1^%N)^n does not seem less accurate than _1^n%N _1&^31%32 _0.99518472667219682j0.098017140329560826 (_1^%32)^31 _0.99518472667219682j0.098017140329560826 _1^%32 0.99518472667219693j0.098017140329560604
--- Den fre 31/12/10 skrev Marshall Lochbaum <[email protected]>: > Fra: Marshall Lochbaum <[email protected]> > Emne: Re: [Jprogramming] Convolution using FFT > Til: "'Programming forum'" <[email protected]> > Dato: fredag 31. december 2010 04.41 > Okay, so this is not the method roots > uses. However, using this method > returns the same results: > > a=.r.1p1%32 > 5!:5<'a' > 0.99518472667219693j0.098017140329560604 > b=.a^31 > 5!:5<'b' > _0.99518472667219682j0.098017140329560826 > a-b > 1.99037j_2.22045e_16 > > Marshall > > -----Original Message----- > From: [email protected] > [mailto:[email protected]] > On Behalf Of Henry Rich > Sent: Thursday, December 30, 2010 4:24 PM > To: Programming forum > Subject: Re: [Jprogramming] Convolution using FFT > > 1p1 * 31%32 > > takes 5 more high-order bits than > > 1p1 * 1%32 > > so the absolute error in representation is that much > greater, and when you > take the sin, you should expect to lose at least that many > low-order bits. > That seems to be approximately the difference you see > here. > > IEEE floating point, with 53 bits of fraction, is at best > accurate to within > 1 part in 1e16, so we are getting almost all the accuracy > allowed by the > precision of the number system. > > Henry Rich > > > On 12/30/2010 3:35 PM, Marshall Lochbaum wrote: > > 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 > > > ---------------------------------------------------------------------- > 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
