# What is so special about Prime?

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.

-
I love the questions people come up with. +1 – Mr.Wizard Mar 21 '12 at 22:53
@Mr.Wizard Thanks for your support. I really have many questions but not too much time to formulate them. – Artes Mar 21 '12 at 22:57
Really a nice question. I doubt any answer (and I have none) will deserve more upvotes than the question. – Dr. belisarius Mar 22 '12 at 0:36
@belisarius Thank You ! I believe there is a need for Primes tag since M contains quite a good functionality in this field and there could appear many interesting and related questions. – Artes Mar 22 '12 at 0:48
An interesting observation is that Prime calls PrimePi many (namely, 1,013,381) times when given an argument of your OmegaPrime: nums = Reap[InternalInheritedBlock[{PrimePi}, Unprotect[PrimePi]; pp:PrimePi[n_] /; (Sow[n]; True) := pp; Protect[PrimePi]; Prime[7783516045221]]][[2, 1]]; ListLogLogPlot[nums, MaxPlotPoints -> 1000, Joined -> True] gives nearly a straight line. What this means, if anything, I have no idea, but it shows at least some concrete relationship between the two functions. – Oleksandr R. Mar 22 '12 at 3:19

## 2 Answers

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

-
However when I evaluate with a fresh kernel PrimePi@249999999909367;// AbsoluteTiming//First I get slightly bigger number than for Prime@7783516045221;// AbsoluteTiming//First  – Artes Mar 22 '12 at 10:19
In my case Prime took slightly longer: 153.04 vs 151.9 for PrimePi. Prime certainly calls on PrimePi so all this means is that AbsoluteTiming is unreliable when the difference is so small. – Andrzej Kozlowski Mar 22 '12 at 11:06
Thank you for these interesting references. Apparently Prime and PrimePi are internally related. However it is not clear how Prime works when its arguments are greater than OmegaPrime. I suspect there is communication between NextPrime and/or PrimeQ, since these functions have domains upper-bounded only by $MaxNumber. Have you got any idea why OmegaPrime is just 7783516045221 ? – Artes Mar 23 '12 at 10:15 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 ---

-
I would love to hear from you more about those internal implementation issues. Have you got any idea why the value of OmegaPrime is 7783516045221 ? There seems to be no strict relation between implementation of Prime and PrimePi. Does Prime use NextPrime (or do they both use primality tests as in PrimeQ) when arguments exceed OmegaPrime? – Artes Mar 23 '12 at 10:35
I have run into this too and wondered why Mathematica had this limit. A way to circumvent it would be some sort of mathlink/library link and C++ mixture. – nixeagle Mar 25 '12 at 6:54
@Artes See edited response. It's not much, but it's all I have to offer right now. – Daniel Lichtblau Mar 29 '12 at 22:18
@DanielLichtblau Thank You for the edit. The constant you've mentioned OmegaPrime1 is such that Prime[OmegaPrime1]= 249999999999977 is the largest prime smaller than OmegaPrimePi, so it is also such that PrimePi[2.5*10^14]=OmegaPrime1, however I have OmegaPrime = 7783516045221 in Mathematica ver 7 as well as in ver. 8 (under Windows) so it is still a bit more puzzling why it is different. – Artes Mar 29 '12 at 22:39
@DanielLichtblau I have to emphasize that the domain of Prime is extensible during a session (I mentioned that in the post). I mean that Prime[7783516108362] returns the same with the fresh kernel. However after some playing around with e.g. NextPrime` I can evaluate it getting what I wrote before in my comment. – Artes Mar 29 '12 at 22:49