Re: padic multiplication bug

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

[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/