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.

UPGRADE version 12.1

Domains of Prime as well as PrimePi have been significantly extended and now

OmegaPrimePi = 2^63 - 1

 9223372036854775807

It is a much larger domain than formerly:

OmegaPrimePi/(25 10^13 - 1) // N

36893.5

Now OmegaPrime is simply the value of PrimePi on OmegaPrimePi, namely

OmegaPrime = 216289611853439384;

and the domain of Prime is extensible as before, i.e. evaluating on the fresh kernel

Prime[OmegaPrime + 1]

 Prime[216289611853439385]

Prime[OmegaPrime]

9223372036854775783

this took about $6$ minutes to evaluate on my computer. Prime[OmegaPrime] is the gratest prime number below OmegaPrimePi.

Prime[OmegaPrime + 1]

9223372036854775837

Since version 12.1 PrimePi works with several algorithms which can be chosen with Method.

Answers by Andrzej Kozłowski and Daniel Lichtblau were both helpful, although they haven't explained the origin of OmegaPrime.
Now with the last upgrade this issue becomes clear.

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 ---

Answer 3

I tried to reproduce this problem on my iMac and I did get somewhat a different result.

I entered

Prime[10^13]

and got

like you did.

I did make use of the offered assistance by Mathematica:

And got ref/message/Prime/largp for Prime, PrimPi and ZetaZero. This contained the message

The largest supported argument in Prime is typically about 2^42.

So I entered

Prime[2^42]

and got

138666449011757

But some little check shows that the given limiting number is greater than the number with the error message. The upper limit given is about 4.39805 the number with the fail message of Your question posed. The example number in the largp message is 10^20 on my iMac.

I iterated the problem and found that the problem occurs between 7 10^12 and 8 10^12. That confirms the answer of @daniel-lichtblau roughly.

My results are:

• seems that Prime uses memoization. That hints toward recursive calculations.
• seems that fresh kernel is a meaningful hint for objective judgements.
• seems that this is due to performance choices hidden before the users.
• seems to be related with timing settings, but implicit as largp shows.
• this might change is Mathematica flags the computer it is installed on as capable for parallelization, or the type of parallel computing that is implemented in Mathematica.

To my astonishment on this question I got:

AbsoluteTiming[Prime[10^12*7]]

(* {9.*10^-6, 224066643897983} *)

Oh I am really fast!!!