获取具有不同部分的整数分区数的有效算法（分区函数 Q）

Question

我需要创建一个函数，它将接受一个参数int 并输出int，它表示输入整数分区的不同部分的数量。即，

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

因为我只寻找不同的部分，所以不允许这样的事情：

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

到目前为止，我已经设法提出了一些算法，可以找到所有可能的组合，但它们非常缓慢且仅在 n=100 左右之前有效。 而且由于我只需要 number 个组合而不是组合本身 Partition Function Q 应该可以解决问题。 有人知道如何有效地实现这一点吗？

有关问题的更多信息：OEIS、Partition Function Q

编辑：

为避免任何混淆，DarrylG 答案还包括琐碎（单个）分区，但这不会以任何方式影响其质量。

编辑 2： jodag（已接受的答案）不包括琐碎的分区。

Answer 1

测试了两种算法

1. 简单的递归关系

2. WolframMathword 算法（基于 Georgiadis、Kediaya、Sloane）

两者都使用 LRUCache 实现了记忆化。

结果：WolframeMathword 接近数量级的速度更快。

1.简单的递归关系（带记忆）

Reference

代码

@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

性能

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 算法

（基于 Georgiadis、Kediaya、Sloane）

Reference

代码

# 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

性能

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

结论：该算法比简单递归关系快几个数量级

算法

Reference

Answer 2

我认为解决这个问题的一种直接有效的方法是从原帖中的Wolfram PartitionsQ link 显式计算生成函数的系数。

这是一个非常说明性的示例，说明了如何构造生成函数以及如何使用它们来计算解。首先，我们认识到问题可能如下提出：

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.

我们可以通过构造以下多项式（即生成函数）找到Q(n)

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

解数是项组合成x^n的方式数，即展开多项式后x^n的系数。因此，我们可以通过简单的多项式乘法来解决这个问题。

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]

时间（平均 1000 次迭代）

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

您可以在您为 N 运行时中的二次方链接的 mathematica 文章中记住方程式 8、9 和 10 中的递归式。

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

工作由rec内部函数完成，它需要：

• part - 总和始终等于或小于目标 n 的部分列表
• tgt - 需要添加到 part 的总和才能得到 n 的剩余部分总和
• allowed - 一组仍然允许在完整分区中使用的数字

当tgt = 0 被传递时，这意味着part 的总和如果n 和part 被添加到结果列表中。如果tgt 仍然是正数，则在递归调用中尝试将每个允许的数字作为part 的扩展。