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.

About these ads

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 [...]

  2. 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)) 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)
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 615 other followers

%d bloggers like this: