| 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) =
{
local(l,r,lx);
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));
l;
}
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
interpreter]
> 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].
Cheers,
Karim.
__
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/