Merge pull request JuliaLang#11189 from aiorla/master

Speed up factor(n) + minor tweaks in other "primes.jl" functions

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:https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
 #
 function isprime(n::Integer)
- n == 2 && return true
- (n < 2) | iseven(n) && return false
+ (n < 3 || iseven(n)) && return n == 2
+ n <= 2^16 && return PRIMES[searchsortedlast(PRIMES,n)] == n
  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:https://en.wikipedia.org/wiki/Trial_division
+# http:https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
+# http:https://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 = 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
- 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 = widemul(y,y)%n
+ 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 = 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
  end
- s = isqrt(n)
+ r = 2 * r
+ end
+ G == n && (G = 1)
+ while G == 1
+ 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 = 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

test/numbers.jl

@@ -2092,6 +2092,10 @@ end
 @test 2.0f0^(1//3) == 2.0f0^(1.0f0/3)
 @test 2^(1//3) == 2^(1/3)
 
+# factorization of factors > 2^16
+@test factor((big(2)^31-1)^2) == Dict(big(2^31-1) => 2)
+@test factor((big(2)^31-1)*(big(2)^17-1)) == Dict(big(2^31-1) => 1, big(2^17-1) => 1)
+
 # large shift amounts
 @test Int32(-1)>>31 == -1
 @test Int32(-1)>>32 == -1