Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Incorporate history from base #12

Merged
merged 35 commits into from
Jun 7, 2016
Merged
Changes from 1 commit
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
cd273dc
base/primes.jl: faster primes(n) and isprime(n) functions.
StefanKarpinski May 6, 2013
d6ef64f
isprime(n): better witness choices for various n.
StefanKarpinski May 6, 2013
3e4015d
factor: don't generate as many primes when doing trial division.
StefanKarpinski May 6, 2013
2531fd7
primes: implement and test prime stuff for BigInts [closes #3033]
StefanKarpinski May 6, 2013
62c36c5
factor: precompute and use a list of small primes as before.
StefanKarpinski May 6, 2013
0301c67
add isqrt (integer sqrt), use it to improve factor() (#3035)
JeffBezanson May 7, 2013
6f33c01
isprime: more precise witness choices for Miller-Rabin testing.
StefanKarpinski May 7, 2013
2367182
isprime: fix potential mulmod overflow exhibited by isprime(Uint64)
StefanKarpinski May 7, 2013
7313f43
isprime: support 128-bit integers via either demontion or promotion.
StefanKarpinski May 8, 2013
4f64dcc
isprime: fix isprime(typemin(Int128)+17)
StefanKarpinski May 8, 2013
173e1bc
isprime: fix precedence bug.
StefanKarpinski Jun 4, 2013
2f43407
primes: add comments identifying algorithms and sources.
StefanKarpinski Aug 1, 2013
ffd283c
removing function names from error messages
davidssmith Dec 4, 2013
87fc481
fix #5210, factor() failing for some integer types
JeffBezanson Dec 24, 2013
dd5add8
base/primes.jl: generalize primes/mask to BitArray or Array{Bool}.
StefanKarpinski Feb 16, 2014
7e5a179
deprecate old dict syntax and change uses in base and tests
JeffBezanson Oct 6, 2014
90938de
rename Uint => UInt (closes #8905)
StefanKarpinski Nov 5, 2014
076f520
itrunc -> trunc, etc, improve iround
simonbyrne Nov 23, 2014
ece0b08
in factor() check if remaining number is prime before looking for mor…
jlapeyre Jan 5, 2015
0cda3a0
Try and remove as many generic error() messages from base as possible
jakebolewski Jan 14, 2015
6e574fb
begin removing lowercase conversion functions (#1470)
JeffBezanson Mar 8, 2015
61b5b75
adds_license_headers [skip ci]
May 1, 2015
0d2243c
Speed up factor(n) + minor tweaks in other "primes.jl" functions
aiorla May 7, 2015
3ecf94b
Revert to previous header of primes(n) + 2 tests for big factorizations.
aiorla May 7, 2015
8f684cf
Redundant oftypes removes + Use of a more eficient searchsorted.
aiorla May 26, 2015
f28ef11
Fixing types conflict in pollardfactors!
aiorla May 28, 2015
8f48a46
Skipping the typemax check for BigInts
aiorla May 28, 2015
7756ed3
replace Union( ) with Union{ } everywhere
JeffBezanson Jun 15, 2015
efb8233
Replace http links with https links to suppress warning during doc li…
yuyichao Jul 27, 2015
c85ba8f
Use a wheel sieve to improve performance and export primesmask.
pabloferz Jul 6, 2015
0725830
Simplify primes and primesmask [fixes #14652]
pabloferz Jan 12, 2016
8324752
Replace unsafe_(get|set)index with inbounds macro
mbauman Feb 6, 2016
dcdf586
eachindex: use separate indexes for each array
timholy Mar 10, 2016
75582d9
rename base/primes.jl to src/Primes.jl
rfourquet May 29, 2016
a4de93f
update src/Primes.jl to master's version
rfourquet May 29, 2016
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 76 additions & 8 deletions src/Primes.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,13 @@
# This file is a part of Julia. License is MIT: http://julialang.org/license
__precompile__()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should keep this, just replace "is a part" with "was formerly a part"

module Primes

if VERSION >= v"0.5.0-dev+4340"
if isdefined(Base,:isprime)
import Base: isprime, primes, primesmask, factor
else
export isprime, primes, primesmask, factor
end
end

# Primes generating functions
# https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
@@ -9,8 +18,14 @@ const wheel = [4, 2, 4, 2, 4, 6, 2, 6]
const wheel_primes = [7, 11, 13, 17, 19, 23, 29, 31]
const wheel_indices = [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7]

@inline wheel_index(n) = ( (d, r) = divrem(n - 1, 30); return 8d + wheel_indices[r+2] )
@inline wheel_prime(n) = ( (d, r) = ((n - 1) >>> 3, (n - 1) & 7); return 30d + wheel_primes[r+1] )
@inline function wheel_index(n)
(d, r) = divrem(n - 1, 30)
8d + wheel_indices[r+2]
end
@inline function wheel_prime(n)
(d, r) = ((n - 1) >>> 3, (n - 1) & 7)
30d + wheel_primes[r+1]
end

function _primesmask(limit::Int)
limit < 7 && throw(ArgumentError("limit must be at least 7, got $limit"))
@@ -53,7 +68,12 @@ function _primesmask(lo::Int, hi::Int)
return sieve
end

# Sieve of the primes from lo up to hi represented as an array of booleans
"""
primesmask([lo,] hi)

Returns a prime sieve, as a `BitArray`, of the positive integers (from `lo`, if specified)
up to `hi`. Useful when working with either primes or composite numbers.
"""
function primesmask(lo::Int, hi::Int)
0 < lo ≤ hi || throw(ArgumentError("the condition 0 < lo ≤ hi must be met"))
sieve = falses(hi - lo + 1)
@@ -76,14 +96,19 @@ primesmask(limit::Int) = primesmask(1, limit)
primesmask(n::Integer) = n <= typemax(Int) ? primesmask(Int(n)) :
throw(ArgumentError("requested number of primes must be ≤ $(typemax(Int)), got $n"))

"""
primes([lo,] hi)

Returns a collection of the prime numbers (from `lo`, if specified) up to `hi`.
"""
function primes(lo::Int, hi::Int)
lo ≤ hi || throw(ArgumentError("the condition lo ≤ hi must be met"))
list = Int[]
lo ≤ 2 ≤ hi && push!(list, 2)
lo ≤ 3 ≤ hi && push!(list, 3)
lo ≤ 5 ≤ hi && push!(list, 5)
hi < 7 && return list
sizehint!(list, floor(Int, hi / log(hi)))
sizehint!(list, 5 + floor(Int, hi / (log(hi) - 1.12) - lo / (log(max(lo,2)) - 1.12*(lo > 7))) ) # http://projecteuclid.org/euclid.rmjm/1181070157
sieve = _primesmask(max(7, lo), hi)
lwi = wheel_index(lo - 1)
@inbounds for i = 1:length(sieve) # don't use eachindex here
@@ -95,10 +120,20 @@ primes(n::Int) = primes(1, n)

const PRIMES = primes(2^16)

# Small precomputed primes + Miller-Rabin for primality testing:
# https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
#
"""
isprime(x::Integer) -> Bool

Returns `true` if `x` is prime, and `false` otherwise.

```jldoctest
julia> isprime(3)
true
```
"""
function isprime(n::Integer)
# Small precomputed primes + Miller-Rabin for primality testing:
# https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
#
(n < 3 || iseven(n)) && return n == 2
n <= 2^16 && return PRIMES[searchsortedlast(PRIMES,n)] == n
s = trailing_zeros(n-1)
@@ -116,6 +151,22 @@ function isprime(n::Integer)
return true
end

"""
isprime(x::BigInt, [reps = 25]) -> Bool

Probabilistic primality test. Returns `true` if `x` is prime with high probability (pseudoprime);
and `false` if `x` is composite (not prime). The false positive rate is about `0.25^reps`.
`reps = 25` is considered safe for cryptographic applications (Knuth, Seminumerical Algorithms).

```jldoctest
julia> isprime(big(3))
true
```
"""
isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr{BigInt}, Cint), &x, reps) > 0



# Miller-Rabin witness choices based on:
# http://mathoverflow.net/questions/101922/smallest-collection-of-bases-for-prime-testing-of-64-bit-numbers
# http://primes.utm.edu/prove/merged.html
@@ -142,6 +193,21 @@ isprime(n::Int128) = n < 2 ? false :
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
# http://maths-people.anu.edu.au/~brent/pub/pub051.html
#
"""
factor(n) -> Dict

Compute the prime factorization of an integer `n`. Returns a dictionary. The keys of the
dictionary correspond to the factors, and hence are of the same type as `n`. The value
associated with each key indicates the number of times the factor appears in the
factorization.

```jldoctest
julia> factor(100) # == 2*2*5*5
Dict{Int64,Int64} with 2 entries:
2 => 2
5 => 2
```
"""
function factor{T<:Integer}(n::T)
0 < n || throw(ArgumentError("number to be factored must be ≥ 0, got $n"))
h = Dict{T,Int}()
@@ -204,3 +270,5 @@ function pollardfactors!{T<:Integer,K<:Integer}(n::T, h::Dict{K,Int})
end
end
end

end # module