# 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 is too large
for this implementation. >>
Prime


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 is too large for
this implementation. >>
{Prime, 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 Mar 21, 2012 at 22:53
• @Mr.Wizard Thanks for your support. I really have many questions but not too much time to formulate them. Mar 21, 2012 at 22:57
• Really a nice question. I doubt any answer (and I have none) will deserve more upvotes than the question. Mar 22, 2012 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. Mar 22, 2012 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]][[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. Mar 22, 2012 at 3:19

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  Mar 22, 2012 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. Mar 22, 2012 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 ? Mar 23, 2012 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? Mar 23, 2012 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. Mar 25, 2012 at 6:54
• @Artes See edited response. It's not much, but it's all I have to offer right now. Mar 29, 2012 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. Mar 29, 2012 at 22:39
• @DanielLichtblau oh! How much I do enjoy reading your answers :) May 24, 2013 at 20:45

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

• Since when is 2^42 about 4.39805 times larger than 10^13?? May 4 at 16:51
• The results are cached (so 9.*10^-6 is wrong) and calculated upwards if you do +1, also on my windows 11 with Intel it does not fail on 10^13. In fact, it prints 323780508946331. 10^14 also works. And 10^15. And even 10^16. Did not check further. May 25 at 17:08
• My statement Prime` uses memoization is another wording for if run a second time, it is lightning-fast because is memorizes the numbers. This is for sure a property, an attribute of the local machine that hosts Mathematica. May 26 at 18:54