Re: 2.2.10 four times slower than 2.2.8

["Karim Belabas" <Karim.Belabas@math.u-bordeaux1.fr>]

* 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]