Interprocedural optimization

Question

Lets say I want to build a list of the first 1000000 primes. This is simple enough, I can just do:

AbsoluteTiming[plist = Table[Prime[i], {i, 1000000}];] and it takes 4.8 seconds.

I can also find out of the number of primes less than 15485863 (the million$^{th}$ prime) using :

AbsoluteTiming[PrimePi[15485863]] which takes negligible time.

My impression is that there is no way to get a precise PrimePi[n] result without actually finding all the primes less than n. This means that to create a list of the first 1000000 primes, it should take negligible time.

I believe though that instead of reusing the intermediaries generated from previous Prime[n] calls, it instead starts from scratch every time, and this is why it takes so long to run.

While I can write my own function to memoize this, I don't see how this could not be built in. Additionally, while in this case I can memoise the process, that is something that I cannot do in general.

I assume that in many cases involving table there would be significant savings from saving partial terms from the previous iterations, so I assume that there is a way to do this.

Is this a way to enable these optimizations, or is it just wishful thinking?

Update: Few things: - Regardless of how Prime or PrimePi is implemented, taking 4+ seconds to generate the first million primes is rather ridiculous (I have a simple little sieve in C# that can do just over 50 million primes in that time) and I think that it stands to reason that this is because it is doing the same work over and over an over again.

Perhaps this is a better example:
AbsoluteTiming[list1 = Table[Sum[i,i,10000000+q],{q,1,10000}];]
This took 95 seconds on my computer. That is more than a little excessive. The only possible way that I can think of that it took that long was that for every iteration, it optimized it (hopefully just once, but less hopefully, every time) to $1/2 (1000000 + q) (1000001 + q)$ and then at every time evaluated it.

If I instead do AbsoluteTiming[ list1 = Table[1/2 (1000000 + q) (1000001 + q), {q, 1, 10000}];] it takes just .06 seconds (much better). Is there a way to get it to always turn the first into the second? Additionally, I would tend to doubt that it realizes that the different in terms forms a linear equation and that it would be easier to do an two additions, than two divisions, three multiplications and a division.

Answer 1

The implementation of Prime and PrimePi are not as you describe.

The implementation notes say:

Prime and PrimePi use sparse caching and sieving. For large n, the Lagarias-Miller-Odlyzko algorithm for PrimePi is used, based on asymptotic estimates of the density of primes, and is inverted to give Prime.

I don't know the details of this algorithm but you can get a sense of how a prime-counting function works (Wikipedia).

Surely not all primes are calculated in order.

---

Regarding your update you seem to be asking more than one question.

Prime is not set up for rapid generation of sequential primes, but instead fairly rapid generation of arbitrarily indexed primes.

The second question seems to be about optimizing evaluation order. I don't believe there is a general solution for automatic optimization, but in this case a well placed Evaluate works, causing the Sum to first be evaluated symbolically:

AbsoluteTiming[
 list1 = Table[Evaluate[Sum[i, i, 10000000 + q]], {q, 1, 10000}];
]

{0.0160009, Null}