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.
In number theory and combinatorics, a partition of a positive integer n, also called an integer partition, is a way of writing n as a sum of positive integers. Two sums that differ only in the order of their summands are considered the same partition. (If order matters, the sum becomes a composition.)
How many partitions will be formed for the integer 3? Explanation: We need to find the combinations of positive integers which give 3 as their sum. These will be {3}, {2,1}, {1,1,1}. Thus the correct answer is 3.
Srinivasa Ramanujan first discovered that the partition function has nontrivial patterns in modular arithmetic, now known as Ramanujan's congruences. For instance, whenever the decimal representation of n ends in the digit 4 or 9, the number of partitions of n will be divisible by 5.
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With