From cd273dc5cf25d2251928033b4ac801b506ab0345 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Mon, 6 May 2013 00:02:45 -0400 Subject: [PATCH 01/35] base/primes.jl: faster primes(n) and isprime(n) functions. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Sieve of Atkin for generating primes: http://en.wikipedia.org/wiki/Sieve_of_Atkin Code very loosely based on this: http://thomasinterestingblog.wordpress.com/2011/11/30/generating-primes-with-the-sieve-of-atkin-in-c/ http://dl.dropboxusercontent.com/u/29023244/atkin.cpp Miller-Rabin for primality testing: http://en.wikipedia.org/wiki/Miller–Rabin_primality_test # Conflicts: # base/intfuncs.jl # base/sysimg.jl # test/numbers.jl --- base/primes.jl | 84 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 base/primes.jl diff --git a/base/primes.jl b/base/primes.jl new file mode 100644 index 0000000..d492755 --- /dev/null +++ b/base/primes.jl @@ -0,0 +1,84 @@ +function primesmask(n::Int) + s = falses(n) + n < 2 && return s; s[2] = true + n < 3 && return s; s[3] = true + r = ifloor(sqrt(n)) + for x = 1:r + xx = x*x + for y = 1:r + yy = y*y + i, j, k = 4xx+yy, 3xx+yy, 3xx-yy + i <= n && (s[i] $= (i%12==1)|(i%12==5)) + j <= n && (s[j] $= (j%12==7)) + 1 <= k <= n && (s[k] $= (x>y)&(k%12==11)) + end + end + for i = 5:r + s[i] && (s[i*i:i*i:n] = false) + end + return s +end +primes(n::Int) = find(primesmask(n)) + +function isprime(n::Int) + n == 2 && return true + n <= 2 | iseven(n) && return false + s = trailing_zeros(n-1) + d = (n-1) >>> s + for a in (2,3,5,7,11,13,17) # TODO: more checks for larger n + a < n || break + x = powermod(a,d,n) + x == 1 && continue + t = s + while x != n-1 + (t-=1) <= 0 && return false + x = x*x % n + x == 1 && return false + end + end + return true +end + +# TODO: replace this factorization routine + +function factor{T<:Integer}(n::T) + if n <= 0 + error("factor: number to be factored must be positive") + end + h = (T=>Int)[] + if n == 1 return h end + local p::T + s = ifloor(sqrt(n)) + P = primes(n) + for p in P + if p > s + break + end + if n % p == 0 + while n % p == 0 + h[p] = get(h,p,0)+1 + n = div(n,p) + end + if n == 1 + return h + end + s = ifloor(sqrt(n)) + end + end + p = P[end]+2 + while p <= s + if n % p == 0 + while n % p == 0 + h[p] = get(h,p,0)+1 + n = div(n,p) + end + if n == 1 + return h + end + s = ifloor(sqrt(n)) + end + p += 2 + end + h[n] = 1 + return h +end From d6ef64f803a08731cb8586ab853abbb44d7791d5 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Mon, 6 May 2013 00:47:47 -0400 Subject: [PATCH 02/35] isprime(n): better witness choices for various n. http://mathoverflow.net/questions/101922/smallest-collection-of-bases-for-prime-testing-of-64-bit-numbers http://primes.utm.edu/prove/merged.html http://miller-rabin.appspot.com # Conflicts: # test/numbers.jl --- base/primes.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index d492755..bcafadb 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -20,12 +20,12 @@ function primesmask(n::Int) end primes(n::Int) = find(primesmask(n)) -function isprime(n::Int) +function isprime(n::Integer) n == 2 && return true n <= 2 | iseven(n) && return false s = trailing_zeros(n-1) d = (n-1) >>> s - for a in (2,3,5,7,11,13,17) # TODO: more checks for larger n + for a in witnesses(n) a < n || break x = powermod(a,d,n) x == 1 && continue @@ -38,6 +38,14 @@ function isprime(n::Int) end return true end +@eval begin + witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? + $(map(int32,(2,3))) : + $(map(int32,(2,7,61))) + witnesses(n::Union(Uint64,Int64)) = n < 341550071728321 ? + $(map(int64,(2,3,5,7,11,13,17))) : + $(map(int64,(2,325,9375,28178,450775,9780504,1795265022))) +end # TODO: replace this factorization routine From 3e4015d77c066c423357118682c057cc422bf0fc Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Mon, 6 May 2013 15:38:32 -0400 Subject: [PATCH 03/35] factor: don't generate as many primes when doing trial division. --- base/primes.jl | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index bcafadb..9966749 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -50,26 +50,21 @@ end # TODO: replace this factorization routine function factor{T<:Integer}(n::T) - if n <= 0 - error("factor: number to be factored must be positive") - end + 0 < n || error("factor: number to be factored must be positive") h = (T=>Int)[] - if n == 1 return h end - local p::T + n == 1 && return h + n <= 3 && (h[n] = 1; return h) + local s::T, p::T s = ifloor(sqrt(n)) - P = primes(n) + P = primes(s) for p in P - if p > s - break - end + p <= s || break if n % p == 0 while n % p == 0 h[p] = get(h,p,0)+1 n = div(n,p) end - if n == 1 - return h - end + n == 1 && return h s = ifloor(sqrt(n)) end end From 2531fd77b8443308ac81d0db510af73a3693f02c Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Mon, 6 May 2013 16:03:14 -0400 Subject: [PATCH 04/35] primes: implement and test prime stuff for BigInts [closes #3033] It is not recommended to use this code to factor large integers. The primality testing for BigInts calls GMP, so should be pretty good. # Conflicts: # base/gmp.jl # test/numbers.jl --- base/primes.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index 9966749..c4ed42c 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -18,7 +18,12 @@ function primesmask(n::Int) end return s end -primes(n::Int) = find(primesmask(n)) +function primesmask(n::Integer) + n <= typemax(Int) || error("primesmask: you want WAY too many primes ($n)") + primesmask(int(n)) +end + +primes(n::Integer) = find(primesmask(n)) function isprime(n::Integer) n == 2 && return true From 62c36c5d0602446bc8e38a42cd06616df450ca0b Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Mon, 6 May 2013 17:06:30 -0400 Subject: [PATCH 05/35] factor: precompute and use a list of small primes as before. # Conflicts: # base/sysimg.jl --- base/primes.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index c4ed42c..f56ea61 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -54,6 +54,8 @@ end # TODO: replace this factorization routine +const PRIMES = primes(10000) + function factor{T<:Integer}(n::T) 0 < n || error("factor: number to be factored must be positive") h = (T=>Int)[] @@ -61,8 +63,7 @@ function factor{T<:Integer}(n::T) n <= 3 && (h[n] = 1; return h) local s::T, p::T s = ifloor(sqrt(n)) - P = primes(s) - for p in P + for p in PRIMES p <= s || break if n % p == 0 while n % p == 0 @@ -73,7 +74,7 @@ function factor{T<:Integer}(n::T) s = ifloor(sqrt(n)) end end - p = P[end]+2 + p = PRIMES[end]+2 while p <= s if n % p == 0 while n % p == 0 From 0301c67e7532350e68bbb42f8e64d2826e0910f8 Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Tue, 7 May 2013 00:49:11 -0400 Subject: [PATCH 06/35] add isqrt (integer sqrt), use it to improve factor() (#3035) export primes # Conflicts: # base/exports.jl # base/gmp.jl # base/intfuncs.jl # base/mpfr.jl --- base/primes.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index f56ea61..f5f5eef 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -62,7 +62,7 @@ function factor{T<:Integer}(n::T) n == 1 && return h n <= 3 && (h[n] = 1; return h) local s::T, p::T - s = ifloor(sqrt(n)) + s = isqrt(n) for p in PRIMES p <= s || break if n % p == 0 @@ -71,7 +71,7 @@ function factor{T<:Integer}(n::T) n = div(n,p) end n == 1 && return h - s = ifloor(sqrt(n)) + s = isqrt(n) end end p = PRIMES[end]+2 @@ -84,7 +84,7 @@ function factor{T<:Integer}(n::T) if n == 1 return h end - s = ifloor(sqrt(n)) + s = isqrt(n) end p += 2 end From 6f33c0138dfe929f8a39c8c956a4be2544e6a47c Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Tue, 7 May 2013 17:29:18 -0400 Subject: [PATCH 07/35] isprime: more precise witness choices for Miller-Rabin testing. --- base/primes.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index f5f5eef..d395099 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -43,14 +43,14 @@ function isprime(n::Integer) end return true end -@eval begin - witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? - $(map(int32,(2,3))) : - $(map(int32,(2,7,61))) - witnesses(n::Union(Uint64,Int64)) = n < 341550071728321 ? - $(map(int64,(2,3,5,7,11,13,17))) : - $(map(int64,(2,325,9375,28178,450775,9780504,1795265022))) -end +witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3) +witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61) +witnesses(n::Union(Uint64,Int64)) = + n < 1373653 ? (2,3) : + n < 4759123141 ? (2,7,61) : + n < 2152302898747 ? (2,3,5,7,11) : + n < 3474749660383 ? (2,3,5,7,11,13) : + (2,325,9375,28178,450775,9780504,1795265022) # TODO: replace this factorization routine From 2367182def53df4a8b82923ce0b910404e301903 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Tue, 7 May 2013 18:16:21 -0400 Subject: [PATCH 08/35] isprime: fix potential mulmod overflow exhibited by isprime(Uint64) # Conflicts: # test/numbers.jl --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index d395099..7c20486 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -37,7 +37,7 @@ function isprime(n::Integer) t = s while x != n-1 (t-=1) <= 0 && return false - x = x*x % n + x = oftype(n, widemul(x,x) % n) x == 1 && return false end end From 7313f4331453cf6b81b2d5d6b75cd4eefc950d91 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Tue, 7 May 2013 22:49:31 -0400 Subject: [PATCH 09/35] isprime: support 128-bit integers via either demontion or promotion. --- base/primes.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/base/primes.jl b/base/primes.jl index 7c20486..6c689c5 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -52,6 +52,9 @@ witnesses(n::Union(Uint64,Int64)) = n < 3474749660383 ? (2,3,5,7,11,13) : (2,325,9375,28178,450775,9780504,1795265022) +isprime(n::Uint128) = n <= typemax(Uint64) ? isprime(uint64(n)) : isprime(BigInt(n)) +isprime(n::Int128) = n <= typemax(Int64) ? isprime(int64(n)) : isprime(BigInt(n)) + # TODO: replace this factorization routine const PRIMES = primes(10000) From 4f64dccfd0e53ffda903bb3f478a4cbcd2afd0cb Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Wed, 8 May 2013 00:29:04 -0400 Subject: [PATCH 10/35] isprime: fix isprime(typemin(Int128)+17) --- base/primes.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 6c689c5..e07a912 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -52,8 +52,10 @@ witnesses(n::Union(Uint64,Int64)) = n < 3474749660383 ? (2,3,5,7,11,13) : (2,325,9375,28178,450775,9780504,1795265022) -isprime(n::Uint128) = n <= typemax(Uint64) ? isprime(uint64(n)) : isprime(BigInt(n)) -isprime(n::Int128) = n <= typemax(Int64) ? isprime(int64(n)) : isprime(BigInt(n)) +isprime(n::Uint128) = + n <= typemax(Uint64) ? isprime(uint64(n)) : isprime(BigInt(n)) +isprime(n::Int128) = n < 2 ? false : + n <= typemax(Int64) ? isprime(int64(n)) : isprime(BigInt(n)) # TODO: replace this factorization routine From 173e1bc7df91258df03ab84421396d90a3a96d17 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Tue, 4 Jun 2013 15:12:04 -0400 Subject: [PATCH 11/35] isprime: fix precedence bug. --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index e07a912..24b73c3 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -27,7 +27,7 @@ primes(n::Integer) = find(primesmask(n)) function isprime(n::Integer) n == 2 && return true - n <= 2 | iseven(n) && return false + (n < 2) | iseven(n) && return false s = trailing_zeros(n-1) d = (n-1) >>> s for a in witnesses(n) From 2f43407818cdafbf92367432b53982f1cdc77a98 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Thu, 1 Aug 2013 14:38:17 -0400 Subject: [PATCH 12/35] primes: add comments identifying algorithms and sources. --- base/primes.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index 24b73c3..873f842 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -1,3 +1,9 @@ +# Sieve of Atkin for generating primes: +# http://en.wikipedia.org/wiki/Sieve_of_Atkin +# Code very loosely based on this: +# http://thomasinterestingblog.wordpress.com/2011/11/30/generating-primes-with-the-sieve-of-atkin-in-c/ +# http://dl.dropboxusercontent.com/u/29023244/atkin.cpp +# function primesmask(n::Int) s = falses(n) n < 2 && return s; s[2] = true @@ -25,6 +31,9 @@ end primes(n::Integer) = find(primesmask(n)) +# Miller-Rabin for primality testing: +# http://en.wikipedia.org/wiki/Miller–Rabin_primality_test +# function isprime(n::Integer) n == 2 && return true (n < 2) | iseven(n) && return false @@ -43,6 +52,12 @@ function isprime(n::Integer) end return true end + +# 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 +# http://miller-rabin.appspot.com +# witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3) witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61) witnesses(n::Union(Uint64,Int64)) = @@ -57,7 +72,7 @@ isprime(n::Uint128) = isprime(n::Int128) = n < 2 ? false : n <= typemax(Int64) ? isprime(int64(n)) : isprime(BigInt(n)) -# TODO: replace this factorization routine +# TODO: faster factorization algorithms? const PRIMES = primes(10000) From ffd283cf5e5b0a1d558e37ab55eea0b76bf5d1e5 Mon Sep 17 00:00:00 2001 From: davidssmith Date: Wed, 4 Dec 2013 14:45:17 -0600 Subject: [PATCH 13/35] removing function names from error messages # Conflicts: # base/pcre.jl # base/reduce.jl # base/reflection.jl # base/statistics.jl # base/util.jl --- base/primes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 873f842..26a1ed7 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -25,7 +25,7 @@ function primesmask(n::Int) return s end function primesmask(n::Integer) - n <= typemax(Int) || error("primesmask: you want WAY too many primes ($n)") + n <= typemax(Int) || error("you want WAY too many primes ($n)") primesmask(int(n)) end @@ -77,7 +77,7 @@ isprime(n::Int128) = n < 2 ? false : const PRIMES = primes(10000) function factor{T<:Integer}(n::T) - 0 < n || error("factor: number to be factored must be positive") + 0 < n || error("number to be factored must be positive") h = (T=>Int)[] n == 1 && return h n <= 3 && (h[n] = 1; return h) From 87fc481b40a266093a4aa0e7fcc58a1755ae5621 Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Tue, 24 Dec 2013 01:26:02 -0500 Subject: [PATCH 14/35] fix #5210, factor() failing for some integer types # Conflicts: # test/numbers.jl --- base/primes.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 26a1ed7..21c4df1 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -84,11 +84,14 @@ function factor{T<:Integer}(n::T) local s::T, p::T s = isqrt(n) for p in PRIMES - p <= s || break + if p > s + h[n] = 1 + return h + end if n % p == 0 while n % p == 0 h[p] = get(h,p,0)+1 - n = div(n,p) + n = oftype(n, div(n,p)) end n == 1 && return h s = isqrt(n) @@ -99,7 +102,7 @@ function factor{T<:Integer}(n::T) if n % p == 0 while n % p == 0 h[p] = get(h,p,0)+1 - n = div(n,p) + n = oftype(n, div(n,p)) end if n == 1 return h From dd5add8189334289107b1796f27bfe8cecadd684 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Sun, 16 Feb 2014 18:48:11 -0500 Subject: [PATCH 15/35] base/primes.jl: generalize primes/mask to BitArray or Array{Bool}. For smaller numbers of primes, the Sieve of Eratosthenes still seems to be quite a bit faster. Perhaps we ought to use a polyalgorithm here, only using the Sieve of Atkin for large N. --- base/primes.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 21c4df1..a8898ae 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -4,8 +4,8 @@ # http://thomasinterestingblog.wordpress.com/2011/11/30/generating-primes-with-the-sieve-of-atkin-in-c/ # http://dl.dropboxusercontent.com/u/29023244/atkin.cpp # -function primesmask(n::Int) - s = falses(n) +function primesmask(s::AbstractVector{Bool}) + n = length(s) n < 2 && return s; s[2] = true n < 3 && return s; s[3] = true r = ifloor(sqrt(n)) @@ -24,12 +24,11 @@ function primesmask(n::Int) end return s end -function primesmask(n::Integer) - n <= typemax(Int) || error("you want WAY too many primes ($n)") - primesmask(int(n)) -end +primesmask(n::Int) = primesmask(falses(n)) +primesmask(n::Integer) = n <= typemax(Int) ? primesmask(int(n)) : + error("you want WAY too many primes ($n)") -primes(n::Integer) = find(primesmask(n)) +primes(n::Union(Integer,AbstractVector{Bool})) = find(primesmask(n)) # Miller-Rabin for primality testing: # http://en.wikipedia.org/wiki/Miller–Rabin_primality_test From 7e5a17915aa2ffd3d137db52a4d1454a63d3077d Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Mon, 6 Oct 2014 13:40:28 -0400 Subject: [PATCH 16/35] deprecate old dict syntax and change uses in base and tests # Conflicts: # base/LineEdit.jl # base/REPL.jl # base/cartesian.jl # base/client.jl # base/combinatorics.jl # base/datafmt.jl # base/dates/adjusters.jl # base/dates/io.jl # base/dates/query.jl # base/dict.jl # base/latex_symbols.jl # base/loading.jl # base/multi.jl # base/pkg/generate.jl # base/pkg/github.jl # base/pkg/query.jl # base/pkg/resolve.jl # base/pkg/resolve/interface.jl # base/pkg/resolve/maxsum.jl # base/profile.jl # base/quadgk.jl # base/show.jl # src/jlfrontend.scm # src/julia-parser.scm # test/collections.jl # test/core.jl # test/dates/io.jl # test/dates/query.jl # test/hashing.jl # test/keywordargs.jl # test/lineedit.jl # test/repl.jl # test/resolve.jl # test/sparse.jl # test/spawn.jl # test/strings.jl --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index a8898ae..0d81c8e 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -77,7 +77,7 @@ const PRIMES = primes(10000) function factor{T<:Integer}(n::T) 0 < n || error("number to be factored must be positive") - h = (T=>Int)[] + h = Dict{T,Int}() n == 1 && return h n <= 3 && (h[n] = 1; return h) local s::T, p::T From 90938de2e6b8ab077c069c88af3edf87c78fa6b8 Mon Sep 17 00:00:00 2001 From: Stefan Karpinski Date: Wed, 5 Nov 2014 16:53:54 +0100 Subject: [PATCH 17/35] rename Uint => UInt (closes #8905) # Conflicts: # NEWS.md # base/REPL.jl # base/Terminals.jl # base/abstractarray.jl # base/array.jl # base/ascii.jl # base/base.jl # base/base64.jl # base/bitarray.jl # base/boot.jl # base/broadcast.jl # base/c.jl # base/char.jl # base/client.jl # base/combinatorics.jl # base/constants.jl # base/dSFMT.jl # base/datafmt.jl # base/deepcopy.jl # base/deprecated.jl # base/dict.jl # base/env.jl # base/error.jl # base/expr.jl # base/fftw.jl # base/file.jl # base/float.jl # base/float16.jl # base/floatfuncs.jl # base/fs.jl # base/gmp.jl # base/grisu.jl # base/grisu/bignum.jl # base/grisu/fastfixed.jl # base/grisu/fastprecision.jl # base/grisu/fastshortest.jl # base/grisu/float.jl # base/hashing.jl # base/hashing2.jl # base/int.jl # base/interactiveutil.jl # base/intfuncs.jl # base/intset.jl # base/io.jl # base/iobuffer.jl # base/iostream.jl # base/libc.jl # base/linalg/arpack.jl # base/linalg/bitarray.jl # base/linalg/blas.jl # base/linalg/cholmod.jl # base/linalg/lapack.jl # base/linalg/matmul.jl # base/loading.jl # base/math.jl # base/mmap.jl # base/mpfr.jl # base/multi.jl # base/multidimensional.jl # base/multimedia.jl # base/nullable.jl # base/operators.jl # base/osutils.jl # base/path.jl # base/pcre.jl # base/pkg/types.jl # base/pointer.jl # base/poll.jl # base/precompile.jl # base/printf.jl # base/process.jl # base/profile.jl # base/random.jl # base/range.jl # base/reduce.jl # base/reducedim.jl # base/reflection.jl # base/regex.jl # base/serialize.jl # base/sharedarray.jl # base/show.jl # base/socket.jl # base/stat.jl # base/statistics.jl # base/stream.jl # base/string.jl # base/sysinfo.jl # base/task.jl # base/tuple.jl # base/utf16.jl # base/utf32.jl # base/utf8.jl # base/utf8proc.jl # base/util.jl # base/version.jl # doc/helpdb.jl # doc/images/github_metadata_develbranch.png # doc/images/github_metadata_fork.png # doc/images/github_metadata_pullrequest.png # doc/images/travis-icon.png # doc/manual/calling-c-and-fortran-code.rst # doc/manual/conversion-and-promotion.rst # doc/manual/faq.rst # doc/manual/integers-and-floating-point-numbers.rst # doc/manual/mathematical-operations.rst # doc/manual/methods.rst # doc/manual/networking-and-streams.rst # doc/manual/strings.rst # doc/manual/style-guide.rst # doc/manual/types.rst # doc/stdlib/base.rst # examples/lru.jl # examples/lru_test.jl # examples/plife.jl # examples/typetree.jl # src/file_constants.h # src/init.c # src/julia-parser.scm # test/arrayops.jl # test/bigint.jl # test/bitarray.jl # test/ccall.jl # test/collections.jl # test/core.jl # test/dates/periods.jl # test/file.jl # test/grisu.jl # test/hashing.jl # test/iobuffer.jl # test/math.jl # test/mpfr.jl # test/netload/nettest.jl # test/nullable.jl # test/numbers.jl # test/parallel.jl # test/perf/kernel/actor_centrality.jl # test/perf/kernel/go_benchmark.jl # test/perf/kernel/ziggurat.jl # test/perf/micro/perf.jl # test/perf/shootout/fasta.jl # test/perf/shootout/mandelbrot.jl # test/perf/shootout/meteor_contest.jl # test/perf/shootout/revcomp.jl # test/pollfd.jl # test/random.jl # test/ranges.jl # test/reducedim.jl # test/repl.jl # test/show.jl # test/socket.jl # test/sorting.jl # test/sparse.jl # test/statistics.jl # test/strings.jl # test/unicode.jl --- base/primes.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 0d81c8e..c42b62d 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -57,17 +57,17 @@ end # http://primes.utm.edu/prove/merged.html # http://miller-rabin.appspot.com # -witnesses(n::Union(Uint8,Int8,Uint16,Int16)) = (2,3) -witnesses(n::Union(Uint32,Int32)) = n < 1373653 ? (2,3) : (2,7,61) -witnesses(n::Union(Uint64,Int64)) = +witnesses(n::Union(UInt8,Int8,UInt16,Int16)) = (2,3) +witnesses(n::Union(UInt32,Int32)) = n < 1373653 ? (2,3) : (2,7,61) +witnesses(n::Union(UInt64,Int64)) = n < 1373653 ? (2,3) : n < 4759123141 ? (2,7,61) : n < 2152302898747 ? (2,3,5,7,11) : n < 3474749660383 ? (2,3,5,7,11,13) : (2,325,9375,28178,450775,9780504,1795265022) -isprime(n::Uint128) = - n <= typemax(Uint64) ? isprime(uint64(n)) : isprime(BigInt(n)) +isprime(n::UInt128) = + n <= typemax(UInt64) ? isprime(uint64(n)) : isprime(BigInt(n)) isprime(n::Int128) = n < 2 ? false : n <= typemax(Int64) ? isprime(int64(n)) : isprime(BigInt(n)) From 076f5208c921b4f7781899addc3af232d5948c44 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Sun, 23 Nov 2014 22:14:07 +0000 Subject: [PATCH 18/35] itrunc -> trunc, etc, improve iround # Conflicts: # base/abstractarray.jl # base/bitarray.jl # base/char.jl # base/complex.jl # base/darray.jl # base/dates/ranges.jl # base/deprecated.jl # base/exports.jl # base/float.jl # base/float16.jl # base/floatfuncs.jl # base/grisu/float.jl # base/int.jl # base/interactiveutil.jl # base/intfuncs.jl # base/libc.jl # base/linalg/bidiag.jl # base/linalg/dense.jl # base/linalg/givens.jl # base/linalg/matmul.jl # base/linalg/tridiag.jl # base/mmap.jl # base/mpfr.jl # base/pkg/resolve/maxsum.jl # base/poll.jl # base/printf.jl # base/profile.jl # base/range.jl # base/rational.jl # base/sort.jl # base/sparse/sparsematrix.jl # base/special/gamma.jl # base/statistics.jl # doc/helpdb.jl # doc/manual/arrays.rst # doc/manual/mathematical-operations.rst # doc/stdlib/base.rst # src/intrinsics.cpp # test/arrayperf.jl # test/complex.jl # test/euler.jl # test/file.jl # test/linalg/triangular.jl # test/mpfr.jl # test/numbers.jl # test/random.jl # test/rounding.jl # test/sparse.jl # test/strings.jl --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index c42b62d..fc9bf4e 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -8,7 +8,7 @@ function primesmask(s::AbstractVector{Bool}) n = length(s) n < 2 && return s; s[2] = true n < 3 && return s; s[3] = true - r = ifloor(sqrt(n)) + r = floor(Int,sqrt(n)) for x = 1:r xx = x*x for y = 1:r From ece0b0858049026b70886494502136ab00093cd5 Mon Sep 17 00:00:00 2001 From: John Lapeyre Date: Mon, 5 Jan 2015 12:41:07 +0100 Subject: [PATCH 19/35] in factor() check if remaining number is prime before looking for more factors. # Conflicts: # test/numbers.jl --- base/primes.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index fc9bf4e..f043f43 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -93,6 +93,10 @@ function factor{T<:Integer}(n::T) n = oftype(n, div(n,p)) end n == 1 && return h + if isprime(n) + h[n] = 1 + return h + end s = isqrt(n) end end @@ -103,7 +107,9 @@ function factor{T<:Integer}(n::T) h[p] = get(h,p,0)+1 n = oftype(n, div(n,p)) end - if n == 1 + n == 1 && return h + if isprime(n) + h[n] = 1 return h end s = isqrt(n) From 0cda3a00ddb9b77ce861087727aebd39b323eb73 Mon Sep 17 00:00:00 2001 From: jake bolewski Date: Wed, 14 Jan 2015 14:35:47 -0500 Subject: [PATCH 20/35] Try and remove as many generic error() messages from base as possible * throw more specific Exception types * make error messages more consistent * give more context for the error when possible * update tests # Conflicts: # base/abstractarray.jl # base/array.jl # base/ascii.jl # base/base64.jl # base/bitarray.jl # base/broadcast.jl # base/cartesian.jl # base/client.jl # base/collections.jl # base/combinatorics.jl # base/darray.jl # base/datafmt.jl # base/dict.jl # base/dsp.jl # base/env.jl # base/file.jl # base/fs.jl # base/gmp.jl # base/help.jl # base/intfuncs.jl # base/intset.jl # base/io.jl # base/iobuffer.jl # base/iostream.jl # base/libc.jl # base/linalg/arnoldi.jl # base/linalg/arpack.jl # base/linalg/bidiag.jl # base/linalg/bitarray.jl # base/linalg/bunchkaufman.jl # base/linalg/dense.jl # base/linalg/diagonal.jl # base/linalg/generic.jl # base/linalg/givens.jl # base/linalg/lapack.jl # base/linalg/triangular.jl # base/linalg/tridiag.jl # base/loading.jl # base/math.jl # base/mmap.jl # base/mpfr.jl # base/multidimensional.jl # base/operators.jl # base/osutils.jl # base/path.jl # base/pcre.jl # base/pointer.jl # base/poll.jl # base/printf.jl # base/process.jl # base/profile.jl # base/random.jl # base/range.jl # base/rational.jl # base/reduce.jl # base/reflection.jl # base/regex.jl # base/rounding.jl # base/sharedarray.jl # base/show.jl # base/socket.jl # base/sort.jl # base/statistics.jl # base/stream.jl # base/string.jl # base/subarray.jl # base/tuple.jl # base/utf16.jl # base/utf32.jl # base/utf8.jl # base/version.jl # test/base64.jl # test/broadcast.jl # test/dates/ranges.jl # test/dict.jl # test/file.jl # test/iobuffer.jl # test/mod2pi.jl # test/reduce.jl # test/sets.jl # test/socket.jl # test/spawn.jl # test/statistics.jl --- base/primes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index f043f43..f7c6d32 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -26,7 +26,7 @@ function primesmask(s::AbstractVector{Bool}) end primesmask(n::Int) = primesmask(falses(n)) primesmask(n::Integer) = n <= typemax(Int) ? primesmask(int(n)) : - error("you want WAY too many primes ($n)") + throw(ArgumentError("requested number of primes must be ≤ $(typemax(Int)), got $n")) primes(n::Union(Integer,AbstractVector{Bool})) = find(primesmask(n)) @@ -76,7 +76,7 @@ isprime(n::Int128) = n < 2 ? false : const PRIMES = primes(10000) function factor{T<:Integer}(n::T) - 0 < n || error("number to be factored must be positive") + 0 < n || throw(ArgumentError("number to be factored must be ≥ 0, got $n")) h = Dict{T,Int}() n == 1 && return h n <= 3 && (h[n] = 1; return h) From 6e574fbb3000cd9a71e918230227d5604460acac Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Sun, 8 Mar 2015 16:55:18 -0400 Subject: [PATCH 21/35] begin removing lowercase conversion functions (#1470) not quite done with renaming [ci skip] # Conflicts: # base/LineEdit.jl # base/Makefile # base/Terminals.jl # base/abstractarray.jl # base/array.jl # base/ascii.jl # base/base.jl # base/base64.jl # base/bitarray.jl # base/bool.jl # base/broadcast.jl # base/char.jl # base/client.jl # base/combinatorics.jl # base/complex.jl # base/datafmt.jl # base/dates/conversions.jl # base/dates/ranges.jl # base/deprecated.jl # base/dict.jl # base/dsp.jl # base/env.jl # base/exports.jl # base/fastmath.jl # base/fftw.jl # base/file.jl # base/float.jl # base/float16.jl # base/fs.jl # base/gmp.jl # base/grisu/bignum.jl # base/grisu/fastfixed.jl # base/grisu/fastprecision.jl # base/grisu/fastshortest.jl # base/grisu/float.jl # base/hashing.jl # base/hashing2.jl # base/int.jl # base/intfuncs.jl # base/intset.jl # base/io.jl # base/iobuffer.jl # base/iostream.jl # base/latex_symbols.jl # base/libc.jl # base/linalg.jl # base/linalg/bitarray.jl # base/linalg/dense.jl # base/linalg/factorization.jl # base/linalg/lapack.jl # base/linalg/lu.jl # base/managers.jl # base/math.jl # base/mmap.jl # base/multi.jl # base/number.jl # base/operators.jl # base/path.jl # base/pcre.jl # base/pkg/generate.jl # base/pkg/github.jl # base/pkg/resolve/fieldvalue.jl # base/pkg/resolve/maxsum.jl # base/pkg/resolve/versionweight.jl # base/pointer.jl # base/printf.jl # base/process.jl # base/profile.jl # base/random.jl # base/range.jl # base/rational.jl # base/regex.jl # base/serialize.jl # base/sharedarray.jl # base/show.jl # base/socket.jl # base/sparse/cholmod.jl # base/sparse/cholmod_h.jl # base/sparse/sparsematrix.jl # base/sparse/spqr.jl # base/sparse/umfpack.jl # base/special/bessel.jl # base/special/gamma.jl # base/stat.jl # base/statistics.jl # base/stream.jl # base/string.jl # base/subarray.jl # base/sysinfo.jl # base/utf16.jl # base/utf32.jl # base/utf8.jl # base/utf8proc.jl # base/version.jl # test/arrayops.jl # test/bigint.jl # test/bitarray.jl # test/char.jl # test/combinatorics.jl # test/complex.jl # test/core.jl # test/dates/periods.jl # test/dates/types.jl # test/euler.jl # test/file.jl # test/grisu.jl # test/hashing.jl # test/intfuncs.jl # test/iobuffer.jl # test/llvmcall.jl # test/mod2pi.jl # test/mpfr.jl # test/netload/nettest.jl # test/numbers.jl # test/parallel.jl # test/perf/kernel/ziggurat.jl # test/perf/micro/perf.jl # test/perf/shootout/mandelbrot.jl # test/perf/shootout/meteor_contest.jl # test/perf/shootout/revcomp.jl # test/perf/sparse/getindex.jl # test/pollfd.jl # test/random.jl # test/ranges.jl # test/reduce.jl # test/serialize.jl # test/show.jl # test/socket.jl # test/strings.jl # test/version.jl --- base/primes.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index f7c6d32..769da3d 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -25,7 +25,7 @@ function primesmask(s::AbstractVector{Bool}) return s end primesmask(n::Int) = primesmask(falses(n)) -primesmask(n::Integer) = n <= typemax(Int) ? primesmask(int(n)) : +primesmask(n::Integer) = n <= typemax(Int) ? primesmask(Int(n)) : throw(ArgumentError("requested number of primes must be ≤ $(typemax(Int)), got $n")) primes(n::Union(Integer,AbstractVector{Bool})) = find(primesmask(n)) @@ -67,9 +67,9 @@ witnesses(n::Union(UInt64,Int64)) = (2,325,9375,28178,450775,9780504,1795265022) isprime(n::UInt128) = - n <= typemax(UInt64) ? isprime(uint64(n)) : isprime(BigInt(n)) + n <= typemax(UInt64) ? isprime(UInt64(n)) : isprime(BigInt(n)) isprime(n::Int128) = n < 2 ? false : - n <= typemax(Int64) ? isprime(int64(n)) : isprime(BigInt(n)) + n <= typemax(Int64) ? isprime(Int64(n)) : isprime(BigInt(n)) # TODO: faster factorization algorithms? From 61b5b750befa58e21edcce2230fb06ee36d0922e Mon Sep 17 00:00:00 2001 From: peter1000 Date: Fri, 1 May 2015 15:08:12 -0300 Subject: [PATCH 22/35] adds_license_headers [skip ci] JuliaLang#11073 (comment) JuliaLang#11023 Related pullrequests are: JuliaLang#11079 JuliaLang#11084 # Conflicts: # base/Dates.jl # base/Enums.jl # base/LineEdit.jl # base/REPL.jl # base/REPLCompletions.jl # base/Terminals.jl # base/abstractarray.jl # base/array.jl # base/ascii.jl # base/base.jl # base/base64.jl # base/basedocs.jl # base/bitarray.jl # base/bool.jl # base/boot.jl # base/broadcast.jl # base/build.h # base/c.jl # base/cartesian.jl # base/char.jl # base/client.jl # base/collections.jl # base/combinatorics.jl # base/complex.jl # base/constants.jl # base/dSFMT.jl # base/datafmt.jl # base/dates/accessors.jl # base/dates/adjusters.jl # base/dates/arithmetic.jl # base/dates/conversions.jl # base/dates/io.jl # base/dates/periods.jl # base/dates/query.jl # base/dates/ranges.jl # base/dates/types.jl # base/deepcopy.jl # base/deprecated.jl # base/dict.jl # base/docs.jl # base/dsp.jl # base/emoji_symbols.jl # base/env.jl # base/errno.jl # base/error.jl # base/exports.jl # base/expr.jl # base/fastmath.jl # base/fftw.jl # base/file.jl # base/float.jl # base/float16.jl # base/floatfuncs.jl # base/fs.jl # base/functors.jl # base/gmp.jl # base/grisu.jl # base/hashing.jl # base/hashing2.jl # base/help.jl # base/i18n.jl # base/inference.jl # base/int.jl # base/interactiveutil.jl # base/intfuncs.jl # base/intset.jl # base/io.jl # base/iobuffer.jl # base/iostream.jl # base/iterator.jl # base/latex_symbols.jl # base/libc.jl # base/libdl.jl # base/linalg.jl # base/linalg/arnoldi.jl # base/linalg/arpack.jl # base/linalg/bidiag.jl # base/linalg/bitarray.jl # base/linalg/blas.jl # base/linalg/bunchkaufman.jl # base/linalg/cholesky.jl # base/linalg/dense.jl # base/linalg/diagonal.jl # base/linalg/exceptions.jl # base/linalg/factorization.jl # base/linalg/generic.jl # base/linalg/givens.jl # base/linalg/lapack.jl # base/linalg/ldlt.jl # base/linalg/lu.jl # base/linalg/matmul.jl # base/linalg/rectfullpacked.jl # base/linalg/special.jl # base/linalg/symmetric.jl # base/linalg/triangular.jl # base/linalg/tridiag.jl # base/linalg/uniformscaling.jl # base/loading.jl # base/lock.jl # base/managers.jl # base/markdown/Common/Common.jl # base/markdown/Common/block.jl # base/markdown/Common/inline.jl # base/markdown/GitHub/GitHub.jl # base/markdown/GitHub/table.jl # base/markdown/IPython/IPython.jl # base/markdown/Julia/Julia.jl # base/markdown/Julia/interp.jl # base/markdown/Markdown.jl # base/markdown/parse/config.jl # base/markdown/parse/parse.jl # base/markdown/parse/util.jl # base/markdown/render/html.jl # base/markdown/render/latex.jl # base/markdown/render/plain.jl # base/markdown/render/rich.jl # base/markdown/render/terminal/formatting.jl # base/markdown/render/terminal/render.jl # base/math.jl # base/meta.jl # base/methodshow.jl # base/mmap.jl # base/mpfr.jl # base/multi.jl # base/multidimensional.jl # base/multimedia.jl # base/nullable.jl # base/number.jl # base/operators.jl # base/options.jl # base/ordering.jl # base/osutils.jl # base/path.jl # base/pcre.jl # base/pkg.jl # base/pkg/cache.jl # base/pkg/dir.jl # base/pkg/entry.jl # base/pkg/generate.jl # base/pkg/git.jl # base/pkg/github.jl # base/pkg/query.jl # base/pkg/read.jl # base/pkg/reqs.jl # base/pkg/resolve.jl # base/pkg/resolve/fieldvalue.jl # base/pkg/resolve/interface.jl # base/pkg/resolve/maxsum.jl # base/pkg/resolve/versionweight.jl # base/pkg/types.jl # base/pkg/write.jl # base/pointer.jl # base/poll.jl # base/precompile.jl # base/printf.jl # base/process.jl # base/profile.jl # base/promotion.jl # base/quadgk.jl # base/random.jl # base/range.jl # base/rational.jl # base/reduce.jl # base/reducedim.jl # base/reflection.jl # base/refpointer.jl # base/regex.jl # base/replutil.jl # base/rounding.jl # base/serialize.jl # base/set.jl # base/sharedarray.jl # base/show.jl # base/simdloop.jl # base/socket.jl # base/sort.jl # base/sparse.jl # base/sparse/abstractsparse.jl # base/sparse/cholmod.jl # base/sparse/cholmod_h.jl # base/sparse/linalg.jl # base/sparse/sparsematrix.jl # base/sparse/spqr.jl # base/sparse/umfpack.jl # base/sparse/umfpack_h.jl # base/special/bessel.jl # base/special/erf.jl # base/special/gamma.jl # base/special/log.jl # base/stat.jl # base/statistics.jl # base/stream.jl # base/string.jl # base/subarray.jl # base/subarray2.jl # base/sysimg.jl # base/sysinfo.jl # base/task.jl # base/test.jl # base/tuple.jl # base/utf16.jl # base/utf32.jl # base/utf8.jl # base/utf8proc.jl # base/util.jl # base/version.jl # base/version_git.sh # contrib/build_executable.jl # contrib/build_sysimg.jl # contrib/check-whitespace.sh # contrib/filterArgs.sh # contrib/fixup-libgfortran.sh # contrib/fixup-libstdc++.sh # contrib/install.sh # contrib/julia-config.jl # contrib/mac/app/run-install-name-tool-change.sh # contrib/mac/juliarc.jl # contrib/mac/mac-gtk.sh # contrib/relative_path.sh # contrib/stringreplace.c # contrib/windows/juliarc.jl # contrib/windows/msys_build.sh # contrib/windows/winrpm.sh # examples/bubblesort.jl # examples/clustermanager/0mq/ZMQCM.jl # examples/clustermanager/0mq/broker.jl # examples/clustermanager/0mq/head.jl # examples/clustermanager/0mq/worker.jl # examples/clustermanager/simple/UnixDomainCM.jl # examples/clustermanager/simple/head.jl # examples/clustermanager/simple/test_simple.jl # examples/embedding.c # examples/hpl.jl # examples/juliatypes.jl # examples/lru.jl # examples/lru_test.jl # examples/modint.jl # examples/ndgrid.jl # examples/queens.jl # examples/quine.jl # examples/staged.jl # examples/time.jl # examples/typetree.jl # examples/wordcount.jl # src/alloc.c # src/array.c # src/ast.c # src/builtin_proto.h # src/builtins.c # src/ccall.cpp # src/cgutils.cpp # src/codegen.cpp # src/debuginfo.cpp # src/dlload.c # src/dump.c # src/fenv_constants.h # src/file_constants.h # src/gc.c # src/gf.c # src/init.c # src/interpreter.c # src/intrinsics.cpp # src/jl_uv.c # src/jlapi.c # src/jltypes.c # src/julia.h # src/julia_internal.h # src/llvm-simdloop.cpp # src/llvm-version.h # src/module.c # src/options.h # src/profile.c # src/simplevector.c # src/support/arraylist.c # src/support/arraylist.h # src/support/bitvector.c # src/support/bitvector.h # src/support/dirpath.h # src/support/dtypes.h # src/support/hashing.c # src/support/hashing.h # src/support/htable.c # src/support/htable.h # src/support/int2str.c # src/support/ios.c # src/support/ios.h # src/support/libsupport.h # src/support/libsupportinit.c # src/support/operators.c # src/support/platform.h # src/support/ptrhash.c # src/support/ptrhash.h # src/support/strtod.h # src/support/timefuncs.c # src/support/timefuncs.h # src/support/utf8.h # src/support/utils.h # src/sys.c # src/table.c # src/task.c # src/toplevel.c # src/uv_constants.h # test/TestHelpers.jl # test/arrayops.jl # test/arrayperf.jl # test/backtrace.jl # test/base64.jl # test/bigint.jl # test/bitarray.jl # test/blas.jl # test/broadcast.jl # test/ccall.jl # test/ccalltest.c # test/char.jl # test/choosetests.jl # test/cmdlineargs.jl # test/combinatorics.jl # test/complex.jl # test/copy.jl # test/core.jl # test/dates.jl # test/dates/accessors.jl # test/dates/adjusters.jl # test/dates/arithmetic.jl # test/dates/conversions.jl # test/dates/io.jl # test/dates/periods.jl # test/dates/query.jl # test/dates/ranges.jl # test/dates/types.jl # test/dict.jl # test/docs.jl # test/dsp.jl # test/enums.jl # test/euler.jl # test/examples.jl # test/fastmath.jl # test/fft.jl # test/file.jl # test/float16.jl # test/floatapprox.jl # test/functional.jl # test/functors.jl # test/git.jl # test/gitutils.jl # test/goto.jl # test/grisu.jl # test/hashing.jl # test/i18n.jl # test/intfuncs.jl # test/iobuffer.jl # test/keywordargs.jl # test/libgit2.jl # test/linalg/arnoldi.jl # test/linalg/bidiag.jl # test/linalg/cholesky.jl # test/linalg/diagonal.jl # test/linalg/givens.jl # test/linalg/lapack.jl # test/linalg/lu.jl # test/linalg/pinv.jl # test/linalg/symmetric.jl # test/linalg/triangular.jl # test/linalg/tridiag.jl # test/linalg1.jl # test/linalg2.jl # test/linalg3.jl # test/linalg4.jl # test/lineedit.jl # test/llvmcall.jl # test/markdown.jl # test/math.jl # test/meta.jl # test/misc.jl # test/mod2pi.jl # test/mpfr.jl # test/netload/memtest.jl # test/netload/nettest.jl # test/nullable.jl # test/numbers.jl # test/operators.jl # test/parallel.jl # test/parser.jl # test/path.jl # test/perf/blas/level1.jl # test/perf/blas/level2.jl # test/perf/blas/level3.jl # test/perf/blas/perf.jl # test/perf/cat/perf.jl # test/perf/kernel/actor_centrality.jl # test/perf/kernel/bench_eu.jl # test/perf/kernel/getdivgrad.jl # test/perf/kernel/gk.jl # test/perf/kernel/go_benchmark.c # test/perf/kernel/go_benchmark.jl # test/perf/kernel/indexing.jl # test/perf/kernel/json.jl # test/perf/kernel/laplace.jl # test/perf/kernel/laplace/c_laplace.c # test/perf/kernel/laplace/c_laplace_parallel_update.c # test/perf/kernel/laplace/c_laplace_parallel_update_pointer.c # test/perf/kernel/laplace/cilk_laplace.c # test/perf/kernel/perf.jl # test/perf/kernel/raytracer.jl # test/perf/kernel/simplex.jl # test/perf/kernel/stockcorr.jl # test/perf/kernel/ziggurat.jl # test/perf/lapack/eig.jl # test/perf/lapack/factorizations.jl # test/perf/lapack/perf.jl # test/perf/micro/java/setup.sh # test/perf/micro/perf.c # test/perf/micro/perf.jl # test/perf/perfcomp.jl # test/perf/perfgeneric.jl # test/perf/perfutil.jl # test/perf/report.jl # test/perf/shootout/binary_trees.jl # test/perf/shootout/fannkuch.jl # test/perf/shootout/fasta.jl # test/perf/shootout/k_nucleotide.jl # test/perf/shootout/mandelbrot.jl # test/perf/shootout/meteor_contest.jl # test/perf/shootout/nbody.jl # test/perf/shootout/nbody_vec.jl # test/perf/shootout/perf.jl # test/perf/shootout/pidigits.jl # test/perf/shootout/regex_dna.jl # test/perf/shootout/revcomp.jl # test/perf/shootout/spectralnorm.jl # test/perf/simd/axpy.jl # test/perf/simd/inner.jl # test/perf/simd/perf.jl # test/perf/simd/seismic_fdtd.jl # test/perf/simd/sum_reduce.jl # test/perf/sort/perf.jl # test/perf/sparse/fem.jl # test/perf/sparse/getindex.jl # test/perf/sparse/perf.jl # test/perf/spell/perf.jl # test/pkg.jl # test/pollfd.jl # test/priorityqueue.jl # test/profile.jl # test/random.jl # test/ranges.jl # test/readdlm.jl # test/reduce.jl # test/reducedim.jl # test/reflection.jl # test/regex.jl # test/remote.jl # test/repl.jl # test/replcompletions.jl # test/replutil.jl # test/resolve.jl # test/rounding.jl # test/runtests.jl # test/serialize.jl # test/sets.jl # test/show.jl # test/simdloop.jl # test/socket.jl # test/sorting.jl # test/sparse.jl # test/sparsedir/cholmod.jl # test/sparsedir/sparse.jl # test/sparsedir/spqr.jl # test/sparsedir/umfpack.jl # test/spawn.jl # test/staged.jl # test/statistics.jl # test/strings.jl # test/subarray.jl # test/sysinfo.jl # test/test.jl # test/test_sourcepath.jl # test/testdefs.jl # test/tuple.jl # test/unicode.jl # test/version.jl --- base/primes.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/base/primes.jl b/base/primes.jl index 769da3d..7bd6a0a 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -1,3 +1,5 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + # Sieve of Atkin for generating primes: # http://en.wikipedia.org/wiki/Sieve_of_Atkin # Code very loosely based on this: From 0d2243c0d45fedb0178745750326d1994623ae1c Mon Sep 17 00:00:00 2001 From: aiorla Date: Thu, 7 May 2015 19:56:06 +0200 Subject: [PATCH 23/35] Speed up factor(n) + minor tweaks in other "primes.jl" functions --- base/primes.jl | 92 +++++++++++++++++++++++++++++++------------------- 1 file changed, 57 insertions(+), 35 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 7bd6a0a..77bc2ad 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -10,7 +10,7 @@ function primesmask(s::AbstractVector{Bool}) n = length(s) n < 2 && return s; s[2] = true n < 3 && return s; s[3] = true - r = floor(Int,sqrt(n)) + r = isqrt(n) for x = 1:r xx = x*x for y = 1:r @@ -32,16 +32,17 @@ primesmask(n::Integer) = n <= typemax(Int) ? primesmask(Int(n)) : primes(n::Union(Integer,AbstractVector{Bool})) = find(primesmask(n)) -# Miller-Rabin for primality testing: +const PRIMES = primes(2^16) + +# Small precomputed primes + Miller-Rabin for primality testing: # http://en.wikipedia.org/wiki/Miller–Rabin_primality_test # -function isprime(n::Integer) - n == 2 && return true - (n < 2) | iseven(n) && return false +function isprime{T<:Integer}(n::T) + (n < 3 || iseven(n)) && return n == 2 + n <= 2^16 && return length(searchsorted(PRIMES,n)) == 1 s = trailing_zeros(n-1) d = (n-1) >>> s for a in witnesses(n) - a < n || break x = powermod(a,d,n) x == 1 && continue t = s @@ -73,51 +74,72 @@ isprime(n::UInt128) = isprime(n::Int128) = n < 2 ? false : n <= typemax(Int64) ? isprime(Int64(n)) : isprime(BigInt(n)) -# TODO: faster factorization algorithms? - -const PRIMES = primes(10000) +# Trial division of small (< 2^16) precomputed primes + +# Pollard rho's algorithm with Richard P. Brent optimizations +# http://en.wikipedia.org/wiki/Trial_division +# http://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm +# http://maths-people.anu.edu.au/~brent/pub/pub051.html +# function factor{T<:Integer}(n::T) 0 < n || throw(ArgumentError("number to be factored must be ≥ 0, got $n")) h = Dict{T,Int}() n == 1 && return h - n <= 3 && (h[n] = 1; return h) - local s::T, p::T - s = isqrt(n) + isprime(n) && (h[n] = 1; return h) + local p::T for p in PRIMES - if p > s - h[n] = 1 - return h - end if n % p == 0 + h[p] = get(h,p,0)+1 + n = oftype(n, div(n,p)) while n % p == 0 h[p] = get(h,p,0)+1 n = oftype(n, div(n,p)) end n == 1 && return h - if isprime(n) - h[n] = 1 - return h - end - s = isqrt(n) + isprime(n) && (h[n] = 1; return h) end end - p = PRIMES[end]+2 - while p <= s - if n % p == 0 - while n % p == 0 - h[p] = get(h,p,0)+1 - n = oftype(n, div(n,p)) + pollardfactors!(n, h) +end + +function pollardfactors!{T<:Integer}(n::T, h::Dict{T,Int}) + while true + local c::T = rand(1:(n-1)), G::T = 1, r::T = 1, y::T = rand(0:(n-1)), m::T = 1900 + local ys::T, q::T = 1, x::T + while c == n - 2 + c = rand(1:(n-1)) + end + while G == 1 + x = y + for i in 1:r + y = oftype(y, widemul(y,y)%n) + y = oftype(y, (widen(y)+widen(c))%n) end - n == 1 && return h - if isprime(n) - h[n] = 1 - return h + local k::T = 0 + G = 1 + while k < r && G == 1 + for i in 1:(m>(r-k)?(r-k):m) + ys = y + y = oftype(y, widemul(y,y)%n) + y = oftype(y, (widen(y)+widen(c))%n) + q = oftype(y, widemul(q,x>y?x-y:y-x)%n) + end + G = gcd(q,n) + k = k + m end - s = isqrt(n) + r = 2 * r + end + G == n && (G = 1) + while G == 1 + ys = oftype(ys, widemul(ys,ys)%n) + ys = oftype(ys, (widen(ys)+widen(c))%n) + G = gcd(x>ys?x-ys:ys-x,n) + end + if G != n + isprime(G) ? h[G] = get(h,G,0) + 1 : pollardfactors!(G,h) + G2 = oftype(n,div(n,G)) + isprime(G2) ? h[G2] = get(h,G2,0) + 1 : pollardfactors!(G2,h) + return h end - p += 2 end - h[n] = 1 - return h end From 3ecf94b5022104b95d7a028bc76ae098f61b7d8a Mon Sep 17 00:00:00 2001 From: aiorla Date: Thu, 7 May 2015 22:38:18 +0200 Subject: [PATCH 24/35] Revert to previous header of primes(n) + 2 tests for big factorizations. # Conflicts: # test/numbers.jl --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index 77bc2ad..5839fb9 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -37,7 +37,7 @@ const PRIMES = primes(2^16) # Small precomputed primes + Miller-Rabin for primality testing: # http://en.wikipedia.org/wiki/Miller–Rabin_primality_test # -function isprime{T<:Integer}(n::T) +function isprime(n::Integer) (n < 3 || iseven(n)) && return n == 2 n <= 2^16 && return length(searchsorted(PRIMES,n)) == 1 s = trailing_zeros(n-1) From 8f684cfa480c8e29f00f39806e783e73e68eba4d Mon Sep 17 00:00:00 2001 From: aiorla Date: Tue, 26 May 2015 11:14:23 +0200 Subject: [PATCH 25/35] Redundant oftypes removes + Use of a more eficient searchsorted. --- base/primes.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 5839fb9..0094b36 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -39,7 +39,7 @@ const PRIMES = primes(2^16) # function isprime(n::Integer) (n < 3 || iseven(n)) && return n == 2 - n <= 2^16 && return length(searchsorted(PRIMES,n)) == 1 + n <= 2^16 && return PRIMES[searchsortedlast(PRIMES,n)] == n s = trailing_zeros(n-1) d = (n-1) >>> s for a in witnesses(n) @@ -90,10 +90,10 @@ function factor{T<:Integer}(n::T) for p in PRIMES if n % p == 0 h[p] = get(h,p,0)+1 - n = oftype(n, div(n,p)) + n = div(n,p) while n % p == 0 h[p] = get(h,p,0)+1 - n = oftype(n, div(n,p)) + n = div(n,p) end n == 1 && return h isprime(n) && (h[n] = 1; return h) @@ -112,17 +112,17 @@ function pollardfactors!{T<:Integer}(n::T, h::Dict{T,Int}) while G == 1 x = y for i in 1:r - y = oftype(y, widemul(y,y)%n) - y = oftype(y, (widen(y)+widen(c))%n) + y = widemul(y,y)%n + y = (widen(y)+widen(c))%n end local k::T = 0 G = 1 while k < r && G == 1 for i in 1:(m>(r-k)?(r-k):m) ys = y - y = oftype(y, widemul(y,y)%n) - y = oftype(y, (widen(y)+widen(c))%n) - q = oftype(y, widemul(q,x>y?x-y:y-x)%n) + y = widemul(y,y)%n + y = (widen(y)+widen(c))%n + q = widemul(q,x>y?x-y:y-x)%n end G = gcd(q,n) k = k + m @@ -131,13 +131,13 @@ function pollardfactors!{T<:Integer}(n::T, h::Dict{T,Int}) end G == n && (G = 1) while G == 1 - ys = oftype(ys, widemul(ys,ys)%n) - ys = oftype(ys, (widen(ys)+widen(c))%n) + ys = widemul(ys,ys)%n + ys = (widen(ys)+widen(c))%n G = gcd(x>ys?x-ys:ys-x,n) end if G != n isprime(G) ? h[G] = get(h,G,0) + 1 : pollardfactors!(G,h) - G2 = oftype(n,div(n,G)) + G2 = div(n,G) isprime(G2) ? h[G2] = get(h,G2,0) + 1 : pollardfactors!(G2,h) return h end From f28ef1189fd6fbf80f852aeb20979b770fe58765 Mon Sep 17 00:00:00 2001 From: aiorla Date: Thu, 28 May 2015 23:39:33 +0200 Subject: [PATCH 26/35] Fixing types conflict in pollardfactors! --- base/primes.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 0094b36..a3971b7 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -99,12 +99,12 @@ function factor{T<:Integer}(n::T) isprime(n) && (h[n] = 1; return h) end end - pollardfactors!(n, h) + widemul(n-1,n-1) <= typemax(n) ? pollardfactors!(n, h) : pollardfactors!(widen(n), h) end -function pollardfactors!{T<:Integer}(n::T, h::Dict{T,Int}) +function pollardfactors!{T<:Integer,K<:Integer}(n::T, h::Dict{K,Int}) while true - local c::T = rand(1:(n-1)), G::T = 1, r::T = 1, y::T = rand(0:(n-1)), m::T = 1900 + local c::T = rand(1:(n-1)), G::T = 1, r::K = 1, y::T = rand(0:(n-1)), m::K = 1900 local ys::T, q::T = 1, x::T while c == n - 2 c = rand(1:(n-1)) @@ -112,17 +112,17 @@ function pollardfactors!{T<:Integer}(n::T, h::Dict{T,Int}) while G == 1 x = y for i in 1:r - y = widemul(y,y)%n - y = (widen(y)+widen(c))%n + y = (y*y)%n + y = (y+c)%n end - local k::T = 0 + local k::K = 0 G = 1 while k < r && G == 1 for i in 1:(m>(r-k)?(r-k):m) ys = y - y = widemul(y,y)%n - y = (widen(y)+widen(c))%n - q = widemul(q,x>y?x-y:y-x)%n + y = (y*y)%n + y = (y+c)%n + q = (q*(x>y?x-y:y-x))%n end G = gcd(q,n) k = k + m @@ -131,8 +131,8 @@ function pollardfactors!{T<:Integer}(n::T, h::Dict{T,Int}) end G == n && (G = 1) while G == 1 - ys = widemul(ys,ys)%n - ys = (widen(ys)+widen(c))%n + ys = (ys*ys)%n + ys = (ys+c)%n G = gcd(x>ys?x-ys:ys-x,n) end if G != n From 8f48a460f63129eee70864019dccad61a72c8a15 Mon Sep 17 00:00:00 2001 From: aiorla Date: Fri, 29 May 2015 00:21:56 +0200 Subject: [PATCH 27/35] Skipping the typemax check for BigInts --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index a3971b7..4123afa 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -99,7 +99,7 @@ function factor{T<:Integer}(n::T) isprime(n) && (h[n] = 1; return h) end end - widemul(n-1,n-1) <= typemax(n) ? pollardfactors!(n, h) : pollardfactors!(widen(n), h) + T <: BigInt || widemul(n-1,n-1) <= typemax(n) ? pollardfactors!(n, h) : pollardfactors!(widen(n), h) end function pollardfactors!{T<:Integer,K<:Integer}(n::T, h::Dict{K,Int}) From 7756ed3cc746d171e9d2a4bdd67fc1568e029627 Mon Sep 17 00:00:00 2001 From: Jeff Bezanson Date: Sun, 14 Jun 2015 23:15:52 -0400 Subject: [PATCH 28/35] replace Union( ) with Union{ } everywhere # Conflicts: # base/LineEdit.jl # base/REPL.jl # base/REPLCompletions.jl # base/abstractarray.jl # base/array.jl # base/bitarray.jl # base/broadcast.jl # base/c.jl # base/client.jl # base/combinatorics.jl # base/complex.jl # base/constants.jl # base/coreimg.jl # base/datafmt.jl # base/dates/adjusters.jl # base/dates/periods.jl # base/dates/types.jl # base/deepcopy.jl # base/deprecated.jl # base/dict.jl # base/dsp.jl # base/expr.jl # base/fastmath.jl # base/fftw.jl # base/float.jl # base/fs.jl # base/gmp.jl # base/grisu.jl # base/inference.jl # base/int.jl # base/intfuncs.jl # base/libc.jl # base/libdl.jl # base/linalg.jl # base/linalg/bidiag.jl # base/linalg/bitarray.jl # base/linalg/blas.jl # base/linalg/cholesky.jl # base/linalg/dense.jl # base/linalg/eigen.jl # base/linalg/lu.jl # base/linalg/matmul.jl # base/linalg/qr.jl # base/linalg/schur.jl # base/linalg/special.jl # base/linalg/svd.jl # base/linalg/symmetric.jl # base/linalg/triangular.jl # base/linalg/uniformscaling.jl # base/loading.jl # base/markdown/render/html.jl # base/mmap.jl # base/mpfr.jl # base/multi.jl # base/multidimensional.jl # base/multimedia.jl # base/nofloat_hashing.jl # base/nullable.jl # base/operators.jl # base/ordering.jl # base/pkg.jl # base/pkg/entry.jl # base/pkg/generate.jl # base/pkg/query.jl # base/pkg/reqs.jl # base/pkg/resolve/versionweight.jl # base/pointer.jl # base/precompile.jl # base/printf.jl # base/process.jl # base/profile.jl # base/random.jl # base/range.jl # base/reduce.jl # base/reducedim.jl # base/regex.jl # base/replutil.jl # base/rounding.jl # base/sharedarray.jl # base/show.jl # base/socket.jl # base/sort.jl # base/sparse/cholmod.jl # base/sparse/cholmod_h.jl # base/sparse/csparse.jl # base/sparse/sparsematrix.jl # base/sparse/umfpack.jl # base/special/gamma.jl # base/stat.jl # base/statistics.jl # base/stream.jl # base/string.jl # base/subarray.jl # base/utf16.jl # base/utf32.jl # base/utf8.jl # base/utftypes.jl # base/version.jl # test/arrayops.jl # test/bitarray.jl # test/core.jl # test/functional.jl # test/linalg2.jl # test/numbers.jl # test/perf/array/indexing.jl # test/reducedim.jl # test/statistics.jl # test/tuple.jl --- base/primes.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 4123afa..62a2a02 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -30,7 +30,7 @@ primesmask(n::Int) = primesmask(falses(n)) primesmask(n::Integer) = n <= typemax(Int) ? primesmask(Int(n)) : throw(ArgumentError("requested number of primes must be ≤ $(typemax(Int)), got $n")) -primes(n::Union(Integer,AbstractVector{Bool})) = find(primesmask(n)) +primes(n::Union{Integer,AbstractVector{Bool}}) = find(primesmask(n)) const PRIMES = primes(2^16) @@ -60,9 +60,9 @@ end # http://primes.utm.edu/prove/merged.html # http://miller-rabin.appspot.com # -witnesses(n::Union(UInt8,Int8,UInt16,Int16)) = (2,3) -witnesses(n::Union(UInt32,Int32)) = n < 1373653 ? (2,3) : (2,7,61) -witnesses(n::Union(UInt64,Int64)) = +witnesses(n::Union{UInt8,Int8,UInt16,Int16}) = (2,3) +witnesses(n::Union{UInt32,Int32}) = n < 1373653 ? (2,3) : (2,7,61) +witnesses(n::Union{UInt64,Int64}) = n < 1373653 ? (2,3) : n < 4759123141 ? (2,7,61) : n < 2152302898747 ? (2,3,5,7,11) : From efb82334dd77ced2295e9df01e4952a9d58eb725 Mon Sep 17 00:00:00 2001 From: Yichao Yu Date: Mon, 27 Jul 2015 01:54:03 -0400 Subject: [PATCH 29/35] Replace http links with https links to suppress warning during doc linkcheck # Conflicts: # DISTRIBUTING.md # README.md # README.windows.md # doc/devdocs/backtraces.rst # doc/devdocs/object.rst # doc/devdocs/stdio.rst # doc/manual/arrays.rst # doc/manual/calling-c-and-fortran-code.rst # doc/manual/constructors.rst # doc/manual/dates.rst # doc/manual/documentation.rst # doc/manual/functions.rst # doc/manual/integers-and-floating-point-numbers.rst # doc/manual/introduction.rst # doc/manual/linear-algebra.rst # doc/manual/mathematical-operations.rst # doc/manual/metaprogramming.rst # doc/manual/methods.rst # doc/manual/profile.rst # doc/manual/strings.rst # doc/manual/style-guide.rst # doc/manual/types.rst # doc/manual/variables-and-scoping.rst # test/dates/accessors.jl # test/dates/conversions.jl --- base/primes.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 62a2a02..8368ee6 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -1,7 +1,7 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license # Sieve of Atkin for generating primes: -# http://en.wikipedia.org/wiki/Sieve_of_Atkin +# https://en.wikipedia.org/wiki/Sieve_of_Atkin # Code very loosely based on this: # http://thomasinterestingblog.wordpress.com/2011/11/30/generating-primes-with-the-sieve-of-atkin-in-c/ # http://dl.dropboxusercontent.com/u/29023244/atkin.cpp @@ -35,7 +35,7 @@ primes(n::Union{Integer,AbstractVector{Bool}}) = find(primesmask(n)) const PRIMES = primes(2^16) # Small precomputed primes + Miller-Rabin for primality testing: -# http://en.wikipedia.org/wiki/Miller–Rabin_primality_test +# https://en.wikipedia.org/wiki/Miller–Rabin_primality_test # function isprime(n::Integer) (n < 3 || iseven(n)) && return n == 2 @@ -77,8 +77,8 @@ isprime(n::Int128) = n < 2 ? false : # Trial division of small (< 2^16) precomputed primes + # Pollard rho's algorithm with Richard P. Brent optimizations -# http://en.wikipedia.org/wiki/Trial_division -# http://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm +# https://en.wikipedia.org/wiki/Trial_division +# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm # http://maths-people.anu.edu.au/~brent/pub/pub051.html # function factor{T<:Integer}(n::T) From c85ba8f45e44b3af45548d223c239425adb99487 Mon Sep 17 00:00:00 2001 From: pabloferz Date: Mon, 6 Jul 2015 00:04:41 -0500 Subject: [PATCH 30/35] Use a wheel sieve to improve performance and export primesmask. # Conflicts: # NEWS.md # base/docs/helpdb.jl # base/exports.jl # doc/stdlib/numbers.rst # test/numbers.jl --- base/primes.jl | 131 ++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 107 insertions(+), 24 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index 8368ee6..edd55de 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -1,36 +1,119 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license -# Sieve of Atkin for generating primes: -# https://en.wikipedia.org/wiki/Sieve_of_Atkin -# Code very loosely based on this: -# http://thomasinterestingblog.wordpress.com/2011/11/30/generating-primes-with-the-sieve-of-atkin-in-c/ -# http://dl.dropboxusercontent.com/u/29023244/atkin.cpp -# -function primesmask(s::AbstractVector{Bool}) - n = length(s) - n < 2 && return s; s[2] = true - n < 3 && return s; s[3] = true - r = isqrt(n) - for x = 1:r - xx = x*x - for y = 1:r - yy = y*y - i, j, k = 4xx+yy, 3xx+yy, 3xx-yy - i <= n && (s[i] $= (i%12==1)|(i%12==5)) - j <= n && (s[j] $= (j%12==7)) - 1 <= k <= n && (s[k] $= (x>y)&(k%12==11)) +# Primes generating functions +# https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes +# https://en.wikipedia.org/wiki/Wheel_factorization +# http://primesieve.org +# Jonathan Sorenson, "An analysis of two prime number sieves", Computer Science Technical Report Vol. 1028, 1991 +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, 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+1] ) +@inline wheel_prime(n) = ( (d, r) = ((n - 1) >>> 3, (n - 1) & 7); return 30d + wheel_primes[r+1] ) + +function _primesmask(limit::Int) + limit < 7 && throw(ArgumentError("limit must be at least 7, got $limit")) + n = wheel_index(limit) + m = wheel_prime(n) + sieve = ones(Bool, n) + @inbounds for i = 1:wheel_index(isqrt(limit)) + if sieve[i]; p = wheel_prime(i) + q = p * p + j = (i - 1) & 7 + 1 + while q ≤ m + sieve[wheel_index(q)] = false + q = q + wheel[j] * p + j = j & 7 + 1 + end end end - for i = 5:r - s[i] && (s[i*i:i*i:n] = false) + return sieve +end + +function _primesmask(lo::Int, hi::Int) + 7 ≤ lo ≤ hi || throw(ArgumentError("the condition 7 ≤ lo ≤ hi must be met")) + lo == 7 && return _primesmask(hi) + wlo, whi = wheel_index(lo - 1), wheel_index(hi) + m = wheel_prime(whi) + sieve = ones(Bool, whi - wlo) + small_primes = primes(isqrt(hi)) + @inbounds for i in 4:length(small_primes) + p = small_primes[i] + j = wheel_index(2 * div(lo - p - 1, 2p) + 1) + q = p * wheel_prime(j + 1); j = j & 7 + 1 + while q ≤ m + sieve[wheel_index(q)-wlo] = false + q = q + wheel[j] * p + j = j & 7 + 1 + end + end + return sieve +end + +# Sieve of the primes up to limit represented as an array of booleans +function primesmask(limit::Int) + limit < 1 && throw(ArgumentError("limit must be at least 1, got $limit")) + sieve = falses(limit) + limit < 2 && return sieve; sieve[2] = true + limit < 3 && return sieve; sieve[3] = true + limit < 5 && return sieve; sieve[5] = true + limit < 7 && return sieve + wheel_sieve = _primesmask(limit) + @inbounds for i in eachindex(wheel_sieve) + Base.unsafe_setindex!(sieve, wheel_sieve[i], wheel_prime(i)) end - return s + return sieve end -primesmask(n::Int) = primesmask(falses(n)) primesmask(n::Integer) = n <= typemax(Int) ? primesmask(Int(n)) : throw(ArgumentError("requested number of primes must be ≤ $(typemax(Int)), got $n")) -primes(n::Union{Integer,AbstractVector{Bool}}) = find(primesmask(n)) +function primesmask(lo::Int, hi::Int) + 0 < lo ≤ hi || throw(ArgumentError("the condition 0 < lo ≤ hi must be met")) + sieve = falses(hi - lo + 1) + lo ≤ 2 ≤ hi && (sieve[3-lo] = true) + lo ≤ 3 ≤ hi && (sieve[4-lo] = true) + lo ≤ 5 ≤ hi && (sieve[6-lo] = true) + hi < 7 && return sieve + wheel_sieve = _primesmask(max(7, lo), hi) + lsi = lo - 1 + lwi = wheel_index(lsi) + @inbounds for i in eachindex(wheel_sieve) + Base.unsafe_setindex!(sieve, wheel_sieve[i], wheel_prime(i + lwi) - lsi) + end + return sieve +end +primesmask{T<:Integer}(lo::T, hi::T) = lo <= hi <= typemax(Int) ? primesmask(Int(lo), Int(hi)) : + throw(ArgumentError("both endpoints of the interval to sieve must be ≤ $(typemax(Int)), got $lo and $hi")) + +function primes(n::Int) + list = Int[] + n < 2 && return list; push!(list, 2) + n < 3 && return list; push!(list, 3) + n < 5 && return list; push!(list, 5) + n < 7 && return list + sizehint!(list, floor(Int, n / log(n))) + sieve = _primesmask(n) + @inbounds for i in eachindex(sieve) + sieve[i] && push!(list, wheel_prime(i)) + end + return list +end + +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 + sieve = _primesmask(max(7, lo), hi) + lwi = wheel_index(lo - 1) + @inbounds for i in eachindex(sieve) + sieve[i] && push!(list, wheel_prime(i + lwi)) + end + return list +end const PRIMES = primes(2^16) From 0725830bbd3a8e9b562b60fb6bdb7d7d795d70a7 Mon Sep 17 00:00:00 2001 From: pabloferz Date: Tue, 12 Jan 2016 13:22:59 +0100 Subject: [PATCH 31/35] Simplify primes and primesmask [fixes #14652] --- base/primes.jl | 60 ++++++++++++++++---------------------------------- 1 file changed, 19 insertions(+), 41 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index edd55de..a5d85ef 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -7,9 +7,9 @@ # Jonathan Sorenson, "An analysis of two prime number sieves", Computer Science Technical Report Vol. 1028, 1991 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, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7 ] +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+1] ) +@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] ) function _primesmask(limit::Int) @@ -37,37 +37,23 @@ function _primesmask(lo::Int, hi::Int) wlo, whi = wheel_index(lo - 1), wheel_index(hi) m = wheel_prime(whi) sieve = ones(Bool, whi - wlo) - small_primes = primes(isqrt(hi)) - @inbounds for i in 4:length(small_primes) - p = small_primes[i] - j = wheel_index(2 * div(lo - p - 1, 2p) + 1) - q = p * wheel_prime(j + 1); j = j & 7 + 1 - while q ≤ m - sieve[wheel_index(q)-wlo] = false - q = q + wheel[j] * p - j = j & 7 + 1 + hi < 49 && return sieve + small_sieve = _primesmask(isqrt(hi)) + @inbounds for i in eachindex(small_sieve) + if small_sieve[i]; p = wheel_prime(i) + j = wheel_index(2 * div(lo - p - 1, 2p) + 1) + q = p * wheel_prime(j + 1); j = j & 7 + 1 + while q ≤ m + sieve[wheel_index(q)-wlo] = false + q = q + wheel[j] * p + j = j & 7 + 1 + end end end return sieve end -# Sieve of the primes up to limit represented as an array of booleans -function primesmask(limit::Int) - limit < 1 && throw(ArgumentError("limit must be at least 1, got $limit")) - sieve = falses(limit) - limit < 2 && return sieve; sieve[2] = true - limit < 3 && return sieve; sieve[3] = true - limit < 5 && return sieve; sieve[5] = true - limit < 7 && return sieve - wheel_sieve = _primesmask(limit) - @inbounds for i in eachindex(wheel_sieve) - Base.unsafe_setindex!(sieve, wheel_sieve[i], wheel_prime(i)) - end - return sieve -end -primesmask(n::Integer) = n <= typemax(Int) ? primesmask(Int(n)) : - throw(ArgumentError("requested number of primes must be ≤ $(typemax(Int)), got $n")) - +# Sieve of the primes from lo up to hi represented as an array of booleans function primesmask(lo::Int, hi::Int) 0 < lo ≤ hi || throw(ArgumentError("the condition 0 < lo ≤ hi must be met")) sieve = falses(hi - lo + 1) @@ -86,19 +72,9 @@ end primesmask{T<:Integer}(lo::T, hi::T) = lo <= hi <= typemax(Int) ? primesmask(Int(lo), Int(hi)) : throw(ArgumentError("both endpoints of the interval to sieve must be ≤ $(typemax(Int)), got $lo and $hi")) -function primes(n::Int) - list = Int[] - n < 2 && return list; push!(list, 2) - n < 3 && return list; push!(list, 3) - n < 5 && return list; push!(list, 5) - n < 7 && return list - sizehint!(list, floor(Int, n / log(n))) - sieve = _primesmask(n) - @inbounds for i in eachindex(sieve) - sieve[i] && push!(list, wheel_prime(i)) - end - return list -end +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")) function primes(lo::Int, hi::Int) lo ≤ hi || throw(ArgumentError("the condition lo ≤ hi must be met")) @@ -107,6 +83,7 @@ function primes(lo::Int, hi::Int) lo ≤ 3 ≤ hi && push!(list, 3) lo ≤ 5 ≤ hi && push!(list, 5) hi < 7 && return list + sizehint!(list, floor(Int, hi / log(hi))) sieve = _primesmask(max(7, lo), hi) lwi = wheel_index(lo - 1) @inbounds for i in eachindex(sieve) @@ -114,6 +91,7 @@ function primes(lo::Int, hi::Int) end return list end +primes(n::Int) = primes(1, n) const PRIMES = primes(2^16) From 8324752a4279d300afdcddc1d7b743cf9c961cd6 Mon Sep 17 00:00:00 2001 From: Matt Bauman Date: Sat, 6 Feb 2016 13:56:28 -0500 Subject: [PATCH 32/35] Replace unsafe_(get|set)index with inbounds macro This patch does two things: * It all but completely removes `unsafe_getindex` and `unsafe_setindex!` from the standard library. This merges their method definitions into the safe flavors with a `(at)boundscheck checkbounds(...)` clause, and it replaces call-sites with `(at)inbounds A[I...]`. A few uses remain because the inbounds macro doesn't return its value, making the construct a little awkward in some cases. * It changes *where* bounds are checked in SubArrays. Previously, SubArrays would defer bounds checks to their parent arrays; they wouldn't check the indices directly, but rather they'd transform the indices into the proper references into the parent array, and then index into the parent array safely. E.g., previously: julia> A = -.5:.1:.5 S = sub(A, 7:11) S[-4], S[0], S[4] (-0.4,0.0,0.4) julia> S[-6] ERROR: BoundsError: attempt to access -0.5:0.1:0.5 at index [0] This behavior is neither tested nor depended upon in the Base library, but I believe some packages have historically depended upon this. With this patch, the behavior is much more sane; `S[0]` now immediately throws a bounds error from the SubArray itself. I believe this is a strong requirement for making views more prominent (regardless of the syntax). # Conflicts: # base/abstractarray.jl # base/array.jl # base/bitarray.jl # base/linalg.jl # base/linalg/diagonal.jl # base/linalg/symmetric.jl # base/multidimensional.jl # base/number.jl # base/range.jl # base/simdloop.jl # base/subarray.jl --- base/primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/primes.jl b/base/primes.jl index a5d85ef..bcfc3fa 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -65,7 +65,7 @@ function primesmask(lo::Int, hi::Int) lsi = lo - 1 lwi = wheel_index(lsi) @inbounds for i in eachindex(wheel_sieve) - Base.unsafe_setindex!(sieve, wheel_sieve[i], wheel_prime(i + lwi) - lsi) + sieve[wheel_prime(i + lwi) - lsi] = wheel_sieve[i] end return sieve end From dcdf5863bae71d55d1b230365348978011a643e1 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Thu, 10 Mar 2016 05:24:31 -0600 Subject: [PATCH 33/35] eachindex: use separate indexes for each array This also replaces several instances of `for i in eachindex(A); a = A[i]` with `for a in A`. The latter is presumably easier for automatic bounds-checking removal. # Conflicts: # base/abstractarray.jl # base/arraymath.jl # base/combinatorics.jl # base/complex.jl # base/dates/periods.jl # base/floatfuncs.jl # base/io.jl # base/linalg/generic.jl # base/linalg/lapack.jl # base/math.jl # base/multidimensional.jl # base/operators.jl # base/sparse/linalg.jl --- base/primes.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/base/primes.jl b/base/primes.jl index bcfc3fa..52d57d6 100644 --- a/base/primes.jl +++ b/base/primes.jl @@ -39,7 +39,7 @@ function _primesmask(lo::Int, hi::Int) sieve = ones(Bool, whi - wlo) hi < 49 && return sieve small_sieve = _primesmask(isqrt(hi)) - @inbounds for i in eachindex(small_sieve) + @inbounds for i = 1:length(small_sieve) # don't use eachindex here if small_sieve[i]; p = wheel_prime(i) j = wheel_index(2 * div(lo - p - 1, 2p) + 1) q = p * wheel_prime(j + 1); j = j & 7 + 1 @@ -64,7 +64,7 @@ function primesmask(lo::Int, hi::Int) wheel_sieve = _primesmask(max(7, lo), hi) lsi = lo - 1 lwi = wheel_index(lsi) - @inbounds for i in eachindex(wheel_sieve) + @inbounds for i = 1:length(wheel_sieve) # don't use eachindex here sieve[wheel_prime(i + lwi) - lsi] = wheel_sieve[i] end return sieve @@ -86,7 +86,7 @@ function primes(lo::Int, hi::Int) sizehint!(list, floor(Int, hi / log(hi))) sieve = _primesmask(max(7, lo), hi) lwi = wheel_index(lo - 1) - @inbounds for i in eachindex(sieve) + @inbounds for i = 1:length(sieve) # don't use eachindex here sieve[i] && push!(list, wheel_prime(i + lwi)) end return list From 75582d9447dac5924bf1e3dbf92f98627582252d Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sun, 29 May 2016 13:57:19 +0530 Subject: [PATCH 34/35] rename base/primes.jl to src/Primes.jl --- base/primes.jl => src/Primes.jl | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename base/primes.jl => src/Primes.jl (100%) diff --git a/base/primes.jl b/src/Primes.jl similarity index 100% rename from base/primes.jl rename to src/Primes.jl From a4de93f0880ea4085736bc46dfdc3111bb169cc9 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Sun, 29 May 2016 13:58:00 +0530 Subject: [PATCH 35/35] update src/Primes.jl to master's version --- src/Primes.jl | 84 ++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 76 insertions(+), 8 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index 52d57d6..1241201 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -1,4 +1,13 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license +__precompile__() +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,6 +96,11 @@ 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[] @@ -83,7 +108,7 @@ function primes(lo::Int, hi::Int) 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