math - Is there a way to find the approximate value of the nth prime?

Question

Is there a function which will return the approximate value of the n th prime? I think this would be something like an approximate inverse prime counting function. For example, if I gave this function 25 it would return a number around 100, or if I gave this function 1000 it would return a number around 8000. I don't care if the number returned is prime or not, but I do want it to be fast (so no generating the first n prime numbers to return the n th.)

I would like this so that I can generate the first n prime numbers using a sieve (Eratosthenes or Atkin). Therefore, the approximation for n th the would ideally never underestimate the value of the actual n th prime.

(Update: see my answer for a good method of finding the upper bound of the n th prime number.)

Answer 1

Tighter bounds:

static const unsigned short primes_small[] = {0,2,3,5,7,11};

static unsigned long nth_prime_upper(unsigned long n) {
  double fn = (double) n;
  double flogn, flog2n, upper;
  if (n < 6)  return primes_small[n];
  flogn  = log(n);
  flog2n = log(flogn);

  if      (n >= 688383)    /* Dusart 2010 page 2 */
    upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn));
  else if (n >= 178974)    /* Dusart 2010 page 7 */
    upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.95)/flogn));
  else if (n >=  39017)    /* Dusart 1999 page 14 */
    upper = fn * (flogn + flog2n - 0.9484);
  else                    /* Modified from Robin 1983 for 6-39016 _only_ */
    upper = fn * ( flogn  +  0.6000 * flog2n );

  if (upper >= (double) ULONG_MAX) {
     /* Adjust this as needed for your type and exception method */
    if (n <= 425656284035217743UL) return 18446744073709551557UL;
    fprintf(stderr, "nth_prime_upper overflow
"; exit(-1);
  }

  return (unsigned long) ceil(upper);
}

These should not ever be less than the actual nth_prime, should work for any 64-bit input, and be an order of magnitude or more closer than the formula from Robin given earlier (or Wimblik's complicated range-limited formula). For my use I have a slightly larger small primes table so can tighten up the last estimate a bit more. Technically from the formulas we could use floor() instead of ceil() but I worry about precision.

Edit: Another option for improving this a bit is implementing good prime count bounds (e.g. Axler 2014) and doing a binary search on them. My code for this method takes ~10x longer than the above (still running in under a millisecond), but can reduce the error percentage by an order of magnitude.

If you want an estimate for the nth prime, you can do:

• Cipolla 1902 (see Dusart 1999 page 12 or this paper. I find three terms (m=2) plus a third order correction factor to be useful, but with more terms it oscillates too much. The formula shown in the Wikipedia link is this formula (with m=2). Using the two-term inverse li or inverse Riemann R below will give better results.
• Calculate the Dusart 2010 upper and lower bounds and average the results. Not too bad, though I suspect using a weighted average will work better as the bounds aren't equally tight.
• li^{-1}(n) Since li(n) is a decent approximation to the prime count, the inverse is a decent nth_prime approximation. This, and all the rest, can be done fairly quickly as a binary search on the function.
• li^{-1}(n) + li^{-1}(sqrt(n))/4 Closer, since this is getting closer to R(n)
• R^{-1} The inverse Riemann R function is the closest average approximation I know that's reasonable.

Lastly, if you have a very fast prime count method such as one of the LMO implementations (there are three open source implementations now), you can write a fast precise nth_prime method. Computing the 10^10th prime can be done in a few milliseconds, and the 10^13th in a couple seconds (on a modern fast machine). The approximations are extremely fast at all sizes and work for far larger numbers, but everyone has a different idea of what "large" means.