Partition Numbers

April 15, 2011

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.

Pages: 1 2

18 Responses to “Partition Numbers”

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

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)