The power series method would apply to any polynomial whose constant term is a square (or rather, whose lowest degree monomial has the form (square coefficient)*(even power of x). Viewed as a power series such a thing is always a square, as a simple induction constructs a square root. We already have this:
sage: R.<x>=PowerSeriesRing(QQ) sage: f=R(1+2*x+3*x^2) sage: f.sqrt() 1 + x + x^2 - x^3 + 1/2*x^4 + 1/2*x^5 - 3/2*x^6 + 3/2*x^7 + 3/8*x^8 - 29/8*x^9 + 43/8*x^10 - 11/8*x^11 - 155/16*x^12 + 325/16*x^13 - 215/16*x^14 - 393/16*x^15 + 9891/128*x^16 - 10269/128*x^17 - 5901/128*x^18 + 36813/128*x^19 + O(x^20) sage: g=f^2 sage: g.sqrt() 1 + 2*x + 3*x^2 sage: PolynomialRing(QQ,'X')(g.sqrt().list()) 3*X^2 + 2*X + 1 John 2008/8/18 Bill Hart <[EMAIL PROTECTED]>: > > > > On 18 Aug, 17:39, "Joel B. Mohler" <[EMAIL PROTECTED]> wrote: >> On Monday 18 August 2008 11:44:11 am Bill Hart wrote: >> >> > Assuming FLINT is in fact being used for the GCD in Z[x], the >> > implementation is only "fast" up to degree about 250. It depends on >> > your definition of fast though. :-) In a later release of FLINT, GCD >> > will be asymptotically faster, and certainly by degree 1000 will be >> > many times faster than it currently is. >> >> By "fast", I meant faster than my naive implementation of the sqrt. The ZZ >> gcd certainly seems fast enough for any problem I'm likely to have in the >> next 6 months :). The sage I'm using (3.0.6) is using the flint >> fmpz_poly_gcd. > > OK. That is currently quadratic, but with a low constant. :-) > >> >> > If I'm not mistaken, the linear algebra method is nearly cubic >> > complexity. The GCD method should be subquadratic (though only just). >> > The method using power series has the same complexity as Karatsuba >> > multiplication of polynomials, in terms of ring operations (see papers >> > of Paul Zimmermann), and the constant can be made very low by use of >> > the recursive middle product algorithm. >> >> My algorithm is most definitely quadratic (as I had already thought from the >> simple nested for-loops). > > Yes I've had a look now, and I think you are right. > >> Doubling the degree produces a nearly perfect >> quadrupling of run-time. Note that it isn't a linear system of the >> coefficients and the solution is so obvious that it's easy to find one >> unknown at a time. I'm not sure what you are calling the "power series >> method" -- I wonder if it may be what I'm doing. > > No it definitely can't be, since the latter has much better asymptotic > complexity. I'm simply referring to thinking of the polynomial as a > power series, then essentially using Newton iteration to compute the > square root, though there are many variants of this basic idea which > can cut the implied constant by a very considerable factor. > > Another method is to take the logarithm of the power series, multiply > by 0.5 and then exponentiate. > >> >> > Note that if one requires the square root of an exact square of >> > polynomials, if one uses the power series method, one only needs to >> > work to a precision about half the length of the original polynomial >> > since that is how large the square root will be. I think so anyway, >> > unless I'm missing something silly here (altogether likely). >> >> That makes sense. If you know that you have a perfect square you only need >> to >> think about half the coefficients of the square input. Of course, if you >> want to raise an error when something is not a perfect square, you'll need to >> test all the coefficients by comparing with the square of your tentative >> square root. >> > > It depends on how you implement it. As soon as you know one of the > coefficients is wrong for it to be a square, you can bail out and say > it isn't. It may not be necessary to jack the precision right up to n > for that. E.g. suppose n is just below a power of 2 and you use Newton > iteration to double the number of coefficients each time, you may > already know by the time you have just passed half way that you don't > have a square. For that matter, you may know at any time. > > It should also be possible to use a multimodular method to check it is > a square root in time something like n log (nm) where n is the degree > and m is the largest coefficient. This has a significant chance of > bailing out early if it isn't a square root and would certainly be > faster than multiplying out in full. I suspect on average O(1) primes > are needed on average to check if it isn't a square, and so this gives > an asymptotic improvement over multimodular multiplication, which > takes log m primes. > > Obviously over a general ring, you have to multiply using classical or > Karatsuba or multiplication, not all of which needs to be done to get > some of the coefficients. Of course that is a pointless observation in > special rings like Z, Z/pZ and Q where you don't have to use these > algorithms. > > And also note you only need a truncated product, since you know > already the top half of the terms are correct. Asymptotically that > doesn't save anything, but FLINT saves a small factor doing truncated > products anyway, since there is less real work to do, even when it > uses an FFT. You can combine this suggestion with the multimodular > suggestion. > > Another throw away observation is that squaring can be done faster > than a full product. FLINT does this automatically for sufficiently > large polynomials, i.e. length > 20 (FLINT 2.0 will do it for all > polynomials and I've already written some of the code for that). Again > you can use this with the multimodular and truncated product > suggestions. > > Anyhow, I've just outlined why I won't be doing this in FLINT > soon. :-) There's actually quite a lot of tricks to go through, and > there may even be an algorithm for detecting perfect squares which I > know nothing about. > > Bill. > > > --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to [email protected] To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
