Programming Praxis

There are definitely wrong ways and right ways to compute both the sopf and kappa functions.

The wrong way to write the sum of prime factors function is to compute the factors of n, remove duplicates, and sum the remaining factors. Instead, the sum of prime factors for a range of numbers can be computed by sieving; begin with an array initialized to zero, then sieve as in the normal Sieve of Eratosthenes, except that each time a prime is identified by a zero in the next unexamined array location the prime is added to the current array location and all its multiples.

The wrong way to write the kappa function is to simply write the indicated summation as a loop; that’s an exponential time algorithm due to the repeated recursive calls. Instead, cache the results in an array by computing the kappa in increasing order, so smaller included kappa are already known when they are needed; that’s a quadratic-time algorithm, as each kappa up to n must be computed, and each requires a loop on j.

Thus, the computation of the number of prime partitions is simple:

(define (prime-parts n)
  (let ((sopf (make-vector (+ n 1) 0))
        (kappa (make-vector (+ n 1) 0)))
    (vector-set! sopf 0 1)
    (vector-set! kappa 0 1)
    (do ((p 2 (+ p 1))) ((< n p))
      (when (zero? (vector-ref sopf p))
        (do ((i p (+ i p))) ((< n i))
          (vector-set! sopf i
            (+ (vector-ref sopf i) p)))))
    (do ((i 2 (+ i 1)))
        ((< n i) (vector-ref kappa n))
      (vector-set! kappa i
        (do ((j 1 (+ j 1))
             (s (vector-ref sopf i)
                (+ s (* (vector-ref sopf j)
                        (vector-ref kappa (- i j))))))
            ((= i j) (/ s i)))))))

We define the two vectors, handle the special cases at zero at both vectors, perform the sieving in the first outer do loop, and compute successive values of the kappa function in the second outer do loop. Here’s an example:

> (time (prime-parts 1000))
(time (prime-parts 1000))
    3 collections
    356 ms elapsed cpu time, including 0 ms collecting
    401 ms elapsed real time, including 0 ms collecting
    12455200 bytes allocated, including 12647680 bytes reclaimed
48278613741845757

You can run the program at http://programmingpraxis.codepad.org/7ziwpcfj. If you are going to be computing lots of prime partitions, it makes sense to store the two vectors instead of recomputing them each time you call the function.

The problem statement didn’t require it, but it is reasonable to provide a function that gives a list of the prime partitions of a number instead of just a count. That’s an exercise in dynamic programming. For any given number n, first find the partitions of 0, then the partitions of 1, then the partitions of 2, and so on until you reach n.

For each intermediate target, compute the primes not exceeding the target; for instance, if the intermediate target is 10, the primes not exceeding 10 are 2, 3, 5 and 7. Thus the partitions of 10 are each of the partitions of 10-2=8 with 2 added, plus each of the partitions of 10-3=7 with 3 added, plus each of the partitions of 10-5=5 with 5 added, plus each of the partitions of 10-7=3 with 7 added, with all duplicates eliminated. Here’s the code:

(define (primes n)
  (let ((bits (make-vector (+ n 1) #t)))
    (let loop ((p 2) (ps '()))
      (cond ((< n p) (reverse ps))
            ((vector-ref bits p)
              (do ((i (+ p p) (+ i p))) ((< n i))
                (vector-set! bits i #f))
              (loop (+ p 1) (cons p ps)))
            (else (loop (+ p 1) ps))))))

(define (set-cons x xs)
  (if (member x xs) xs
    (cons x xs)))

(define (prime-parts-list n)
  (let ((parts (make-vector (+ n 1) (list))))
    (vector-set! parts 0 (list (list)))
    (do ((i 2 (+ i 1)))
        ((< n i) (vector-ref parts n))
      (vector-set! parts i
        (let ((xs (list)))
          (do ((ps (primes i) (cdr ps))) ((null? ps) xs)
             (do ((yss (vector-ref parts (- i (car ps))) (cdr yss))) ((null? yss))
              (set! xs (set-cons (sort < (cons (car ps) (car yss))) xs)))))))))

To eliminate duplicates, we sort each partition into ascending order to establish a canonical form, then use set-cons to add the new partition only if it differs from those already present. Here’s an example:

> (prime-parts-list 11)
((11) (3 3 5) (2 2 7) (2 2 2 2 3) (2 2 2 5) (2 3 3 3))

Pages: 1 2