IntegerPartition的变体?
IntegerPartitions[n, {3, 10}, Prime ~Array~ 10]
在Mathematica中,这将列出所有的方法来获得n,作为前10个素数中三到十个的总和,允许根据需要重复。
我怎样才能有效地找到等于n的和,让每个元素只能使用一次 ?
使用前十个素数仅仅是一个玩具的例子。 我寻求一个对任意参数有效的解决方案。 在实际情况下,即使使用多项式系数产生所有可能的和,也会占用太多内存。
我忘记包括我正在使用Mathematica 7。
下面将构建一个二叉树,然后分析它并提取结果:
Clear[intParts];
intParts[num_, elems_List] /; Total[elems] < num := p[];
intParts[num_, {fst_, rest___}] /;
fst < num := {p[fst, intParts[num - fst, {rest}]], intParts[num, {rest}]};
intParts[num_, {fst_, rest___}] /; fst > num := intParts[num, {rest}];
intParts[num_, {num_, rest___}] := {pf[num], intParts[num, {rest}]};
Clear[nextPosition];
nextPosition =
Compile[{{pos, _Integer, 1}},
Module[{ctr = 0, len = Length[pos]},
While[ctr < len && pos[[len - ctr]] == 1, ++ctr];
While[ctr < len && pos[[len - ctr]] == 2, ++ctr];
Append[Drop[pos, -ctr], 1]], CompilationTarget -> "C"];
Clear[getPartitionsFromTree, getPartitions];
getPartitionsFromTree[tree_] :=
Map[Extract[tree, #[[;; -3]] &@FixedPointList[nextPosition, #]] &,
Position[tree, _pf, Infinity]] /. pf[x_] :> x;
getPartitions[num_, elems_List] :=
getPartitionsFromTree@intParts[num, Reverse@Sort[elems]];
例如,
In[14]:= getPartitions[200,Prime~Array~150]//Short//Timing
Out[14]= {0.5,{{3,197},{7,193},{2,5,193},<<4655>>,{3,7,11,13,17,19,23,29,37,41},
{2,3,5,11,13,17,19,23,29,37,41}}}
这不是非常快,也许算法可以进一步优化,但至少分区的数量不会像IntegerPartitions
那样快。
编辑:
有趣的是,简单的便笺化在以前使用的示例中将解决方案的速度提高了大约两倍:
Clear[intParts];
intParts[num_, elems_List] /; Total[elems] < num := p[];
intParts[num_, seq : {fst_, rest___}] /; fst < num :=
intParts[num, seq] = {p[fst, intParts[num - fst, {rest}]],
intParts[num, {rest}]};
intParts[num_, seq : {fst_, rest___}] /; fst > num :=
intParts[num, seq] = intParts[num, {rest}];
intParts[num_, seq : {num_, rest___}] :=
intParts[num, seq] = {pf[num], intParts[num, {rest}]};
现在,
In[118]:= getPartitions[200, Prime~Array~150] // Length // Timing
Out[118]= {0.219, 4660}
可以使用整数求解,乘数约束在0和1之间。我将展示一个特定示例(前10个素数,加上100),但很容易为此制作一个通用程序。
primeset = Prime[Range[10]];
mults = Array[x, Length[primeset]];
constraints01 = Map[0 <= # <= 1 &, mults];
target = 100;
Timing[res = mults /.
Solve[Flatten[{mults.primeset == target, constraints01}],
mults, Integers];
Map[Pick[primeset, #, 1] &, res]
]
出[178] = {0.004,{{7,11,13,17,23,29},{5,11,13,19,23,29},{5,7,17,19,23,29} ,2,5,11,13,17,23,29,{2,3,11,13,19,23,29},{2,3,7,17,19,23,29},{ 2,3,5,7,11,13,17,19,23}}}
---编辑---要在版本7中执行此操作,将使用Reduce而不是Solve。 我将把它捆绑在一个函数中。
knapsack[target_, items_] := Module[
{newset, x, mults, res},
newset = Select[items, # <= target &];
mults = Array[x, Length[newset]];
res = mults /.
{ToRules[Reduce[
Flatten[{mults.newset == target, Map[0 <= # <= 1 &, mults]}],
mults, Integers]]};
Map[Pick[newset, #, 1] &, res]]
这是列昂尼德希夫林的例子:
Timing[Length[knapsack[200, Prime[Range[150]]]]]
出[128] = {1.80373,4660}
没有树代码那么快,但仍然(我认为)合理的行为。 至少,并不明显不合理。
---结束编辑---
Daniel Lichtblau Wolfram Research
我想提出一个解决方案,与Leonid的精神类似,但是更短,内存更少。 该代码不是构建树并对其进行后处理,而是在发现时Sow
树并找到解决方案:
Clear[UniqueIntegerParitions];
UniqueIntegerParitions[num_Integer?Positive,
list : {__Integer?Positive}] :=
Block[{f, $RecursionLimit = Infinity},
f[n_, cv_, {n_, r___}] :=
(Sow[Flatten[{cv, n}]]; f[n, cv, {r}];);
f[n_, cv_, {m_, r___}] /; m > n := f[n, cv, {r}];
f[n_, cv_, {m_, r___}] /;
Total[{r}] >= n - m := (f[n - m, {cv, m}, {r}]; f[n, cv, {r}]);
f[___] := Null;
Part[Reap[f[num, {}, Reverse@Union[Cases[list, x_ /; x <= num]]]],
2, 1]]
这段代码比列昂尼德慢
In[177]:=
UniqueIntegerParitions[200, Prime~Array~PrimePi[200]] //
Length // Timing
Out[177]= {0.499, 4660}
但是使用大约6倍以下的内存,从而允许更进一步。
链接地址: http://www.djcxy.com/p/35565.html