The Prime Ant
October 27, 2017
This is a simple matter of following the rules of the game:
(define (prime-ant n . verbose?)
(let ((a (make-vector (+ n 2))))
(do ((i 0 (+ i 1))) ((= (+ n 2) i))
(vector-set! a i (+ i 2)))
(let loop ((n n) (p 0))
(when (and (pair? verbose?) (equal? (car verbose?) #t))
(display a) (display " ") (display p) (newline))
(if (zero? n) p
(if (prime? (vector-ref a p))
(loop (- n 1) (+ p 1))
(let ((q (lpf (vector-ref a p))))
(vector-set! a p (/ (vector-ref a p) q))
(vector-set! a (- p 1) (+ (vector-ref a (- p 1)) q))
(loop (- n 1) (- p 1))))))))
We make use of two auxiliary functions: lpf returns the least prime factor of n, and prime? returns #t if and only if the least prime factor of n is n:
(define (lpf n) ; least prime factor
(let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((f 2) (w 0))
(if (< n (* f f)) n
(if (zero? (modulo n f)) f
(loop (+ f (vector-ref wheel w))
(if (= w 10) 3 (+ w 1))))))))
(define (prime? n) (= (lpf n) n))
Here are two examples:
> (prime-ant 20 #t) #(2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 0 #(2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 1 #(2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 2 #(2 5 2 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 1 #(2 5 2 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 2 #(2 5 2 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 3 #(2 5 2 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 4 #(2 5 2 7 3 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 3 #(2 5 2 7 3 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 4 #(2 5 2 7 3 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 5 #(2 5 2 7 3 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 6 #(2 5 2 7 3 9 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 5 #(2 5 2 7 6 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 4 #(2 5 2 9 3 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 3 #(2 5 5 3 3 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 2 #(2 5 5 3 3 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 3 #(2 5 5 3 3 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 4 #(2 5 5 3 3 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 5 #(2 5 5 3 3 3 4 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 6 #(2 5 5 3 3 5 2 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 5 #(2 5 5 3 3 5 2 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23) 6 6 > (prime-ant 47) 9
If it works forever, the prime ant will convert the entire array A to contain prime numbers. You can run the program at https://ideone.com/dok8t1.
Programming Praxis 
In Python:
import math def is_prime(a): for i in range(2, int(math.sqrt(a))+1): if a%i == 0: return i return True def prime_ant(turns): p = 0 #Starting index A = [] for i in range(turns): A.append(i+2) for n in range(turns): q = is_prime(A[p]) if q == True: p+=1 else: A[p] = A[p] // q A[p-1] = A[p-1]+q p-=1 return p print(prime_ant(47))Output: 9
Update: Cleaned up the code a bit, changed the function name and added some comments for clarity.
import math def divisor(a): # Takes an integer and returns the smallest divisor, or 0 if number is prime for i in range(2, int(math.sqrt(a))+1): if a%i == 0: return i return 0 def prime_ant(turns): p = 0 #Starting index A = list(range(2,turns)) #Generate initial array A for n in range(turns): q = divisor(A[p]) if q: A[p] = A[p] // q A[p-1] = A[p-1]+q p-=1 else: p+=1 return p print(prime_ant(47))Output: 9
In Python. The array for the ant grows when needed.
def smallest_factor(n): for f in accumulate(chain((2, 1, 2, 2), cycle((4, 2, 4, 2, 4, 6, 2, 6)))): if not n % f: return f if f * f > n: return n def prime_ant(): p, A, lenA = 0, [2], 1 while 1: yield p q = smallest_factor(A[p]) if q == A[p]: if lenA == p+1: # need to grow array? A.append(p+2) lenA += 1 p += 1 else: A[p] //= q A[p-1] += q p -= 1def is_prime(n, i=3): if int(n)&1 is 0 and n>2: return False; while i<n: if n%i == 0: return False i += 2 return True; def div(n, i=2): # return smallest divisor of n, assumes a non-prime number while n%i: i+=1 return i def prime_ant(n, p=0, i=0): A = {} # use an empty dict while i < n: if p not in A: A[p] = 2+p # set cell item if not set if is_prime(A[p]): p += 1 else: q = div(A[p]) A[p] = A[p]/q if p-1 not in A: A[p-1] = 2+(p-1) A[p-1] += q p -= 1 i += 1 return pUpdate:
def div(n, i=2): if n is 2: return 0 while n%i and i < n**0.5: i+=1 return i if n%i is 0 else 0 def prime_ant(n): A, p, i = {}, 0, 0 while i < n: if p not in A: A[p] = 2+p q = div(A[p]) if q is 0: p += 1 else: A[p] //= q if p-1 not in A: A[p-1] = 2+(p-1) A[p-1] += q p -= 1 i += 1 return pHi,
Can you explain how you came up with the wheel, and why it works?
Thanks
@B: See https://programmingpraxis.com/2009/05/08/wheel-factorization/.
@programmingpraxis, thanks!