Karim BELABAS on Tue, 18 Apr 2000 14:02:48 +0200 (MET DST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: padic multiplication bug |
[Dan Gordon:] > I've run into an odd bug multiplying p-adic polynomials in newer > versions of GP. Let > > f = (1 + O(2^20))*y^9 + (2^14 + O(2^34))*y^7 + (2^54 + O(2^74)) > g = (1 + O(2^20))*y^12 + (2^92 + O(2^102)) > > The correct answer is: > > f*g = (1 + O(2^20))*y^21 + (2^14 + O(2^34))*y^19 + > (2^54 + O(2^74))*y^12 + > (2^92 + O(2^102))*y^9 + (2^106 + O(2^116))*y^7 + (2^146 + O(2^156)) > > In 2.0.19 (and versions at least as far back as 2.0.16), I get: > > f*g = (1 + O(2^20))*y^21 + (2^14 + O(2^34))*y^19 + > O(2^20)*y^14 + O(2^34)*y^12 + > (2^92 + O(2^102))*y^9 + (2^106 + O(2^116))*y^7 + (2^146 + O(2^156)) > > with wrong coefficients of y^14 and y^12. The coefficients are not wrong [e.g we indeed have 2^54 + O(2^74) = O(2^20)], although they are not as precise as naive multiplication could provide > Can anyone tell me what's going on? Karatsuba multiplication (kicking in when both polynomials have degree > 9): in blockwise notation, (A X + B) * (C X + D) is computed as (A*C) X^2 + ((A+B)*(C+D) - (A*C) - (B*D)) X + (B*D) Computation of the middle term involve cancellations which decrease the precision. It is possible (though a bit hackish, not documented, and in general not advised) to explicitly disable Karatsuba multiplication in library mode. You can do it in the following way under GP: install(setmulpol, vL); install(setsqpol,vL); /* use naive multiplication and squarings if degree < n */ setlim(n) = setmulpol(n); setmulpol(n) (13:16) gp > P(n) = Pol(vector(n,i,random)) (13:16) gp > a = P(1000); b = P(1000); (13:16) gp > a * b; time = 624 ms. (13:16) gp > setlim(10^6) (13:16) gp > a * b; time = 2,243 ms. the gain gets more important as the coefficients ring gets more complicated: (13:16) gp > P(n) = Pol(vector(n,i,random + O(2^100))) (13:16) gp > a = P(1000); b = P(1000); time = 3,412 ms. (13:16) gp > setlim(10^6) (13:16) gp > a * b; time = 17,605 ms. [ btw. this demonstrates that implementing padic operations using PARI padic type is a bad idea: t_PADIC is highly inefficient compared to t_INT or even t_INTMOD, which does mostly the same job (except it doesn't allow denominators). ] > With later versions, factorpadic easily handles this and other examples > we had. Is there any documentation about what improvements have > been made to factorpadic? Xavier already answered that one, but I have a side remark: the file CHANGES in the distribution toplevel (tersely) documents most changes [and so does "cvs log" for all changes made after the cvs server started (1999/09/16)] Searching that file incrementally for "factorpadic" (starting from top), the first item I find is XR 56- updated factorpadic to use new Round 4 [section: Done for version 2.0.18.beta (released 20/12/1999):] (cvs makes it easy to check the precise differences between the versions, esp. comments in the code...). The AUTHORS file provides a way of contacting the author of the changes [no author = done by the maintainer]. This is not to say that such matters shouldn't be discussed on the list, on the contrary ! Cheers, Karim. __ Karim Belabas email: Karim.Belabas@math.u-psud.fr Dep. de Mathematiques, Bat. 425 Universite Paris-Sud Tel: (00 33) 1 69 15 57 48 F-91405 Orsay (France) Fax: (00 33) 1 69 15 60 19 -- PARI/GP Home Page: http://www.parigp-home.de/