Programming Praxis

Partition Numbers

April 15, 2011 9:00 AM

The partition-number function P(n) gives the number of ways of writing n as the sum of positive integers, irrespective of the order of the addends. For instance P(4) = 5, since 4 = 4 = 3 + 1 = 2 + 2 = 2 + 1 + 1 = 1 + 1 + 1 + 1. Sloane’s A000041 gives the first ten partition numbers as 1, 2, 3, 5, 7, 11, 15, 22, 30, and 42; the numerous references to that sequence point to many fascinating facts about partition numbers, including their close association with pentagonal numbers. By convention, P(0) = 1 and P(n) = 0 for negative n. Partition numbers are normally calculated by the formula, which was discovered by Leonhard Euler:

$P(n) = \sum_{k=1}^n (-1)^{k+1} \left[ P\left( n - \frac{k\left(3k-1\right)}{2} \right) + P\left( n - \frac{k\left(3k+1\right)}{2} \right) \right]$

Your task is to write a function that computes P(n), and to calculate P(1000). When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Posted by programmingpraxis

Categories: Exercises

Tags:

« Older Newer »

22 Responses to “Partition Numbers”

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

By Remco Niemeijer on April 15, 2011 at 9:56 AM

[…] today’s Programming Praxis exercise, our goal is to write a function to calculate partition numbers and to […]

By Programming Praxis – Partition Numbers « Bonsai Code on April 15, 2011 at 9:56 AM

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)

By Dave Webb on April 15, 2011 at 10:19 AM

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)

By arturasl on April 15, 2011 at 4:13 PM

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)

By Graham on April 15, 2011 at 5:19 PM

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)

By Graham on April 15, 2011 at 5:32 PM

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?

By slabounty on April 15, 2011 at 8:59 PM
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

By Elvis Montero on April 15, 2011 at 11:31 PM

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

By Graham on April 16, 2011 at 3:26 AM

inspired by you graham

By rohit on April 16, 2011 at 9:58 AM

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

By Khanh Nguyen on April 16, 2011 at 5:13 PM

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

By Jussi Piitulainen on April 17, 2011 at 7:01 AM

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

By Dan Harasty on April 19, 2011 at 9:17 PM

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

By Dan Harasty on April 19, 2011 at 9:25 PM

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)

By Jussi Piitulainen on April 20, 2011 at 5:50 PM

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

By Vikas Tandi on May 7, 2011 at 7:51 AM

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)

By Eric on May 30, 2011 at 3:48 AM

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)
```
By Eric on May 30, 2011 at 4:51 AM

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);
}

By danaj on October 30, 2013 at 3:48 PM

Here is a short Haskell solution that implicitly memoizes without using any non-Prelude data structures, but achieves quite good performances:

integerPartitions :: [Integer]
integerPartitions =
let
p = go (let kg k = (2*k):k:kg (k+1) in 0:kg 1) [] 0

go :: [Int] -> [[Integer]] -> Int -> [Integer]
go ks a !c = s:(if c>0 then go ks b (c-1) else (let (d:kt)=ks in go kt b d))
where
(s,b) = sf a
sf ((a1:t1):(a2:t2):r) = let (t,tr)=sf r in (a1+a2-t,t1:t2:tr)
sf ((a1:t1):_) = if c>0 then (a1,t1:[]) else (a1+1,t1:p:[])
sf _ = if c>0 then (0,[]) else (1,p:[])
in
1:p

By Carl Edman (@CarlEdman) on August 9, 2015 at 1:49 PM
Forgive me for screwing up the formatting in the previous post. This is better, I hope: http://codepad.org/wprw8CLb . FWIW, this code calculates integerPartitions!!100000 in a little over 10 seconds on my ancient workstation.

By Carl Edman (@CarlEdman) on August 9, 2015 at 1:58 PM
[…] I’ve written a Haskell function to count the number of partitions of an integer – it’s basically an implementation of the same algorithm as this one, using Euler’s Pentagonal Number Theorem: […]

By Calculating the count of integer partitions of a number (uses stateful vectors internally) - PhotoLens on February 1, 2022 at 1:30 AM

Mobile Site | Full Site

Get a free blog at WordPress.com Theme: WordPress Mobile Edition by Alex King.

Programming Praxis

Partition Numbers

Share this:

Related

22 Responses to “Partition Numbers”

Leave a Reply