PrimeSieve

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.