What is so special about Prime?

Question

When we try to evaluate Prime on big numbers (e.g. 10^13) we encounter the following issue:

Prime[10^13]

Prime::largp: Argument 10000000000000 in Prime[10000000000000] is too large 
for this implementation. >>
Prime[10000000000000]

Following this message, we can read in the documentation that the largest supported argument in Prime is typically about $2^{42}$. With a kind of divide and conquer approach, we can figure out that the maximal argument of Prime is:

OmegaPrime = 7783516045221;

1. What determines this number ? A hardware/software and/or conceptual/mathematical issue or maybe an arbitrary system cut-off ?

The problem seems to be a bit more obscure, since one encounters something like this:

Prime @ {# + 1, #, # + 1} & @ OmegaPrime

Prime::largp: Argument 7783516045222 in Prime[7783516045222] is too large for
this implementation. >>
{Prime[7783516045222], 249999997909357, 249999997909367}

(It takes more than two minutes to evaluate.)

An analog of OmegaPrime is OmegaPrimePi for PrimePi :

OmegaPrimePi = 25 10^13 -1;

I can find even bigger primes with Prime if I evaluate for example:

Select[Range[249999997909357, 25 10^13], PrimeQ] // Length

63142

Prime @ ( OmegaPrime + 63142 )

250000000000043

However I cannot evaluate PrimePi for numbers greater than OmegaPrimePi. It appears that Prime has a dynamically extensible domain while PrimePi does not.

2. How do I detect this property in advance from the system ?

I mean not to play around with e.g. Select[Range[a,b], PrimeQ], but for example to read it from Attributes or anything else.

Answer 1

Actually, I believe the issue reduced to that of implementing PrimePi[]. It is easy to implement Prime[] using PrimePi[] and FindRoot[] — in fact this is done on page 134 of Bressoud and Wagon, "A Course in Computational Number Theory". So all you need is to have a fast implementation of PrimePi[].

The first efficient way was found by Legendre in 1808. The modern approach of Lagarias, Miller and Odlyzko (1985) gives PrimePi[] for $10^{20}$, which is larger than the Mathematica implementation. All this is discussed in detail in the Bressound and Wagon book. Curiously they include a Mathematica package that implements the Lagarias, Miller and Odlyzko method, but it appears that (somewhat surprisingly) it has not been included in the Mathematica kernel.

---

J. C. Lagarias, V. S. Miller, and A. M. Odlyzko, Computing $\pi(x)$*: The Meissel-Lehmer method*, Math. Comp. 44 (1985), 537–560. MR 86h:11111

Answer 2

It's due to an implementation-dependent issue. We should try to improve on it. Has not been much clamor to do so, therefore it has not been a high priority.

--- edit ---

I've had a look at the code. It is quite intentional that the largest is around what you state (I see the constant being set to $7.783516108362\times 10^{12}$). It has to do with this being PrimePi[2.5*10^14] or thereabouts, and that is due to (unknown to me) implementation limitations. Basically it's a case of "Thar be dragons". Some day maybe I'll get a chance to look into this morass. (Some day I'll no doubt wind up in the La Brea Tar Pits.)

--- end edit ---