Re: Fast substitutions in Pari

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

On Mon, 24 Feb 2003, Philippe Elbaz-Vincent wrote:
> On Mon, 24 Feb 2003, Franck MICHEL wrote:
>> \ps 16;
>> u=sum(i=1,15,(-1)^i*i^2*x^i);
>> a = 3741/2996+34323/2630*I;
>> b =
>> 2771/3873+1653/2366*I;
>> c = 1264/1653+1883/2922*I;
>> d =
>> -2159/2174+91/776*I;
>> phi=a*t*(1-t*b)*I/(1-t)/(1-c*t)/(1-d*t);
>> subst(u,x,tayl
>> or(phi,t));
>> taylor(subst(u,x,phi),t);
>>
>> On my PC, subst(u,x,taylor(phi,t)) is done in 1502ms, and
>> taylor(subst(u,x,phi),t) is done in 552ms.
>
> But what version of PARI do you use ? (2.1.x I guess, could you confirm)
>
> with 2.2.5 (development CHANGES-1.646)
> I got (on a P4 3Ghz and no load!)
>
> taylor(subst(u,x,phi),t);
> time = 1,370 ms.
> subst(u,x,taylor(phi,t));
> time = 560 ms.
>
> nevertheless there seems to be something "weird" here.
> With 2.1.3 on a P2/450Mhz
>
> I got
> (18:49) gp > taylor(subst(u,x,phi),t);
> time = 880 ms.
> (18:49) gp > subst(u,x,taylor(phi,t));
> time = 2,140 ms.
>
> Any explanations ?

Yes, this is an unforeseen side effect of

2.2.3 A-2: gcd for Gaussian integers

While operating on rational functions, we compute a huge number of polynomial
contents for simplifications [ PARI is not a computer algebra system:
everything is always completely simplified¹, which is very inefficient in the
generic cases, where few simplifications occur ]. Contents are computed using
generic gcd() operations, which are now done over Q[i] in the example above,
not over Q !

I would do something like this:

  /* h^degpol(P) P(x / h) */
  rescale(P, h) =
  { local(G,H, d);
    d = poldegree(P); H = 1;
    Pol( vector(d + 1, i, G = H; H = G * h; polcoeff(P, d-(i-1)) * G) );
  }
  h = taylor(denominator(phi), t);
  z = subst(rescale(u, h), x, taylor(numerator(phi), t));
  z / h^poldegree(u);   \\ not needed in most applications !!!

On my slow PIII (1GHz), I get:

(19:16) gp > taylor(subst(u,x,phi),t); \\ regression as above !!!
time = 2,400 ms.
(19:16) gp > subst(u,x,taylor(phi,t)); \\ slow
time = 1,010 ms.

(19:16) gp > h = taylor(denominator(phi), t);
time = 0 ms.
(19:16) gp > z = subst(rescale(u, h), x, taylor(numerator(phi), t));
time = 10 ms.
(19:16) gp > z / h^poldegree(u);  \\ fast
time = 20 ms.

If this operation is really crucial, a better approach could be to multiply
by Mod(1, p) for some primes = 1 mod 4, not dividing any denominator in
sight, and reconstruct using Chinese remaindering. But you'll need a priori
bounds on the heights of the final coefficients.

Hope this helps,

    Karim.

P.S: it is not trivial to build this in, since if we substitute a genuine
rational function [ not a power series as above ], I now get

(19:19) gp > h = denominator(phi);
time = 0 ms.
(19:22) gp > z = subst(rescale(u, h), x, numerator(phi));
time = 10 ms.
(19:22) gp > z / h^poldegree(u);
  ***   user interrupt after 10mn, 10,390 ms.

Non-trivial multivariate arithmetic is inefficient in PARI :-(. Due to
constant oversimplification as explained above, and (slooow) generic
subresultant GCDs where everything should be modular [ possibly here in the
guise of evaluation / interpolation. ]. But modular algorithms are not
triggered for complex coefficients...


¹: this is not quite true. Polynomials of degree 0 are not promoted to
scalars.
-- 
Karim Belabas                    Tel: (+33) (0)1 69 15 57 48
Dép. de Mathématiques, Bât. 425  Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud             Email: Karim.Belabas@math.u-psud.fr
F-91405 Orsay (France)           http://www.math.u-psud.fr/~belabas/
--
PARI/GP Home Page: http://www.parigp-home.de/