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

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

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