A variation of IntegerPartition?
IntegerPartitions[n, {3, 10}, Prime ~Array~ 10]
In Mathematica this will give a list of all the ways to get n as the sum of from three to ten of the first ten prime numbers, allowing duplicates as needed.
How can I efficiently find the sums that equal n, allowing each element to only be used once?
Using the first ten primes is only a toy example. I seek a solution开发者_如何转开发 that is valid for arbitrary arguments. In actual cases, generating all possible sums, even using polynomial coefficients, takes too much memory.
I forgot to include that I am using Mathematica 7.
The following will build a binary tree, and then analyze it and extract the results:
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]];
For example,
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}}}
This is not insanely fast, and perhaps the algorithm could be optimized further, but at least the number of partitions does not grow as fast as for IntegerPartitions
.
Edit:
It is interesting that simple memoization speeds the solution up about twice on the example I used before:
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}]};
Now,
In[118]:= getPartitions[200, Prime~Array~150] // Length // Timing
Out[118]= {0.219, 4660}
Can use Solve over Integers, with multipliers constrained between 0 and 1. I'll show for a specific example (first 10 primes, add to 100) but it is easy to make a general procedure for this.
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]
]
Out[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}}}
---edit--- To do this in version 7 one would use Reduce instead of Solve. I'll bundle this in one function.
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]]
Here is Leonid Shifrin's example:
Timing[Length[knapsack[200, Prime[Range[150]]]]]
Out[128]= {1.80373, 4660}
Not as fast as the tree code, but still (I think) reasonable behavior. At least, not obviously unreasonable.
---end edit---
Daniel Lichtblau Wolfram Research
I would like to propose a solution, similar in spirit to Leonid's but shorter and less memory intensive. Instead of building the tree and post-processing it, the code walks the tree and Sow
s the solution when found:
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]]
This code is slower than Leonid's
In[177]:=
UniqueIntegerParitions[200, Prime~Array~PrimePi[200]] //
Length // Timing
Out[177]= {0.499, 4660}
but uses about >~ 6 times less memory, thus allowing to go further.
精彩评论