Skip to content

Commit 54a6be0

Browse files
author
oscarddssmith
committed
only store min factor for odd numbers
1 parent f83d703 commit 54a6be0

File tree

1 file changed

+31
-25
lines changed

1 file changed

+31
-25
lines changed

src/Primes.jl

Lines changed: 31 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -130,27 +130,26 @@ function primes(lo::Int, hi::Int)
130130
end
131131
primes(n::Int) = primes(1, n)
132132

133-
# internal function for generating the minimum factor of a relatively small number
134-
function _min_factor(n)
135-
n == 1 && return 1
136-
for i in 2:isqrt(n)
137-
n%i == 0 && return i
133+
function _generate_min_factors(limit)
134+
function min_factor(n)
135+
n < 4 && return n
136+
for i in 3:2:isqrt(n)
137+
n%i == 0 && return i
138+
end
139+
return n
138140
end
139-
return n
140-
end
141-
142-
function _min_factors(limit)
143141
res = Int[]
144-
for i in 1:limit
145-
m = _min_factor(i)
142+
for i in 3:2:limit
143+
m = min_factor(i)
146144
push!(res, m==i ? 1 : m)
147145
end
148146
return res
149147
end
150148

151-
# `MIN_FACTOR[i]` is `1` if `i` is `1` or prime
152-
# otherwise, `MIN_FACTOR[i]` is the first number which factors `i`
153-
const MIN_FACTOR = UInt8.(_min_factors(2^16))
149+
const N_SMALL_FACTORS = 2^16
150+
const _MIN_FACTOR = UInt8.(_generate_min_factors(N_SMALL_FACTORS))
151+
# _min_factor(n) = 1 if n is prime, otherwise the minimum factor of n
152+
_min_factor(n) = _MIN_FACTOR[n>>1]
154153

155154
"""
156155
isprime(n::Integer) -> Bool
@@ -166,11 +165,13 @@ function isprime(n::Integer)
166165
# Small precomputed primes + Miller-Rabin for primality testing:
167166
# https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
168167
# https://github.com/JuliaLang/julia/issues/11594
169-
if n < length(MIN_FACTOR)
170-
return n > 1 && MIN_FACTOR[n] == 1
168+
n < 2 && return false
169+
trailing_zeros(n) > 1 && return n==2
170+
if n < N_SMALL_FACTORS
171+
return _min_factor(n) == 1
171172
end
172-
for m in (2, 3, 5, 7, 11, 13, 17, 19, 23)
173-
n % m == 0 && return n == m
173+
for m in (3, 5, 7, 11, 13, 17, 19, 23)
174+
n % m == 0 && return false
174175
end
175176
s = trailing_zeros(n - 1)
176177
d = (n - 1) >>> s
@@ -270,22 +271,27 @@ function factor!(n::T, h::AbstractDict{K,Int}) where {T<:Integer,K<:Integer}
270271
elseif n == 0
271272
h[n] = 0
272273
return h
273-
elseif n <= length(MIN_FACTOR)
274+
elseif trailing_zeros(n) > 1
275+
tz = trailing_zeros(n)
276+
increment!(h, tz, 2)
277+
n >>= tz
278+
end
279+
if n <= N_SMALL_FACTORS
274280
while true
275281
n == 1 && return h
276-
if MIN_FACTOR[n]==1
282+
if _min_factor(n)==1
277283
return increment!(h, 1, n)
278284
else
279-
increment!(h, 1, MIN_FACTOR[n])
280-
n = div(n, MIN_FACTOR[n])
285+
increment!(h, 1, _min_factor(n))
286+
n = div(n, _min_factor(n))
281287
end
282288
end
283289
elseif isprime(n)
284290
return increment!(h, 1, n)
285291
end
286292
local p::T
287-
for p in 2:length(MIN_FACTOR)
288-
MIN_FACTOR[p] == 1 || continue
293+
for p in 3:2:N_SMALL_FACTORS
294+
_min_factor(p) == 1 || continue
289295
num_p = 0
290296
while true
291297
q, r = divrem(n, T(p)) # T(p) so julia <1.9 uses fast divrem for `BigInt`
@@ -297,7 +303,7 @@ function factor!(n::T, h::AbstractDict{K,Int}) where {T<:Integer,K<:Integer}
297303
if num_p > 0
298304
increment!(h, num_p, p)
299305
# if n is small, then recursing will hit the fast path.
300-
n < length(MIN_FACTOR) && return factor!(n, h)
306+
n < N_SMALL_FACTORS && return factor!(n, h)
301307
end
302308
p*p > n && break
303309
end

0 commit comments

Comments
 (0)