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

You can run the program at 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).

    My Haskell solution (see for a version with comments):

    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. Adam said

    Thanks a lot for this formula. It was a huge problem for me to understand the Euler transformation and I think that without Your formula I’d give up. :D

Leave a Reply

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

You are commenting using your 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


Get every new post delivered to your Inbox.

Join 678 other followers

%d bloggers like this: