获取具有不同部分的整数分区数的有效算法(分区函数 Q)
Efficient algorithm for getting number of partitions of integer with distinct parts (Partition function Q)
我需要创建一个函数,它将采用一个参数 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
(接受的答案)不包括平凡的分区。
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
被传递时,如果n
,则意味着part
的总和,并且part
被添加到结果列表中。如果 tgt
仍然是正数,则在递归调用中尝试将每个允许的数字作为 part
的扩展。
您可以在您为 N 运行时中的二次方程链接的 Mathematica 文章中记住方程式 8、9 和 10 中的递归。
测试了两种算法
简单递归关系
WolframMathword 算法(基于 Georgiadis、Kediaya、Sloane)
Both implemented with Memoization using LRUCache.
Results: WolframeMathword approach orders of magnitude faster.
1.简单递归关系(带记忆)
代码
@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)
代码
# 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
Conclusion: This algorithm is orders of magnitude faster than the simple recurrence relationship
算法
我认为解决这个问题的一种直接有效的方法是根据原始 post.
中的 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
我需要创建一个函数,它将采用一个参数 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
(接受的答案)不包括平凡的分区。
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
被传递时,如果n
,则意味着part
的总和,并且part
被添加到结果列表中。如果 tgt
仍然是正数,则在递归调用中尝试将每个允许的数字作为 part
的扩展。
您可以在您为 N 运行时中的二次方程链接的 Mathematica 文章中记住方程式 8、9 和 10 中的递归。
测试了两种算法
简单递归关系
WolframMathword 算法(基于 Georgiadis、Kediaya、Sloane)
Both implemented with Memoization using LRUCache.
Results: WolframeMathword approach orders of magnitude faster.
1.简单递归关系(带记忆)
代码
@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)
代码
# 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
Conclusion: This algorithm is orders of magnitude faster than the simple recurrence relationship
算法
我认为解决这个问题的一种直接有效的方法是根据原始 post.
中的 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