Here is an updated primes package: https://github.com/jlapeyre/PrimeSieve.jl
I didn't realize that the author of libprimesieve had recently released
a library libprimecount, parts of which supersede the older library. I
compared a bit with Mathematica 8, Maxima, and Numbers.jl. (don't have
access to Maple, or other systems)
The prime number functions in Mathematica are not competitive with
PrimeSieve.jl. Since Wolfram has an infinite supply of money, I guess
this is because market research shows there is no demand for
improvement, or maybe they have backward compatibility issues.
Here are some trials: (I guess this markdown will not be rendered, but
its already there!)
### Prime counting function
```
In[3]:= Timing[PrimePi[10^14]]
Out[3]= {44.4268, 3204941750802}
```
```julia
julia> @time primepi(10^14) # not fair, because 10^14 is in a table
elapsed time: 8.695e-6 seconds (120 bytes allocated)
3204941750802
julia> @time primepi(10^14; alg = :dr) # force no tables: Deleglise-Rivat
elapsed time: 0.45482745 seconds (176 bytes allocated) # bytes not
including library routines
3204941750802
```
PrimePi only accepts arguments < 2.5 * 10^14.
primepi accepts an argument up to 10^27. Deleglise first reported
primepi(10^19) in 1996 using 40 hours of processor time.
We (libprimecount) do it in 5 minutes.
There are two algorithms implemented. My logic for choosing the faster
is far from perfect, but could be made solid with effort or cleverness.
Also a third algorithm of complementary, but restricted, applicability
is easy, but i did not implement it.
### Find nth prime
```
In[1]:= Timing[Prime[10^12]]
Out[1]= {19.1772, 29996224275833}
In[2]:= Timing[Prime[10^13];]
Prime::largp: Argument 10000000000000 in Prime[10000000000000]
is too large for this implementation.
```
```julia
julia> @time nthprime(10^12)
elapsed time: 0.26973466 seconds (112 bytes allocated)
29996224275833
julia> @time nthprime(10^15)
elapsed time: 17.464598193 seconds (112 bytes allocated)
37124508045065437
```
### Next prime
I translated and modified the [Maxima code to
Julia](https://github.com/jlapeyre/PrimeSieve.jl/blob/master/src/nextprime.jl)
I compared in various loops to [Hans Borchers'
Numbers.jl](https://github.com/hwborchers/Numbers.jl/blob/master/src/primes.jl#L70-L98)
The former is faster, but usually (not always) by less than a factor of
two. But, the former is
is more complicated. With BigInt's, Numbers.jl code is faster, no matter
the value of the argument. Don't know why.
So, I use his routine for BigInts:
```
In[10]:= Timing[NextPrime[10^100];]
Out[10]= {0.016001, Null} (* consistently gives approx this time *)
```
```julia
julia> @time @bigint(nextprime(10^100)); # Hans' code is faster than
Mma for this value
elapsed time: 0.001856079 seconds (11624 bytes allocated)
```
### generate primes
Here are typical numbers:
```julia
julia> @time Numbers.primesieve(10^8);
elapsed time: 1.352471761 seconds (452346720 bytes allocated, 1.70% gc time)
julia> @time Base.primes(10^8);
elapsed time: 4.06491833 seconds (58631152 bytes allocated)
julia> @time PrimeSieve.genprimes(10^8);
elapsed time: 0.081639303 seconds (46091944 bytes allocated)
```
Numbers.jl is faster than Base, but allocates a lot.
genprimes(a,b) is much faster than Han's code. I can't find a
recommended efficient way to generate primes in Mma. I tried a few
suggestions; they are much slower. I have two algorithms for genprimes.
Which is better depends on a and b. The logic for choosing could be
greatly improved.
### note
libprimesieve can usually be interrupted without a crash, but may leave
Julia in
a bad state. libprimecount always segfaults on interrupt. I could do
disable_sigint,
but, I think it also has a memory leak, so better to allow stopping it.
Not ready for prime-time.
(I'll be here all week.)
--John
