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/