I finished rewriting and cleaning up my code. There is now a function mpfr_poly_mul(mpfr_poly_t res, mpfr_poly_t pol1, mpfr_poly_t pol2, ulong guard_bits).
If pol1 and pol2 are computed to sufficient precision in the first place (i.e. a bit more than the precision you want in the end) and guard_bits is set to some small quantity, usually 8 is enough even for diabolical examples, then this function handles all the "interesting" cases, including the challenges presented above. Certainly it can do much better than classical multiplication, both speed wise (it'll use the FFT if the problem gets big enough) and accuracy-wise (in a precisely defined way, which I can't be bothered going into here). The guard_bits actually does some filtering on the output of the classical and FFT multiplications which are used internally, removing spurious data below some threshold. This allows it to handle the examples rjf provided. Basically the algorithm automatically chooses a newton box and scaling for each polynomial and uses the first few tricks in Joris' paper. All the infrastructure is now there for anyone who wants to implement Joris' algorithm in full (which should now be quite straightforward), although I don't know any interesting examples where it is needed at this point. Certainly the code now does everything *I* want it to do. I've cleaned up all the test code and profiling code too, and committed to my flint2 repo. Again if anyone is interested in implementing any other interesting RR[x] algorithms, let me know. As for speeding it up further, Kronecker Segmentation, and even direct use of an integer FFT could be used to speed it up further. Also karatsuba and toom cook algorithms could be implemented. But I have no plans to work on those improvements just now. These improvements could be added quite easily to the inner routine _mpfr_poly_mul_inplace which finally does the actual multiplication once all the other jazz has been taken care of. It is defined in mpfr_poly/mul.c. Bill. On May 3, 2:43 am, Bill Hart <goodwillh...@googlemail.com> wrote: > That should say arbitrary exponents, not arbitrary precision. > > On May 3, 2:36 am, Bill Hart <goodwillh...@googlemail.com> wrote: > > > > > > > This thread is about multiple precision floating point arithmetic. > > What have machine floats got to do with it? > > > I'm using mpfr, which is what Sage uses. It has guaranteed rounding > > for *arbitrary precision* floats with essentially arbitrary precision > > exponents (there is a limit of course). > > > There's no need to even think about what happens when you switch from > > machine doubles to multiple precision, because the writers of the mpfr > > library already thought it through already. Their mpfr_t type just > > works. > > > Bill. > > > On May 3, 12:06 am, rjf <fate...@gmail.com> wrote: > > > > On May 2, 9:02 am, Bill Hart <goodwillh...@googlemail.com> wrote: > > > > > On May 2, 4:14 pm, rjf <fate...@gmail.com> wrote: > > > > > > I repeat, > > > > > > The interesting cases are obvious those which are not covered. > > > > > Sorry, I don't know what you mean. Are you saying that by definition > > > > they are interesting because they are not covered by Joris' algorithm, > > > > whatever they may be? > > > > I haven't looked at Joris' algorithm, and if all you are doing is > > > copying > > > what he has done, that might be better than making up something to do > > > something that you haven't defined. I assumed Joris defined what > > > he is doing. > > > > > > I don't know what your fix is, nor do I especially care, but I gather > > > > > that, now, at least your "stable" word is meant to indicate something > > > > > like a small bound in the maximum over all coefficients of the > > > > > difference in relative error between the true result and the computed > > > > > result. > > > > > That sounds reasonable as a definition to me. However it isn't > > > > precisely the measure Joris defines. > > > > > > I have no reason to believe this is an especially relevant measure, > > > > > since some coefficients (especially the first and the last) are > > > > > probably far more important and, incidentally, far easier to compute. > > > > > > Here are some more cases. > > > ... snip.. > > > > > Here is another > > > > > > p=1.7976931348623157E308; > > > > > q= 10*x > > > > > > What do you do when the coefficients overflow? > > > > > I actually don't understand what you mean. Why would there be an > > > > overflow? > > > > there would be an overflow if you are using machine floating point > > > numbers, > > > since p is approximately the largest double-float, and 10*p cannot be > > > represented > > > in a machine double-float. > > > > I'm missing something important here. I'm using floating > > > > > point and the exponents can be 64 bits or something like that. There > > > > should be no overflow. > > > > Really? So you are not using IEEE double-floats? > > > > What, then, do you do if the number exceeds whatever bounds you have > > > for your floats? > > > > ... snip... > > > > In fact, what does Sage do? Probably you can't say, because it > > > doesn't do the same thing > > > in various circumstances. For example, you could cause a floating > > > point overflow in the > > > midst of some computation with Maxima. > > > In that case it might depend on which Lisp that Maxima was running in. > > > Or in what operating system, or what hardware. > > > > And if the overflow occurred somewhere else, perhaps in Python or C, > > > maybe something > > > yet different. > > > > RJF > > > > -- > > > 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
