Find The Nearest Prime
August 25, 2017
The trick to this exercise is to realize that it’s a problem about binary search rather than a problem about prime numbers. Enumerate the first hundred primes in an array, then perform binary search. Here’s a simple Sieve of Eratosthenes to compute the primes less than n:
(define (primes n)
(let ((sieve (make-vector n #t)) (ps (list)))
(do ((p 2 (+ p 1))) ((= n p))
(when (vector-ref sieve p)
(set! ps (cons p ps))
(do ((i (* p p) (+ i p))) ((<= n i))
(vector-set! sieve i #f))))
(list->vector (reverse ps))))
For our application, we need the primes to 541, which is the hundredth prime:
(define ps (primes 542))
Normal binary search fails if the desired item isn’t in the input. We need a version of binary search that finds the nearest item if the requested item isn’t in the array:
(define (bsearch x xs)
(let loop ((lo 0) (hi (- (vector-length xs) 1)))
(let ((mid (+ lo (quotient (- hi lo) 2))))
(cond ((< hi lo)
(if (< (- (vector-ref xs hi) x)
(- x (vector-ref xs lo)))
(vector-ref xs lo)
(vector-ref xs hi)))
((< x (vector-ref xs mid)) (loop lo (- mid 1)))
((< (vector-ref xs mid) x) (loop (+ mid 1) hi))
(else (vector-ref xs mid))))))
Now we binary search the array of primes for the desired value:
(define (find-prime x) (bsearch x ps))
We take the lower number in case of a tie. Here are some examples:
> (find-prime 123.7) 127 > (find-prime 4) 3 > (find-prime 1000) error: out of bounds
You can run the program at https://ideone.com/xPH7v6.
In Python. This function returns the last prime, if x > last prime. bisect_left is from the bisect module.
def closest_prime(x, prs): index = bisect_left(prs, x) return min(prs[max(index-1, 0):index+1], key=lambda p: abs(p-x))I’m glad to see you providing a pragmatic solution. To be even more pragmatic, though, I’d just Google up a list of primes and copy and paste the first hundred of them.
Implementation in Julia. If it doesn’t look right for some reason, let me know and I’ll post a screenshot instead.
global P = primes(1000)[1:100] # first 100 primes
function NearestPrime(x::Int64)
if x > P[end]
return P[end]
else
d = abs(P – x)
return P[indmin(d)]
end
end
Quite a nifty little problem. It would be more of a challenge if the scope of it were larger though. Thanks
@Zack: I’m not familiar with Julia. Is `primes` a built-in function? Or part of the standard library? Does Julia provide integer factorization?
An alternative solution using a lazy prime generator. Generate primes until the distance between the pime and target increases.
def primegen(): yield from (2, 3, 5, 7, 11, 13) ps = primegen() p = ps.next() and ps.next() q, sieve, n = p**2, {}, 13 while True: if n not in sieve: if n < q: yield n else: next, step = q + 2*p, 2*p while next in sieve: next += step sieve[next] = step p = ps.next() q = p**2 else: step = sieve.pop(n) next = n + step while next in sieve: next += step sieve[next] = step n += 2 def closest_prime(x, limit=542): prs = lazyprimegen() last_prime = next(prs) dist = abs(x - last_prime) for prime in prs: d = abs(x - prime) if d >= dist or prime > limit: return last_prime dist, last_prime = d, prime