Partition Numbers

April 15, 2011

This exercise is trickier than it looks. The problem is that the recursive calls to P(n) make the naive implementation take exponential time:

(define (partition n) ; naive version
  (if (negative? n) 0
    (if (zero? n) 1
      (let loop ((k 1) (sum 0))
        (if (< n k) sum
          (loop (+ k 1)
                (+ sum (* (expt -1 (+ k 1))
                          (+ (partition (- n (* k (- (* 3 k) 1) 1/2)))
                             (partition (- n (* k (+ (* 3 k) 1) 1/2))))))))))))

To calculate P(1000), we need to memoize the function by storing previously-computed values of the function. There are many ways to store those values, including hash tables, treaps, red-black trees, a simple, unbalanced binary search tree, a cuckoo hash, and functional arrays; for variety, we choose a growing arrays approach that keeps locally, inside the function closure, a vector that doubles in length whenever it becomes too small. Although it looks expensive to copy vector elements, timing tests show that it seems to work well enough, since the number of doublings is small:

(define partition ; memoized version
  (let ((len 1) (pv (vector 1)))
    (lambda (n)
      (do () ((< n len))
        (let* ((new-len (+ len len)) (new-pv (make-vector new-len #f)))
          (do ((i 0 (+ i 1))) ((= i len))
            (vector-set! new-pv i (vector-ref pv i)))
          (set! len new-len) (set! pv new-pv)))
      (cond ((negative? n) 0) ((zero? n) 1)
            ((and (< n len) (vector-ref pv n)) => (lambda (x) x))
            (else (let loop ((k 1) (sum 0))
                    (if (< n k) (begin (vector-set! pv n sum) sum)
                      (loop (+ k 1) (+ sum (* (expt -1 (+ k 1))
                        (+ (partition (- n (* k (- (* 3 k) 1) 1/2)))
                           (partition (- n (* k (+ (* 3 k) 1) 1/2))))))))))))))

The line (do () ((< n len)) (grow)) is the standard Scheme idiom for a while loop. The body of the while performs the growing function by creating a new vector, copying elements one-by-one from the old vector to the new vector, then updating len and pv to point to the new vector. The cond clause (and (< n len) (vector-ref pv n)) simultaneously determines if the partition number has already been calculated (because the vector is filled with #f when it is initialized) and returns its value if it has; the => operator passes on the value as the single argument of a local function, which in this case is the identity function (lambda (x) x) that passes on its input to its output. Here is partition in action:

> (partition 1000)
24061467864032622473692149727991

You can run the program at http://programmingpraxis.codepad.org/YO13k0sw.

About these ads

Pages: 1 2

19 Responses to “Partition Numbers”

  1. My Haskell solution (see http://bonsaicode.wordpress.com/2011/04/15/programming-praxis-partition-numbers/ for a version with comments):

    import Data.IntMap ((!), fromList, insert, findWithDefault)
    
    partition :: Int -> Integer
    partition x = if x < 0 then 0 else foldl p (fromList [(0,1)]) [1..x] ! x where
        p s n = insert n (sum [(-1)^(k+1) * (r pred + r succ) | k <- [1..n],
            let r f = findWithDefault 0 (n - div (k * f (3 * k)) 2) s]) s
    
    main :: IO ()
    main = print $ partition 1000 == 24061467864032622473692149727991
    
  2. [...] today’s Programming Praxis exercise, our goal is to write a function to calculate partition numbers and to [...]

  3. Dave Webb said

    A Python solution:

    class Memoize:
        """decorator to memoise a function"""
        def __init__(self, f):
            self.f = f
            self.cache = {}
    
        def __call__(self, *args):
            if not args in self.cache:
                self.cache[args] = self.f(*args)
            return self.cache[args]
    
    @Memoize
    def bigp(n):    
        if n < 0:
            return 0
        if n == 0:
            return 1
        
        # run k from n to 1 to avoid excessive recursion depth
        return sum((-1) ** (k + 1) * (bigp(n - k * (3 * k - 1) / 2) + bigp(n - k * (3 * k + 1) / 2)) for k in range(n,0,-1))
    
    print bigp(1000)
    
  4. arturasl said

    My non-recursive solution in Python:

    #!/usr/bin/python
    
    aMem = {0 : 1}
    aStackNeedto = []
    
    def Update(nParent, nSum):
    	if ((aStackNeedto[nParent]['Num'] * 2 - aStackNeedto[nParent]['Count'] - 1) // 2) % 2 != 0:
    		nSum *= -1
    
    	aStackNeedto[nParent]['Count'] += 1
    	aStackNeedto[nParent]['Sum']   += nSum
    	aStackNeedto.pop()
    
    def Partition(n):
    	aStackNeedto.append({'Num' : n, 'Count' : 0, 'Sum' : 0, 'Parent' : 0})
    
    	while aStackNeedto:
    		nNewParent = len(aStackNeedto) - 1
    		aCur       = aStackNeedto[nNewParent]
    
    		if aCur['Num'] < 0:
    			Update(aCur['Parent'], 0)
    		elif aCur['Num'] in aMem:
    			Update(aCur['Parent'], aMem[aCur['Num']])
    		elif aCur['Count'] == aCur['Num'] * 2:
    			aMem[aCur['Num']] = aCur['Sum']
    			Update(aCur['Parent'], aMem[aCur['Num']])
    		else:
    			for k in range(1, aCur['Num'] + 1):
    				aStackNeedto.append(
    				     {'Num' : (aCur['Num'] - (k * (3 * k - 1) / 2)), 'Count' : 0, 'Sum' : 0, 'Parent' : nNewParent}
    				)
    				aStackNeedto.append(
    				     {'Num' : (aCur['Num'] - (k * (3 * k + 1) / 2)), 'Count' : 0, 'Sum' : 0, 'Parent' : nNewParent}
    				)
    
    	return aMem[n]
    
    print Partition(1000)
    
  5. Graham said

    My Python solution, similar to Dave Webb’s:

    #!/usr/bin/env python
    
    from functools import update_wrapper
    
    def memoize(func):
        func.memo = {}
    
        def memoizer(arg):
            if arg in func.memo:
                return func.memo[arg]
            else:
                func.memo[arg] = result = func(arg)
                return result
    
        return update_wrapper(memoizer, func)
    
    
    @memoize
    def partitions(n):
        if n in xrange(11):
            return (1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42)[n]
        else:
            return sum(pow(-1, k + 1) * (partitions(n - k * (3 * k - 1) // 2) +
                partitions(n - k * (3 * k + 1) // 2)) for k in xrange(n, 0, -1))
    
    if __name__ == "__main__":
        print partitions(1000)
    
  6. Graham said

    A version that runs faster on my laptop, based on a formula given on the Mathworld page:

    #!/usr/bin/env python
    
    from functools import update_wrapper
    
    def memoize(func):
        func.memo = {}
    
        def memoizer(arg):
            if arg in func.memo:
                return func.memo[arg]
            else:
                func.memo[arg] = result = func(arg)
                return result
    
        return update_wrapper(memoizer, func)
    
    @memoize
    def sigma_1(n):
        return sum(filter(lambda k: n % k == 0, xrange(1, n + 1)))
    
    @memoize
    def partitions(n):
        if n in xrange(11):
            return (1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42)[n]
        else:
            return sum(sigma_1(n - k) * partitions(k) for k in xrange(n)) // n
    
    if __name__ == "__main__":
        print partitions(1000)
    
  7. slabounty said

    Here’s a ruby version. The PN array is used for caching. I actually tried using inject, but couldn’t get the caching working with it. I’d love to see someone else’s version with that in there.

    PN = []
    def p(n)
        return 1 if n == 0
        return 0 if n < 0
        return PN[n] if PN[n]
    
        sum = 0
        (1..n).each do |k|
            sum += ((-1)**(k+1)) * (p(n-(k*(3*k-1))/2) + p(n-(k*(3*k+1))/2))
            PN[n] = sum
        end
        sum
    end
    
    puts "Partition of 1000: #{p(1000)}"
    

    Also @Dave’s python version memoizing the function is pretty cool. I haven’t tried anything like that in Ruby. Anyone know if it’s doable?

  8. Here’s my Ruby implementation:

    $partitions = [1,1,2,3,5,7,11,15,22,30,42]

    def partition(n)
    return 0 if n < 0
    return $partitions[n] if $partitions[n] != nil

    partition = 0
    1.upto(n) do |k|
    partition += (-1) ** (k + 1) * (partition(n – (k * (3 * k – 1)) / 2) + partition(n – (k * (3 * k + 1)) / 2))
    end
    $partitions[n] = partition
    return $partitions[n]
    end

    p partition(1000)

    Here’s the output:

    24061467864032622473692149727991

  9. Graham said

    I was inspired by @Remco’s slick Haskell solutions to try my own. After little
    luck with IntMaps (perhaps due to overflow?) I came up with the following. Note:
    I’m indebted to @Remco for nearly all of this, either in this exercise or
    previous ones.

    import Data.List (foldl', nub, sort, subsequences)
    import Data.Map ((!), findWithDefault, fromList, insert)
    import Data.Numbers.Primes (primeFactors)
    
    sigma1 :: (Integral a) => a -> a
    sigma1 = sum . nub . map product . subsequences . primeFactors
    
    partitions :: (Integral a) => a -> a
    partitions n = foldl' p' (fromList [(0, 1)]) [1 .. n] ! n where
        p' m i = insert i (sum [sigma1 (i - j) * p'' j | j <- [0 .. i - 1],
            let p'' k = findWithDefault 1 k m] `div` i) m
    
    main :: IO ()
    main = print $ partitions 1000
    
  10. rohit said

    inspired by you graham

  11. Khanh Nguyen said
    open System
    open System.Collections.Generic
    open Microsoft.FSharp.Collections
    open Microsoft.FSharp.Math
    
    let memoize (f: 'a -> 'b) = 
        let dict = new Dictionary<'a,'b>()
    
        let memoizedFunc (input : 'a) =
            match dict.TryGetValue(input) with
            | true, x -> x
            | false, _ ->
                let ans = f input
                dict.Add(input, ans)
                ans
        memoizedFunc
    let rec P  = 
        let aux (n: int) = 
            match n with
            | x when x < 0 -> 0.0
            | 0 -> 1.0
            | _ -> 
                [1..n] 
                |> List.fold (fun acc k -> 
                                let a = n - (k*(3*k-1))/2
                                let b = n - (k*(3*k+1))/2
                                acc + float(pown -1 (k+1))*((P a) + (P b))) 0.0
        memoize aux
            
    P 1000
    
  12. Jussi Piitulainen said

    Straight Scheme with delayed and forced promises.

    (define (generate-vector n g)
      (let ((v (make-vector n)))
        (do ((k 0 (+ k 1))) ((= k n)) (vector-set! v k (g v k)))
        v))
    
    (define world
      (generate-vector
       1001
       (lambda (world n)  ;delayed numbers                                                                                       
         (delay
           (let sum ((k 1) (partial 0))
             (if (<= k n)
                 (sum (+ k 1) ((if (odd? k) + -)  ;                                                                              
                               partial
                               (+ (P (- n (* k (- (* 3 k) 1) 1/2)))
                                  (P (- n (* k (+ (* 3 k) 1) 1/2))))))
                 partial))))))
    
    (define (P n)         ;forces the delayed numbers                                                                            
      (cond ((negative? n) 0)
    	((zero? n) 1)
    	((positive? n) (force (vector-ref world n)))))
    
    (write (P 1000)) (newline)
    

    (The promise at 0 remains an empty sum. The rest are good.)

  13. Dan Harasty said

    Hey there! New to the site, but already loving it. Here’s my Python solution. Much like the above, with one VERY BIG time saving feature. (I’m not sure if the “regular crowd” here is in to optimizations… but this one seemed worthy of note.)

    @memoize
    def partition_number_opt1(n):
    if n < 0: return 0
    if n == 0: return 1
    sum = 0
    for k in range(1,n+1):
    # slight optimization:
    # once the inner expressions get below zero, then all remaining
    # partition numbers (in the rest of the sum) are zero,
    # thus they contribute nothing to the sum.
    # so: break out of the loop
    if n-(k*(3*k-1)//2) < 0: break
    sum += ((-1) ** (k+1)) * (partition_number_opt1(n-(k*(3*k-1)//2))+partition_number_opt1(n-(k*(3*k+1)//2)))
    return sum

  14. Dan Harasty said

    Just trying my post once again, after I “read the manual”. Mea culpa!

    @memoize
    def partition_number_opt1(n):
    	if n < 0: return 0
    	if n == 0: return 1
    	sum = 0
    	for k in range(1,n+1):
    		# slight optimization:
    		# once the inner expressions get below zero, then all remaining
    		# partition numbers (in the rest of the sum) are zero,
    		# thus they contribute nothing to the sum.
    		# so: break out of the loop
    		if n-(k*(3*k-1)//2) < 0: break
    		sum += ((-1) ** (k+1)) * (partition_number_opt1(n-(k*(3*k-1)//2))+partition_number_opt1(n-(k*(3*k+1)//2))) 
    	return sum
    
  15. Jussi Piitulainen said

    Using Dan’s observation as an excuse for a little refactoring that I wished
    I did in the first place: making sum a procedure of its own. And yes, it
    does feel faster with the early completion of the sums.

    (define (generate-vector n g)
      (do ((k 0 (+ k 1)) (v (make-vector n) v))
          ((= k n) v) (vector-set! v k (g v k))))
    
    (define (sum b e g)   ;sum (g k) from b until (e k)                                                                          
      (do ((k b (+ k 1)) (s 0 (+ s (g k)))) ((e k) s)))
    
    (define world
      (generate-vector
       1001
       (lambda (world n)  ;delayed numbers                                                                                       
         (delay
           (sum 1 (lambda (k) (< n (* k (- (* 3 k) 1) 1/2)))
                (lambda (k)
                  (* (if (odd? k) +1 -1)
                     (+ (P (- n (* k (- (* 3 k) 1) 1/2)))
                        (P (- n (* k (+ (* 3 k) 1) 1/2)))))))))))
    
    (define (P n)         ;forces the delayed numbers                                                                            
      (cond ((negative? n) 0)
            ((zero? n) 1)
            ((positive? n) (force (vector-ref world n)))))
    
    (write (P 1000)) (newline)
    
  16. Vikas Tandi said

    The recursive solution is very slow because redundant calculation of partition numbers.
    So, The trick is to remember the previously calculated value of partition numbers.
    Checkout my implementation in C using red-black trees.
    http://codepad.org/vKEivBGq

  17. Eric said
    def p(n, L = [[1]]):
        for i in range(1, n + 1):
            L.append([L[x][0] for x in range(i - 1, -1, -1)])
            for j in range((i + 2)/2, i):
                L[j][0] = L[j].pop(0) + L[j][0]
        return sum(L[n])
    
    print p(1000)
    
  18. Eric said

    Here is my non-recursive python solution, simple and efficient!
    Since Euler’s formula is complex for calculation, I used another algorithm,
    instead, which can easily be discovered from the process of partition.

    def p(n, L = [[1]]):
        for i in range(1, n + 1):
            L.append([L[x][0] for x in range(i - 1, -1, -1)])
            for j in range((i + 2)/2, i):
                L[j][0] = L[j].pop(0) + L[j][0]
        return sum(L[n])
    
    print p(1000)
    
  19. danaj said

    I found this quite interesting, as there is an enormous range of performance between different solutions. The Perl version below is combinatorial and, at least for me, is quite a bit faster and less memory than the previously shown ones (barring the C double solution which isn’t correct for large values). However it is about 300 times slower than the GMP combinatorial solution below. In turn, that is about 1,000 times slower than the Rademacher formula used in Pari or SymPy. Pari in turn is quite a bit slower than Jonathan Bober’s MPFR Rademacher solution. The current state of the art implementation seems to be Fredrik Johansson’s FLINT/MPFR/Arb implementation which looks to be over 1000x faster than Pari, computing p(10^12) in under 12s. The growth rates differ (especially between the combinatorial and Rademacher solutions), so the “x faster” is simplistic.

    In Perl, computes p(10000) in about 8s.

    #!/usr/bin/env perl
    use warnings;
    use strict;
    use Math::BigInt try => "GMP";
    
    my $n = shift || die "Usage: $0 <n>    Returns the integer partition number of <n>\n";
    print partitions($n), "\n";
    
    sub partitions {
      my $n = shift;
      return 1 if defined $n && $n <= 0;
      my $d = int(sqrt($n+1));
      my @part = (Math::BigInt->bone);
      my @pent = (1, map { (($_*(3*$_+1))>>1, (($_+1)*(3*$_+2))>>1) } 1 .. $d);
      foreach my $j (scalar @part .. $n) {
        my ($psum1, $psum2, $k) = (Math::BigInt->bzero, Math::BigInt->bzero, 1);
        foreach my $p (@pent) {
          last if $p > $j;
          if ((++$k) & 2) { $psum1->badd( $part[ $j - $p ] ); }
          else            { $psum2->badd( $part[ $j - $p ] ); }
        }
        $part[$j] = $psum1 - $psum2;
      }
      return $part[$n];
    }
    

    In C+GMP, computes p(200000) in about 5s. Compile with “gcc -O partitions.c -lgmp -lm”.

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <gmp.h>
    
    static void partitions(unsigned long n, mpz_t result)
    {
      if (n <= 0) {
        mpz_set_ui(result, 1);
      } else if (n <= 3) {
        mpz_set_ui(result, n);
      } else {
        unsigned long i, j, k, d = (unsigned long) sqrt(n+1);
        unsigned long* pent = (unsigned long*) malloc((2*d+2) * sizeof(unsigned long));
        mpz_t* part = (mpz_t*) malloc((n+1) * sizeof(mpz_t));
        mpz_t psum;
    
        pent[0] = 0;
        pent[1] = 1;
        for (i = 1; i <= d; i++) {
          pent[2*i  ] = ( i   *(3*i+1)) / 2;
          pent[2*i+1] = ((i+1)*(3*i+2)) / 2;
        }
        mpz_init_set_ui(part[0], 1);
        mpz_init(psum);
        for (j = 1; j <= n; j++) {
          mpz_set_ui(psum, 0);
          for (k = 1; pent[k] <= j; k++) {
            if ((k+1) & 2) mpz_add(psum, psum, part[ j - pent[k] ]);
            else           mpz_sub(psum, psum, part[ j - pent[k] ]);
          }
          mpz_init_set(part[j], psum);
        }
        mpz_clear(psum);
        free(pent);
        mpz_set(result, part[n]);
        for (i = 0; i <= n; i++)
          mpz_clear(part[i]);
        free(part);
      }
    }
    
    int main(int argc, char **argv)
    {
      unsigned long n;
      mpz_t res;
    
      if (argc != 2) {
        printf("Usage: %s <n>    Returns the integer partition number of <n>\n", argv[0]);
        exit(-1);
      }
      n = strtoull(argv[1], 0, 10);
      mpz_init(res);
      partitions(n, res);
      gmp_printf("%Zd\n", res);
      mpz_clear(res);
      exit(0);
    }
    

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 598 other followers

%d bloggers like this: