Prime numbers again

[John Lapeyre]

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

[Hans W Borchers]

John,

please consider that my Numbers.jl package is something I did as my first trial 

in working with Julia -- and left it behind after a few days. It is certainly 

not in a state nor has the intention to be comparable in speed or completeness 

to other such systems.

As I said in another thread, when number theoretic computations are involved, 

PARI/GP is another system to look at closely, not only Mathematica or Maple. 

For instance, returning the next prime for really large input values takes only 

milliseconds ("See Ma, no tables!"):

    gp> nextprime(10^100)

    time = 9 ms.

    %1 = 100000000000000000000000000000000000000000000000000\

          00000000000000000000000000000000000000000000000267

There is no prime counting function. Probably as there is no direct need for 

this functionality for mathematicians, as you have already suspected. So the 

following trick counts primes:

    gp> n = 0;

    gp> forprime (p = 1841378967856, 1850000000000, n++);

    time = 1min, 49,432 ms.

    gp> n

    %4 = 305232179

Of course, a prime counting function based on tables cannot be beaten. That is 

why I'm quite glad you have provided this library for Julia.

If someone is thinking about a more extended number theory module for Julia, I 

would be glad to join and be of help.

[lapeyre....@gmail.com]

Hi Hans,

These things are learning exercises for me too. I agree that PARI/GP (and Gap)
are interesting. I copied some algorithms from PARI for PermutationsA.jl.

Still, I find it interesting that the nextprime functions from (your) Numbers.jl, Mathematica,
and PARI are very close in speed for some large arguments. I tried nextprime(10^1000)
in Numbers.jl and PARI/gp and they both took about 900ms (PARI was 10% faster). No
tables used anywhere. Mathematica took about 500ms for the same calculation.

I don't know how to count primes in an interval in PARI, but it does have a primepi(x)
function. The doc says  "Uses checkpointing and a naive O(x) algorithm".  If you look, this means
it finds the closest value in a table and then uses, in essence, nextprime. Maybe you are
right and they don't feel much need for something faster:

? primepi(10^9 + 10^8)       # looks like this is getting a close table value
time = 273 ms.
%29 = 55662470
? primepi(2*10^9)         # looks like no close table value. But I think the tables are configurable somehow
time = 2,407 ms.
%30 = 98222287

But, PrimeSieve.jl does not need to use tables:

julia> @time primepi(10^9 + 10^8 + 10^6)  # It has rather dense tables in this range of values. This is a table value.
elapsed time: 6.89e-6 seconds (168 bytes allocated)
55710422

julia> @time primepi(10^9 + 10^8 + 10^6; alg = :sieve)    # This uses no table, just the sieve.
elapsed time: 0.0700969 seconds (224 bytes allocated)
55710422

julia> @time primepi(10^9 + 10^8 + 10^6; alg = :dr)   # This uses no table, just a modern version of the Legendre formula
elapsed time: 0.002998701 seconds (224 bytes allocated)
55710422

Anyway, thanks for putting your code up. It's been useful for comparing and understanding sieves and nextprime.

--John