## Prime Partitions

### October 19, 2012

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

### 10 Responses to “Prime Partitions”

1. […] today’s Programming Praxis exercise, our goal is to calculate the number of prime partitions for a given […]

2. The formula in the exercise is wrong: the second call to sopf should be sopf(j), not sopf(n).

```import Data.List
import qualified Data.Map as M
import Data.Numbers.Primes

sopf :: Integer -> Integer
sopf = sum . nub . primeFactors

k :: Integer -> Integer
k n = snd \$ foldl' calc M.empty [1..n] M.! n where
calc dict i = M.insert i (s, div (s + sum (map (\j -> (fst \$ dict M.! j) *
(snd \$ dict M.! (i - j))) [1..i-1])) i) dict where s = sopf i

main :: IO ()
main = print \$ k 1000 == 48278613741845757
```
3. ardnew said

small typo: sopf(42)=12, not 7

4. programmingpraxis said

@ardnew: Fixed. Thanks.

5. Paul said

A version in Python. It uses memoization to speed up the calculations.

```from ma.primee import td_factors as factors

sopfs = {}
def sopf(n):
if n not in sopfs:
sopfs[n] = sum(set(factors(n)))
return sopfs[n]

kappas = {1:0}
def kappa(n):
if n not in kappas:
kmax = max(kappas.keys())
for i in range(kmax + 1, n + 1):
kappas[i] = (sopf(i) + sum(sopf(j) * kappas[i - j]
for j in range(1, i))) // i
return kappas[n]

print kappa(1000)
```
6. […] we’re back into the mathy sort of problems from Programming Praxis, tasked with calculating the number of prime partitions for a given number–essentially, how many different […]

7. JP said

Here’s my solution in Racket (along with a memoization macro that makes the runtime bearable):
Prime Partitions
“Memoization in Racket

( Also, there’s still a typo in the equation. That took me a fair bit to notice until I saw Remco Niemeijer’s comment. :) )

8. JP said

As a follow up, I wrote up a way to actually generate the prime partitions which might be interesting:
Prime Partitions II: The Listing

The interesting points at least from my perspective are that it uses generators and `for*/list` to control the recursion, making the code relatively clean and simple.

9. […] Pages: 1 2 […]