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

Data types

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.

Functions

Most of the following functions are vectorized.

genprimes

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")

primepi

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

mfactor

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.

countprimes

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, prevprime

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.

Iterators

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

nthprime()

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.

nprimes

Return an array of the first n primes >= start.

Usage

nprimes(n,[start=1])

single threaded versions

Usage

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

printprimes

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.

legendrephi

The legendre sum or phi function

legendre(x,a)

The arguments are converted to Int64.

primeLi

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