*Prime number and factoring functions*

This package provides functions related to prime numbers and factoring,
mostly for generating and counting primes and factoring integers. This package
uses some of the fastest open-source libraries for these functions.
The name `Primes`

would be better, but that might cause collisions.

Note: Dec 5, 2016. Bitrot rests for no man. Quite a bit still works, but quite a bit not.
Broken tests are commented out, so tests suite should pass. If you want to fix something and
make a PR, please do. The most challenging problem is the C wrapper to make `libsmsieve.so`

from a standalone program no longer works for all inputs.

See LICENSE.md for links to the authors of the tables and the libraries used in this package.

I am unaware of binaries of libprimesieve and libprimecount for Windows and OSX, so there is no easy installation for these platforms.

This is not a registered package, and it has a non-registered dependency. You can install it (at least on Unix) with

```
Pkg.clone("https://github.com/jlapeyre/DeepConvert.jl")
Pkg.clone("https://github.com/jlapeyre/PrimeSieve.jl")
Pkg.build("PrimeSieve")
```

Some functions in this package

`genprimes(a,b)`

generate array of primes between`a`

and`b`

, inclusive`genprimes(b)`

generate array of primes between 2 and`b`

`mfactor(n)`

factor integers up to about 100 decimal digits.`primepi(n)`

the prime counting function: number of primes < n`countprimes(a,b)`

number of primes between a and b`nextprime(n)`

,`prevprime(n)`

first prime greater (or smaller) than n`someprimes(n1,n2)`

iterator. all primes from n1 through n2`someprimes(n2)`

all primes from 2 through n2`allprimes(n)`

iterator. all primes >= n`allprimes()`

all primes`nthprime(n)`

the nth prime`nprimes(n,start)`

generate array of the first n primes > start`isprime(z)`

primality test for Gaussian integers`randprime(a,b)`

random prime between`a`

and`b`

This package uses the following tables and libraries.

