Karim BELABAS on Thu, 25 May 2000 19:57:23 +0200 (MET DST)

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: This *can't* be missing, can it?

[Tom Womack:]
> What are the pari names for pi(x) [the inverse function for prime(x)] and
> pimax [the largest value for which prime(x) is a precomputed prime]?

None of these two exist... It's an easy task to program a simple version in
library mode, but quite an annoying one under GP. I'll add it to the TODO
file [more generally, coding the algorithms of Deléglise and Rivat should be
feasible and could be quite useful].

> It's trivial to write a little binary-search routine
> {
> pith(x) =
>  local(t,l,r);
>  t=1;
>  while(prime(t)<=x,t=t*2);
>  l=t/2;r=t;
>  while (r-l>1,m=(l+r)/2;if (prime(m)<=x,l=m,r=m));
>  l;
> }
> but that works badly if X is near enough to primelimit for the first while
> loop to overshoot, and feels as if it ought to be a primitive function in
> any case.

I couldn't resist hacking a wee little bit on your function. The following
one works unless x is _really_ close to primelimit [and exhibits some
deplorable hacks permitted by recent versions of the interpreter]

  myprime(m) = prime(m);

  pith(x) =

    if (x <= 2, return (x==2));
    lx = log(x);
    l = floor(x/lx);
    r = floor(l * (1 + 3/(2*lx)));
    while (r-l>1,
      m = (l+r)>>1;
      if (myprime(m)<=x, l=m, r=m));

Of course prime(n) should return a conventional answer (e.g 0) instead of
stopping the script with an error. Since it doesn't...

  /* return something > x if prime(m) overflows */
  myprime(m) = trap(,1e50, prime(m))

Now pith() works up to primelimit and possibly a little beyond.

[I'll let somebody else write a GP script dumping the relevant C code to a
file, compiling it as a shared module, and providing the function to the

> I'm trying to implement the Buhler-Gross clever algorithm for computing
> L^(r)(1) (which requires only sqrt(n) storage to do a sum up to a_n, and so
> should converge reasonably for conductors up to O(available storage ^ 4) ),
> with the eventual aim of adding an ellanalyticrank(e) function to Pari
> [would this be useful functionality?]

I definitely think it would [it might also be very useful in teaching
elllseries how to use a decent amount of memory].


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/