Efficient algorithm for getting number of partitions of integer with distinct parts (Partition function Q)

Question

I need to create function which will take one argument int and output int which represents the number of distinct parts of input integer's partition. Namely,

input:3 -> output: 1 -> {1, 2}
input:6 -> output: 3 -> {1, 2, 3}, {2, 4}, {1, 5}
...

Since I am looking only for distinct parts, something like this is not allowed:

4 -> {1, 1, 1, 1} or {1, 1, 2}

So far I have managed to come up with some algorithms which would find every possible combination, but they are pretty slow and effective only until n=100 or so. And since I only need number of combinations not the combinations themselves Partition Function Q should solve the problem. Does anybody know how to implement this efficiently?

More information about the problem: OEIS, Partition Function Q

EDIT:

To avoid any confusion, the DarrylG answer also includes the trivial (single) partition, but this does not affect the quality of it in any way.

EDIT 2: The jodag (accepted answer) does not include trivial partition.

Answer 1

Tested two algorithms

1. Simple recurrence relation

2. WolframMathword algorithm (based upon Georgiadis, Kediaya, Sloane)

Both implemented with Memoization using LRUCache.

Results: WolframeMathword approach orders of magnitude faster.

1. Simple recurrence relation (with Memoization)

Reference

Code

@lru_cache(maxsize=None)
def p(n, d=0):
  if n:
    return sum(p(n-k, n-2*k+1) for k in range(1, n-d+1))
  else:
    return 1

Performance

n    Time (sec)
10   time elapsed: 0.0020
50   time elapsed: 0.5530
100  time elapsed: 8.7430
200  time elapsed: 168.5830

2. WolframMathword algorithm

(based upon Georgiadis, Kediaya, Sloane)

Reference

Code

# Implementation of q recurrence
# https://mathworld.wolfram.com/PartitionFunctionQ.html
class PartitionQ():
  def __init__(self, MAXN):
    self.MAXN = MAXN
    self.j_seq = self.calc_j_seq(MAXN)

  @lru_cache
  def q(self, n):
    " Q strict partition function "
    assert n < self.MAXN
    if n == 0:
      return 1

    sqrt_n = int(sqrt(n)) + 1
    temp = sum(((-1)**(k+1))*self.q(n-k*k) for k in range(1, sqrt_n))

    return 2*temp + self.s(n)

  def s(self, n):
    if n in self.j_seq:
      return (-1)**self.j_seq[n]
    else:
      return 0

  def calc_j_seq(self, MAX_N):
    """ Used to determine if n of form j*(3*j (+/-) 1) / 2 
        by creating a dictionary of n, j value pairs "
    result = {}
    j = 0
    valn = -1
    while valn <= MAX_N:
      jj = 3*j*j
      valp, valn = (jj - j)//2, (jj+j)//2
      result[valp] = j
      result[valn] = j
      j += 1

    return result

Performance

n    Time (sec)
10   time elapsed: 0.00087
50   time elapsed: 0.00059
100  time elapsed: 0.00125
200  time elapsed: 0.10933

Conclusion: This algorithm is orders of magnitude faster than the simple recurrence relationship

Algorithm

Reference

Answer 2

I think a straightforward and efficient way to solve this is to explicitly compute the coefficient of the generating function from the Wolfram PartitionsQ link in the original post.

This is a pretty illustrative example of how to construct generating functions and how they can be used to count solutions. To start, we recognize that the problem may be posed as follows:

Let m_1 + m_2 + ... + m_{n-1} = n where m_j = 0 or m_j = j for all j.

Q(n) is the number of solutions of the equation.

We can find Q(n) by constructing the following polynomial (i.e. the generating function)

(1 + x)(1 + x^2)(1 + x^3)...(1 + x^(n-1))

The number of solutions is the number of ways the terms combine to make x^n, i.e. the coefficient of x^n after expanding the polynomial. Therefore, we can solve the problem by simply performing the polynomial multiplication.

def Q(n):
    # Represent polynomial as a list of coefficients from x^0 to x^n.
    # G_0 = 1
    G = [int(g_pow == 0) for g_pow in range(n + 1)]
    for k in range(1, n):
        # G_k = G_{k-1} * (1 + x^k)
        # This is equivalent to adding G shifted to the right by k to G
        # Ignore powers greater than n since we don't need them.
        G = [G[g_pow] if g_pow - k < 0 else G[g_pow] + G[g_pow - k] for g_pow in range(n + 1)]
    return G[n]

Timing (average of 1000 iterations)

import time
print("n    Time (sec)")
for n in [10, 50, 100, 200, 300, 500, 1000]:
    t0 = time.time()
    for i in range(1000):
        Q(n)
    elapsed = time.time() - t0
    print('%-5d%.08f'%(n, elapsed / 1000))

n    Time (sec)
10   0.00001000
50   0.00017500
100  0.00062900
200  0.00231200
300  0.00561900
500  0.01681900
1000 0.06701700

Answer 3

You can memoize the recurrences in equations 8, 9, and 10 in the mathematica article you linked for a quadratic in N runtime.

Answer 4

def partQ(n):

    result = []

    def rec(part, tgt, allowed):
        if tgt == 0:
            result.append(sorted(part))
        elif tgt > 0:
            for i in allowed:
                rec(part + [i], tgt - i, allowed - set(range(1, i + 1)))

    rec([], n, set(range(1, n)))

    return result

The work is done by the rec internal function, which takes:

• part - a list of parts whose sum is always equal to or less than the target n
• tgt - the remaining partial sum that needs to be added to the sum of part to get to n
• allowed - a set of number still allowed to be used in the full partitioning

When tgt = 0 is passed, that meant the sum of part if n, and the part is added to the result list. If tgt is still positive, each of the allowed numbers is attempted as an extension of part, in a recursive call.