## What is so special about Prime?

58

9

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.

14I love the questions people come up with. +1 – Mr.Wizard – 2012-03-21T22:53:50.983

1@Mr.Wizard Thanks for your support. I really have many questions but not too much time to formulate them. – Artes – 2012-03-21T22:57:21.213

2Really a nice question. I doubt any answer (and I have none) will deserve more upvotes than the question. – Dr. belisarius – 2012-03-22T00:36:39.517

1@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 – 2012-03-22T00:48:03.670

One should certainly not expect Attributes to answer such a question, as that function deals with various general-type properties. My guess is that the only possibility beyond experimentation is looking at the documentation. And if that doesn't help, ask Wolfram tech support. – murray – 2012-03-22T00:55:12.163

@murray I gave Attributes here only for instance of what one could expect since I had checked this before and it hadn't enlightened me more. – Artes – 2012-03-22T01:01:38.160

7An 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. – 2012-03-22T03:19:17.850

@OleksandrR. There is an obvious mathematical relationship between these functions : they are "almost" inverse, I mean e.g. : Prime@PrimePi@Range[101, 105] returns 'Range[101, 105], i.e. identity, but notPrimePi@Prime. There are also some issues of their internal implementations which are certainly more obscure. – Artes – 2012-03-22T09:11:54.770

26

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 – 2012-03-22T10:19:22.167

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 – 2012-03-22T11:06:57.293

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 – 2012-03-23T10:15:28.663 23 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 ---

@Artes See edited response. It's not much, but it's all I have to offer right now. – Daniel Lichtblau – 2012-03-29T22:18:31.953

@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 – 2012-03-29T22:39:21.580

@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 – 2012-03-29T22:49:44.390

@DanielLichtblau I believe that there is a helpful comment by Oleksandr R. http://mathematica.stackexchange.com/questions/3327/what-is-so-special-about-prime#comment9022_3327. I suppose, that Prime works having defined in the system a kind of a sparse map of the values of PrimePi, and to evaluate Prime for an argument between any two consecutive values in that map it has to work in the same way as NextPrime. However I may be completely wrong with this.

– Artes – 2012-03-29T23:04:05.903