## 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.

Pages: 1 2

### 3 Responses to “Two Factoring Games”

1. […] today’s Programming Praxis exercise, our goal is to implement two algorithms related to prime factors; one […]

```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)) 1
```
3. Graham said

My 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)
```