[T. Oliveira's tables of the prime counting function] (http://www.ieeta.pt/~tos/primes.html)

Prime number sieve library [libprimesieve] (http://primesieve.org/) and prime counting function library libprimecount

Integer factoring libraries msieve and gmp-ecm

The tables are encoded in Int128. The native type of the sieve (libprimesieve) is
Uint64. The input/output type of the fastest primepi algorithm in libprimecount, the Deleglise Rivat
algorithm, is Int128. There is a risk of overflow when constructing and giving
arguments to functions in this package. The easiest way to avoid this
is to put arguments in quotes: eg `countprimes("10^19","10^19+100")`

.
Also available are `@bigint`

and `@int128`

from DeepConvert.

Most of the following functions are vectorized.

`genprimes(start,stop)`

Return an array of all primes `>= start`

and `<= stop`

.

`genprimes(stop)`

Return an array of all primes between 1 and `stop`

`genprimes(start,stop; alg = algorithm)`

Generate primes using a specified algorithm. The algorithm must be
one of `:auto`

(the default), `:sieve`

, or `:next`

. Which algorithm is
more efficient depends on the parameters. In general, `:sieve`

is
better for larger intervals, and `:next`

is better for larger values
of `start`

. The keyword `:sieve`

uses a very fast sieve (libprimesieve), and
`:next`

uses the function `nextprime`

.

If you exceed the upper limit for argument to the sieve, then `:next`

is chosen automatically.

```
julia> @bigint genprimes(10^20, 10^20+1000)
24-element Array{BigInt,1}:
100000000000000000039 ...
```

This could also have been written `genprimes(bi"10^20", bi"10^20+1000")`

Computes the prime counting function.

`primepi(x; alg = algorithm)`

The efficient algorithms (or methods) are :auto (the default), :dr, and :tabsieve. The default, :auto, tries to choose the faster between :dr and :tabsieve (see Notes below). The other algorithms are slower in all cases. They are: :legendre, :lehmer, :meissel, :lmo, :sieve. The algorithm :dr uses an efficient parallel Deleglise Rivat method. The algorithm :tabsieve uses a combination of tables and a sieve and is more efficient when x is not too much greater than a table entry. (Note: Below, 10^14+10^8 is not too much greater than 10^14.) For example

```
julia> @time primepi(10^14+10^10; alg = :tabsieve)
elapsed time: 6.622672664 seconds (216 bytes allocated)
3205251958942
julia> @time primepi(10^14+10^10; alg = :dr) # Deleglise Rivat is faster
elapsed time: 0.495413145 seconds (208 bytes allocated)
3205251958942
julia> @time primepi(10^14+10^8; alg = :dr)
elapsed time: 0.505796298 seconds (208 bytes allocated)
3204944853481
julia> @time primepi(10^14+10^8; alg = :tabsieve) # Table and sieve is faster
elapsed time: 0.08235147 seconds (216 bytes allocated)
3204944853481
```

Factor an integer.

Example:

```
julia> @time mfactor("2^251-1")
elapsed time: 29.709989827 seconds (13283880 bytes allocated)
Dict{Int128,Int64} with 5 entries:
12070396178249893039969681 => 1
178230287214063289511 => 1
61676882198695257501367 => 1
503 => 1
54217 => 1
```

`mfactor(n; ecm = true)`

uses ecm
to search for factors larger than 15 digits, rather than less than 15 digits.

`mfactor(n; deadline = m)`

aborts the factoring after `m`

minutes.

`mfactor(n; logfile = "filename.log")`

writes information to a log file.

`mfactor(n; info = true)`

prints log information to the terminal.

`mfactor( @bigint [ 2^100 + i for i in -5:5] ) # returns an array of factorizations.`

`mfactor`

calls libmsieve and libecm.

Count the number of primes (or
prime tuplets in an interval. This
looks up the largest value in the table that is smaller than the
requested value and computes the remaining values. Note that `primepi`

is
logically equivalent to countprimes with `start=1`

. For `start=1`

,
The function `primepi`

is often much faster than, and is never slower than `countprimes`

```
countprimes(stop) # count the number of primes less than or equal to stop
countprimes(start,stop) # count the number of primes >= start and <= stop
countprimes([start], stop, tuplet=n) # Count prime n-tuplets
countprimes(start, stop, alg = algorithm) # Count prime n-tuplets
```

The default value of start is 1. The optional keyword argument 'tuplet' may take values between 1 and 6, that is primes, through prime sextuplets. Tables are implemented only for 'tuplet' equal to one, that is for primes, but not tuplets.

The optional keyword argument alg may be one of :tabsieve (the default),
:next, :nexta, or :sieve (:sieve will always be slower than :tabsieve).
As above, `:tabsieve`

uses a combination of tables and a fast sieve.
:next and :nexta are two different variants of `next_prime`

.

Examples

```
countprimes(100,1000) # the number of primes x satisfying 100 <= x <= 1000
143
countprimes(100,tuplet=3) # the number of prime triplets between 1 and 100
8
countprimes(10,1000,tuplet=6) # the number of prime sextuplets between 100 and 1000
1
```

If you quote the arguments (either as an expression or a string), they will be converted to Int128. This prevents overflow.

```
countprimes("10^19+10^9")
234057667299198865
```

If you use BigInt's, then the method :nexta will be chosen automatically. For example

```
julia> @bigint countprimes(10^50, 10^50+1000)
7
```

`nextprime(n)`

returns the smallest prime greater than n.

`prevprime(n)`

returns the largest prime less than n.

`nextprime(n,k)`

returns the kth prime following n.

`prevprime(n,k)`

returns the kth prime preceeding n.

Several algorithms are used. Finding the optimal one (of the available) is partially automated. nextprime1 and prevprime1 use an alternate algorithm coded by H W Borcher.

`someprimes(n2)`

All primes n, 2 <= n <= n2

`someprimes(n1,n2)`

All primes n, n1 <= n <= n2

`allprimes(n1)`

All primes n, n > n1

`allprimes()`

All primes

For example, here is the primorial function defined using an iterator:

```
julia> primorial(n) = prod(someprimes(n))
julia> @bigint primorial(100)
2305567963945518424753102147331756070
```

Returns the nth prime using a fast algorithm from libprimecount. The argument is converted to Int64.

`nthprime(n; alg = :sieve)`

uses the older algorithm from
libprimesieve, which is much slower.

Return an array of the first `n`

primes `>= start`

.

Usage

`nprimes(n,[start=1])`

Usage

`scountprimes([start],stop, tuplets=1)`

Print all primes (or prime n-tuplets) that are `>= start`

and `<= stop`

Usage

`printprimes([start],stop, [tuplet=1])`

The default value of 'start' is 1. The optional keyword argument 'tuplet' may take values between 1 and 6.

The legendre sum or phi function

`legendre(x,a)`

The arguments are converted to Int64.

The offset logarithmic integral. The argument is converted to Int64.

The inverse Li function. The argument is converted to Int64.

`randprime(a,b)`

choose a random prime between `a`

and `b`

. All primes in the range
are equally likely to be chosen.

`randprime(b)`

choose a random prime between 2 and `b`

.

`randprime(a,b; lim=n)`

find a maximum of `n`

random composite numbers
before giving up.

`randprime(a,b,n1,n2,...)`

return a n1xn2x... array of random primes

Returns `true`

if the
Gaussian integer `z`

is prime.

Get, set the sieve size in kilobytes. (setting does not seem to work)
`sz`

must satisfy `1 <= sz <= 2048`

Usage

```
primesievesize()
primesievesize(sz)
```

Get, set the number of threads used in the parallel sieve. By default, the number of cores is used.

Usage

```
primesieve_num_threads()
primesieve_num_threads(numthreads)
```

Get, set the number of threads used in the parallel primepi. By default, the number of cores is used.

Usage

```
primepi_num_threads()
primepi_num_threads(numthreads)
```

Return the largest value (as a `Uint64`

) that can be passed as the parameter
stop in the sieve.

Usage

`primemaxstop()`

Function that returns the largest allowed argument to `primepi`

when using the :dr algorithm.

Run a test of the sieve algorithm.

Usage

`primetest()`

Run a test of the primepi algorithms

Some of the above functions use tables, but this is completley hidden from the user. But it is also possible to access them directly. The tables work like this:

```
julia> @time countprimes(10^17 + 10^14 + 10^10)
elapsed time: 3.729049749 seconds (168 bytes allocated)
2626112053757377
```

To see what happened, we can look in the tables:

```
julia> primelookup(10^17 + 10^14 + 10^10)
(14,(2626111798288135,100100000000000000,10000000000))
```

The 14th table was used. The value of prime pi for `10^17+10^14`

,
`2626111798288135`

is in the table, and the primes in an
interval of length `10^10`

must be found with the sieves.

Print information about the prime pi tables.

```
julia> primetableinfo()
Tables of π(x). Listed are: table number, increment in x (and first value of x),
number of entries in the table, largest x in table.
table incr tab len max x
1 10^1 10^4 10^5
2 10^2 10^4 10^6
3 10^3 10^4 10^7
4 10^4 10^4 10^8
5 10^5 10^4 10^9
6 10^6 10^4 10^10
7 10^7 10^4 10^11
8 10^8 10^4 10^12
9 10^9 10^4 10^13
10 10^10 10^4 10^14
11 10^11 10^4 10^15
12 10^12 10^4 10^16
13 10^13 10^4 10^17
14 10^14 10^4 10^18
15 10^15 10^4 10^19
16 10^16 10^4 10^20
17 10^17 10^3 10^20
18 10^18 10^2 10^20
19 10^19 10^2 10^21
20 10^20 10^2 10^22
21 10^21 10^2 10^23
22 10^22 10^1 10^23
```

Look up a value of the prime pi function in the tables. This is only provided
to aid in understanding the behavior of `countprimes`

.

Usage

`primelookup(x)`

A tuple of a single element and another tuple of three elements is returned:

`(j,(p,y,rem))`

`j`

is the number of the best table found`y`

is the largest index satisfying`y<x`

found.`p`

is the value of prime pi at`y`

`rem`

is`x-y`

The array of type `Array{PrimeTable,1}`

containing the prime tables.
See tables.jl for the format.

Example

`show(map(length,primetables)) # see the number of tables and their lengths`

Function returning the path to the file containing the prime pi tables. The tables are loaded when the package is loaded.

For `x>typemax(Int)`

, you need to explicitly ask for a bigger data type.
For instance,

```
julia> countprimes(int128(10)^23)
1925320391606803968923
```

This example returned a value from the table. The argument was larger than than primemaxstop().

With any of the routines, you can quote the arguments and they will be converted to the appropriate type.

```
julia> countprimes(:(10^23))
1925320391606803968923
julia> countprimes("10^19 + 10^9")
234057667299198865
```

Routines that use the tables will convert the arguments to Int128. This is because
some indices in the tables are greater than `typemax(Uint64)`

. Routines that
only use the sieve will be converted to `Uint64`

, which is the data type that
the sieve routines use.

The largest stop value that may be given is `2^64 - 10 * 2^32`

.
The largest start value that may be given is `2^64 - 11 * 2^32`

.
The sieve works with the `Uint64`

data type. But conversions are done depending
on the types of start, stop, and n.

`countprimes`

returns `Int128`

, because it uses tables and sieves
The other routines only support smaller data types.

Return data type of tables. This should be Int128. The largest values cannot be used together with the sieve.

Return the native prime sieve type. This should be Uint64. libprimesieve
returns the data in various integer formats. These are chosen by the Julia
interface by the type of the `start`

parameter.

Return element type of values in table.

This is a only curiosity. It is supposed to be an optimized C (C++) function, but doing the same thing in a Julia loop is much faster.

`apopcount(arr)`

gives the number of 1's in the binary representation of
the array `arr`

. The length of the array is truncated to a multiple of 8.

Note that this treats the contents of the array as a bits type. In
particular, if `arr`

is not an array of bits type, then the number of 1's
in the pointers in the array are counted.

In the comments, Kim Walisch says this about the algorithm,

This algorithm counts the number of 1 bits (population count) in an array using 64-bit tree merging. To the best of my knowledge this is the fastest integer arithmetic bit population count algorithm, it uses only 8 operations for 8 bytes on 64-bit CPUs

Example:

```
julia> aa = Array(Uint64,100000000);
julia> fill!(aa,typemax(Uint64));
julia> @time apopcount(aa)
elapsed time: 0.140950344 seconds (96 bytes allocated) # test shows overhead is negligible
6400000000
julia> lpopcount(x) = ( s = 0; for i in 1:length(x) s += count_ones(x[i]) end; s)
julia> @time lpopcount(aa)
elapsed time: 0.106306266 seconds (96 bytes allocated) # Julia version is faster !
6400000000
```

The algorithm is by Cédric Lauradoux

The algorithms used by `:auto`

in `genprimes`

and `primepi`

are
simple and do not always choose the best method. They could be improved.
However, the routines are still fast. In the worst case, `genprimes`

is more than an order of magnitude faster than `Base.primes`

.

Interrupting a call to the sieves usually does not cause a memory error. But, libprimesieve apparently has some static state, such that, after the interrupt, subsequent sieving runs much slower, and may not give correct results.

Interrupting a call to libprimecount, results in a segfault. We could
`disable_sigint`

, but there appear to be memory leaks in libprimecount.
Better to crash Julia than the whole system.

msieve interrupt handler is installed, then all c code uses it. We need to install and then uninstall it on every call.