Karim Belabas on Sat, 27 May 2006 17:19:52 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: 2.2.10 four times slower than 2.2.8 |
* Ilya Zakharevich [2006-05-26 05:25]: > I'm testing the following code, and it is 4.5x slower with current > version than with an ancient. The difference happened between 2.2.8 > and 2.2.10. I thought that new numeric integration code was quickier > than the old one... > > \r J:\test-programs\pari\bruijn.gp > \p15 > 5 > # > bruijn_rho(5) > \q 1) The new code is indeed faster at "high" precision. Try your script at the default precision (on my machine, it becomes 5 times slower) 2) whenever intnum is used intensively using intnuminit() often helps. Here, it doesn't help much. 3) intnumromb is still around, you can use it instead of intnum A few tests at very low precision (10 to 20 digits) indicate that intnum is indeed much slower than intnumromb at these accuracies. Using intnuminit yields a 10-fold speed-up, but it's still 15 times slower to compute trivial integrals: (16:21) gp > \p realprecision = 19 significant digits (15 digits displayed) (16:21) gp > for(j=1,10^2, intnum(i=1,2,i^2)) time = 239 ms. (16:21) gp > tab=intnuminit(1,2); time = 1 ms. (16:21) gp > for(j=1,10^2, intnum(i=1,2,i^2,tab)) time = 33 ms. (16:21) gp > for(j=1,10^2, intnumromb(i=1,2,i^2)) time = 2 ms. The difference is less accute for not-so-trivial function, but still largely in favor of Romberg: (16:21) gp > for(j=1,10,intnumromb(i=2,3,zeta(i))) time = 143 ms. (16:21) gp > for(j=1,10,intnum(i=2,3,zeta(i))) time = 745 ms. On the other hand, behaviour on "non-compact" intervals looks much better now: (16:22) gp > intnum(i=2,[1],1/i^2); time = 7 ms. (16:22) gp > intnumromb(i=2,10^5,1/i^2); time = 4,038 ms. >From \p28 on (e.g at default accuracy) the new code is consistently faster. So it looks like we should use the old Romberg when realprecision is 19 or less AND we have ordinary bounds, and the new code otherwise. Cheers, K.B. P.S: The following script implements Marsaglia's method, as described in Bach and Peralta ( Math Comp. 65 (216), 1996, 1701-1715 ) and is orders of magnitude faster than direct numerical integration rho(u) = { local(n, N, C0, C); if(u<0, return(0)); if(u<=1, return(1)); if(u<=2, return(1 - log(u))); N = ceil( default(realprecision) * log(10)/log(2) ); n = ceil(u); C0 = 1 - log(2); C = vector(N, i, 1. / (i << i)); for (k = 3, n, C = vector(N, i, C0 / (i*k^i) + sum(j = 1, i-1, C[j] / (i*k^(i-j)))); C0 = sum(j = 1, N, C[j]/(j+1)) / (k-1); ); u = n-u; 0. + C0 + sum(i = 1, N, C[i] * u^i) } (16:35) gp > rho(5) time = 20 ms. %1 = 0.0003547247004560397298338945103 -- Karim Belabas Tel: (+33) (0)5 40 00 26 17 Universite Bordeaux 1 Fax: (+33) (0)5 40 00 69 50 351, cours de la Liberation http://www.math.u-bordeaux.fr/~belabas/ F-33405 Talence (France) http://pari.math.u-bordeaux.fr/ [PARI/GP]