n-th prime is too slow. Consider using primecount as a backend

[galanakis]

The n-th prime function of julia is actually really slow.
For example in julia prime(10^7) runs in 14 seconds in my machine.

There is a well known C++ library with a 2-clause BSD license called primecount which can calculate
the prime(10^14) in 2 seconds (of course prime(10^7) is instantaneous.

Mathematica's Prime function is comparable to primecount (it is actually 2 times slower).

So given that, why not use primecount as a backend for the n-th prime function?

I could help with a patch if something like this sounds interesting.

[oscardssmith]

It would be much better (in my opinion) to get a fastish Julia implementation. The current implementation is a loop that calls nextprime, which is very bad. The decent simple implementation would be a segmented prime sieve that counts up as it finishes the segment. This should be relatively simple to implement and fast.

[galanakis]

Thanks for the response. I am interested in implementing the Meissel-Lehmer algorithm. It will take months of part-time work.
Is anyone else working on an implementation of prime-count?

[oscardssmith]

#102 is my 1 hour effort. It computes prime(10^7) in .14 seconds which isn't optimal, but is a lot better.

[galanakis]

Yes, it is certainly much better. However primecount is instantaneous even at 10e14, which is really impressive. I am wondering if there is interest in an implementation of the Messel - Lehmer method.

https://www.ams.org/journals/mcom/1996-65-213/S0025-5718-96-00674-6/S0025-5718-96-00674-6.pdf

Is having best of class algorithms in the scope of this project?

[oscardssmith]

I'd consider it in scope (although finding a reviewer with sufficient knowledge to approve it might be hard). Also, it seems a little premature to me given that we don't even have a segmented prime sieve yet (or efficient primality checking).

[galanakis]

It is true that the algorithm needs a segmented sieve and it is a priority. I could work on that.
Is there a specification or a potential interface for a segmented sieve? Also, it seems that the author of the following project

https://github.com/haampie/FastPrimeSieve.jl

is considering to push it to Primes.jl but is currently lacking an interface.

In general what would be the plan for a segmented prime sieve?

[oscardssmith]

@haampie would it be reasonable to make Primes.jl depend of FastPrimeSieve?

[galanakis]

If I might, could I ask again what would be the objection in implementing an API to the primecount library? It is very actively maintained (>4000 commits in 5 days) and it is faster than mathematica it self. What are the chances that a good implementation in pure julia will be within an order of magnitude from this?

[oscardssmith]

Honestly, I think it is probably too big of a dependency, but if you want to use this in Julia, https://github.com/JuliaBinaryWrappers/primecount_jll.jl already exists as a wrapper library. I don't think there is any technical reason why a Julia version would be slower, although it would require a significant amount of effort. One potentially interesting approach would be to use would be a GPU based segmented sieve (eg https://github.com/curtisseizert/CUDASieve).

[haampie]

#87 notice this PR, but I did not add many tests :( so that's probably why it was not merged

[haampie]

I'm also afraid that the code is too hard to follow, so that if there's ever an issue with it, I'm the one to fix it, but I may not have the time for it

[oscardssmith]

I think it's possible to hit a better point on the speed/complexity graph. I'll take a shot at simplifying this.

[haampie]

The segmented sieve + wheel concept is very nice though :) also the hand macro based unrolled sieving loop. I vaguely remember though that codegen is not optimal compared to primesieve's inner loop (something where LLVM can't be forced to generate a jump table). May have to revisit now that Julia uses LLVM 12

[oscardssmith]

Do you have a copy of the non macro based version? IMO, that would already be a major simplification.

[oscardssmith]

If you replace the first line with function segmented_sieve(limit::Int64, L1D_CACHE_SIZE=128*1024) it's 70x faster. Your version has L1D_CACHE_SIZE as a non-const global variable.