Karim BELABAS on Tue, 25 Feb 2003 19:53:06 +0100 (MET) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Fast substitutions in Pari |
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/