You can access the code here: http://selmer.warwick.ac.uk/gitweb/flint2.git?a=shortlog;h=refs/heads/FHT
You need the latest MPIR release candidate, the latest MPFR and either an x86 or x86_64 machine to run it. set the liib and include paths in the top level makefile, set the LD_LIBRARY_PATH's in flint_env, then do: source flint_env make library make check MOD=mpfr_poly Most of the relevant code is in the mpfr_poly directory. If it doesn't work for you, please understand this is development code only, therefore probably broken in many ways. Information on the Fast Hartley Transform can be found in Joerg Arndt's fxtbook: http://www.jjj.de/fxt/fxtbook.pdf Bill. On Apr 28, 2:45 am, Bill Hart <goodwillh...@googlemail.com> wrote: > Well I got the polynomial convolution working with the Fast Hartley > Transform. It seems to pass my primitive test code. > > As an example of a first timing, 1000 iterations of multiplying > polynomials of length 512 with 100 floating point bits precision per > coefficient takes 17s. > > I think I can make it take about half that time (at least), but that > is a job for tomorrow. > > I haven't looked yet to see how much precision is lost heuristically. > > Bill. > > On Apr 27, 11:49 pm, Bill Hart <goodwillh...@googlemail.com> wrote: > > > > > > > Well, I coded up a mini mpfr_poly module and Fast Hartley Transform in > > flint2 using mpfr's as coefficients. > > > It's hard to know if it is doing the right thing, there's no test code > > yet. But it compiles and doesn't segfault. :-) > > > A little bit more work required to turn it into a convolution, but as > > far as I know, no inverse transform required, as it is its own > > inverse. Maybe by tomorrow I'll have a fast convolution and we'll be > > able to answer some of these questions heuristically. > > > It won't be super efficient, no truncation, no radix 4, no special > > case for double precision coefficients, no cache friendly matrix > > fourier stuff, etc, but I did reuse the twiddle factors as much as > > possible, I think. > > > Bill. > > > On Apr 27, 5:18 pm, Bill Hart <goodwillh...@googlemail.com> wrote: > > > > Also see pages 252 and following of Polynomial and Matrix > > > Computations, Volume 1 by Dario Bini and Victor Pan, which seems to > > > answer your question in detail. > > > > Bill. > > > > On Apr 27, 4:57 pm, Bill Hart <goodwillh...@googlemail.com> wrote: > > > > > Numerical stability is not something I have any experience with. I am > > > > not sure if is is equivalent in some sense to the loss of precision > > > > which occurs when doing FFT arithmetic using a floating point FFT. > > > > > The issue seems to be accumulation of "numerical noise". There are > > > > proven bounds on how many bits are lost to this as the computation > > > > proceeds, but the general consensus has been that for "random > > > > coefficients" the numerical noise is much, much less than the proven > > > > bounds would indicate. > > > > > Zimmermann, Gaudry, Kruppa suggest that the following paper contains > > > > information about proven bounds for the floating point FFT: > > > > > Percival, C. Rapid multiplication modulo the sum > > > > and difference of highly composite numbers. Math. > > > > Comp. 72, 241 (2003), 387–395 > > > > > Paul Zimmermann would almost certainly know what the current State of > > > > the Art is. > > > > > Similar stability issues arise in matrix arithmetic with floating > > > > point coefficient. It seems that one can always prove something when > > > > you know something about the successive minima. But I am not too sure > > > > what the equivalent for an FFT would be. You might be into the area of > > > > research rather than something that is known, i.e. it might just be > > > > conventional wisdom that the FFT behaves well. > > > > > The people to ask about this are definitely Joris van der Hoeven and > > > > Andreas Enge. If you find a very nice algorithm, let me know and we'll > > > > implement it in flint2, which should have a type for polynomials with > > > > floating point coefficients (using MPFR). > > > > > Bill. > > > > > On Apr 27, 2:45 pm, Jason Grout <jason-s...@creativetrax.com> wrote: > > > > > > On 04/26/2010 10:54 PM, Robert Bradshaw wrote: > > > > > > > I should comment on this, as I wrote the code and comments in > > > > > > question. > > > > > > There actually is a fair amount of research out there on stable > > > > > > multiplication of polynomials over the real numbers, but (if I > > > > > > remember > > > > > > right, it was a while ago) there were some results to the effect > > > > > > that no > > > > > > asymptotically fast algorithm has good stability properties over all > > > > > > possible coefficients. > > > > > > If you have a moment and if you remember, could you point out a > > > > > reference or two? I searched for several hours before posting, and > > > > > couldn't find anything that seemed to address the question of > > > > > stability > > > > > satisfactorily. > > > > > > > Of course, in practice one often runs into "well > > > > > > behavied" polynomials in R[x]. (In particular, most of the interest > > > > > > is, > > > > > > understandably, about truncated power series, though there are other > > > > > > uses.) The obvious thing to do is a (non-discrete) FFT. What I was > > > > > > thinking about implementing (though I never got back to it) is > > > > > > something > > > > > > that works like this: > > > > > > > 1) Choose some a and rescale f(x) -> f(alpha*x) so all coefficients > > > > > > have > > > > > > roughly the same size. > > > > > > 2) Multiply the resulting polynomial over Z using the insanely fast > > > > > > code > > > > > > in FLINT. The resulting answer will be exact product of the > > > > > > truncated > > > > > > input, and the precision loss (conversely, required precision for no > > > > > > loss) was fairly easy to work out if I remember right . > > > > > > 3) Cast the result back into R[x]. > > > > > > Thanks! > > > > > > Jason > > > > > > -- > > > > > 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 > > > 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 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
