Two Factoring Games
February 18, 2011
We calculate the home prime of a number with the definitional recursion:
(define (home-prime n)
(let* ((fs (factors n))
(hp (undigits (apply append (map digits fs)))))
(if (prime? hp) hp (home-prime hp))))
A small variant of that function accumulates the intermediate steps:
(define (home-factors n)
(let loop ((n n) (hps '()))
(if (prime? n) (reverse (cons n hps))
(let* ((fs (factors n))
(hp (undigits (apply append (map digits fs)))))
(loop hp (cons fs hps))))))
Here are some examples:
> (home-prime 99)
71143
> (home-factors 8)
((2 2 2) (2 3 37) (3 19 41) (3 3 3 7 13 13) (3 11123771)
(7 149 317 941) (229 31219729) (11 2084656339)
(3 347 911 118189) (11 613 496501723) (97 130517 917327)
(53 1832651281459) (3 3 3 11 139 653 3863 5107)
3331113965338635107)
To calculate the terms of the Euclid-Mullin sequence, we use a variant of trial division that stops after the first found factor to search for small factors, then Pollard rho for larger numbers:
(define (td-least n)
(let loop ((p 3))
(cond ((< 1000000 p) #f)
((zero? (modulo n p)) p)
(else (loop (+ p 2))))))
(define (euclid-mullin)
(let loop ((n 1) (ems '(2)))
(display n) (display " ")
(display (car ems)) (newline)
(let* ((next (+ (apply * ems) 1))
(lpf (td-least next)))
(if lpf (loop (+ n 1) (cons lpf ems))
(let ((fs (factors next)))
(cond ((= n 16)
(loop 17 (cons 30693651606209 ems)))
((= n 27)
(loop 28 (cons 643679794963466223081509857 ems)))
((number? (car fs))
(loop (+ n 1) (cons (car fs) ems)))
((pair? (car fs))
(loop (+ n 1) (cons (caar fs) ems)))
(else (error 'euclid-mullin
"can't complete factorization"))))))))
The product sequence grows very quickly, so quickly that we had to “cheat” at the 16th and 27th terms. For instance, after the first ten factors, 2 × 3 × 7 × 43 × 13 × 53 × 5 × 6221671 × 38709183810571 × 139 + 1 = 208277726667994127981027430331 = 2801 × 2897 × 489241 × 119812279 × 437881957. Here’s our sample output:
> (euclid-mullin)
1 2
2 3
3 7
4 43
5 13
6 53
7 5
8 6221671
9 38709183810571
10 139
11 2801
12 11
13 17
14 5471
15 52662739
16 23003
17 30693651606209
18 37
19 1741
20 1313797957
21 887
22 71
23 7127
24 109
25 23
26 97
27 159227
28 643679794963466223081509857
29 103
30 1079990819
31 9539
32 3143065813
33 29
34 3847
35 89
36 19
37 577
38 223
39 139703
40 457
41 9649
42 61
43 4357
Error in euclid-mullin: can't complete factorization.
We used prime? from the exercise on the Baillie-Wagstaff primality checker, factors from the exercise on Richard Brent’s variant of John Pollard’s rho factorication algorithm, and the digits and undigits functions from the Standard Prelude; you will need some stronger factorization tool to compute further. You can run the program at http://programmingpraxis.codepad.org/TULr3faG.
Programming Praxis 
[…] today’s Programming Praxis exercise, our goal is to implement two algorithms related to prime factors; one […]
My Haskell solution (see http://bonsaicode.wordpress.com/2011/02/18/programming-praxis-two-factoring-games/ for a version with comments):
import Data.List import Data.Numbers.Primes homePrime :: Integer -> Integer homePrime = head . filter isPrime . iterate (read . concatMap show . primeFactors) euclidMullin :: [Integer] euclidMullin = unfoldr (\p -> let a = head (primeFactors $ p + 1) in Just (a, p * a)) 1My solution isn’t very long, but the machinery behind the scenes is rather
extensive. In an effort to save memory, I’ve tried to make heavy use of
iterators—sort of like lazy evaluation of lists in Haskell. Below is what I
wrote for this exercise; I also used an implementation of the Sieve of
Eratosthenes for primes and a Miller-Rabin primality test.
The full code is available here.
def times(x, n): t = 0 while not n % x: t += 1 n //= x return t from itertools import chain def prime_factors(n): pfs = () for p in primes(n): if not n % p: pfs = chain(pfs, [p for i in xrange(times(p, n))]) return pfs def home_prime(n): while not is_prime(n): n = int(''.join(str(p) for p in prime_factors(n))) return n def lpf(n): return next(prime_factors(n)) def euclid_mullin(): a = 2 p = 1 while True: yield a p *= a a = lpf(1 + p)