In fact, if you do exp(2*x) times exp(-x) you'll see the classical algorithm also sucks.
Bill. On May 1, 4:12 pm, Bill Hart <goodwillh...@googlemail.com> wrote: > Actually, in the situation where the magnitude of the terms varies > exponentially like this in a regular way, you can just evaluate your > power series at c*x for some constant c (e.g. 2 in the example you > give), then do the multiplication asymptotically fast, then evaluate > the result at x/c. > > If the magnitude doesn't vary regularly, you are screwed no matter > what algorithm you use. > > Bill. > > On May 1, 4:53 am, Robert Bradshaw <rober...@math.washington.edu> > wrote: > > > > > > > On Apr 30, 2010, at 8:00 PM, Bill Hart wrote: > > > > Finally, the figures for karatsuba with signed coefficients that vary > > > wildly. Hardly distinguishable from the figures for classical > > > multiplication. With that I finish my little study. And in the > > > inimitable words of Adam and Jamie (and with about as much "science" > > > to back it up): MYTH BUSTED!! > > > > So what have we learned: > > > > * The Sage timings for multiplication in RR[x] are not great, (please > > > note however that my original timings of the classical algorithm in > > > flint2 were incorrect due to a bug in my code; flint is actually only > > > 10-20% faster for this - the times for the Hartley Transform were also > > > misreported and should have been 9ms and 125ms for degree 1000 and > > > 10000 multiplications, respectively, still considerably faster than > > > Sage). > > > > * Karatsuba does not display noticeably worse behaviour than classical > > > multiplication. Toom-Cook is also acceptable. > > > > * There is no such thing as Fateman multiplication. > > > > * The "Fateman multiplication" code is broken. > > > > * The Hartley transform is much better than Rader-Brennen. > > > > * For fixed point arithmetic the FFT techniques are fine, but for > > > floating point they will suck when the "magnitude" of the coefficients > > > varies greatly. > > > > * To answer the question asked by the OP, yes, there is a better > > > algorithm; see the algorithm of Joris vdH > > > > * NO algorithm that we have looked at will deal with the kind of > > > cancellation talked about in the docstring which started this whole > > > thread off. > > > May I add another datapoint: rather than looking at random > > polynomials, consider power series expansions. Here I compute the > > expansion for e^x to 100 terms via classical, karatsuba, and fft over > > RR (53-bits) vs. exactly over QQ: > > > def rel_prec(approx, actual): > > return [oo if a==b else -log(abs(a-b)/abs(b), 2) for (a,b) in > > zip(approx, actual)] > > > from numpy.fft import rfft, irfft > > def fft_mul(f, g): > > n = f.degree() + g.degree() + 1 > > fft_f = rfft(list(f) + [0] * (n-f.degree())) > > fft_g = rfft(list(g) + [0] * (n-g.degree())) > > return f.parent()(irfft( fft_f * fft_g )) > > > sage: R.<x> = QQ[] > > sage: f = exp(x+O(x^100)).polynomial() > > sage: ff = f.change_ring(RR) > > > sage: sorted(rel_prec(ff*ff, f*f)) > > [52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, 52, > > 52, +Infinity, +Infinity, +Infinity, +Infinity, +Infinity, +Infinity, > > +Infinity, +Infinity, +Infinity, +Infinity, +Infinity, +Infinity, > > +Infinity, +Infinity, +Infinity, +Infinity, +Infinity, +Infinity, > > +Infinity, +Infinity, +Infinity] > > > sage: sorted(rel_prec(ff._mul_karatsuba(ff), f*f)) > > [7, 9, 10, 11, 12, 12, 12, 14, 15, 16, 20, 21, 21, 25, 25, 26, 26, 28, > > 30, 30, 33, 35, 39, 40, 41, 41, 43, 44, 45, 45, 46, 50, 50, 50, > > +Infinity, +Infinity, +Infinity, +Infinity, +Infinity] > > > sage: sorted(rel_prec(fft_mul(ff, ff), f*f)) > > [-61, -54, -46, -44, -42, -35, -32, -24, -24, -18, -13, -12, -8, -3, > > -1, 0, 4, 7, 10, 13, 17, 24, 27, 29, 31, 33, 35, 37, 40, 43, 45, 48, > > 48, 50, 51, 52, 52, 52, +Infinity] > > > Observation: > > > - Classical -- Not bad, 1 bit loss max. > > - Karatsuba -- Things get pretty bad (especially towards the middle of > > the product). > > - FFT -- Yes, that's a negative amount of precision, i.e. off by > > orders of magnitude, for half the coefficients. It just can't handle a > > range of coefficient magnitudes greater than the working precision. > > > I don't have time to do extensive experiments or draw a conclusion, > > but the code above lends itself easily to more investigations in this > > direction. > > > - Robert > > > -- > > To post to this group, send an email to sage-devel@googlegroups.com > > To unsubscribe from this group, send an email to > > sage-devel+unsubscr...@googlegroups.com > > For more options, visit this group > > athttp://groups.google.com/group/sage-devel > > URL:http://www.sagemath.org > > -- > To post to this group, send an email to sage-devel@googlegroups.com > To unsubscribe from this group, send an email to > sage-devel+unsubscr...@googlegroups.com > For more options, visit this group athttp://groups.google.com/group/sage-devel > URL:http://www.sagemath.org -- To post to this group, send an email to sage-devel@googlegroups.com To unsubscribe from this group, send an email to sage-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URL: http://www.sagemath.org
