10
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.

Mr.Wizard
  • 24,179
  • 5
  • 44
  • 125

3 Answers3

9

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}
Leonid Shifrin
  • 22,449
  • 4
  • 68
  • 100
  • This appears to work nicely. It does not include the parameter for set length, but I am going to play with it a bit. Thank you. – Mr.Wizard Feb 16 '11 at 19:14
  • I wasn't able to easily constrain the length. The problem is that solutions are accumulated from any sub-branch of the tree, not necessarily from the root. The best thing would be to somehow not produce lots of "false" branches at all, plus constrain the lengths of branches at the time we build the tree, but I could not figure out a good way to do this. – Leonid Shifrin Feb 16 '11 at 19:20
  • @Leonid Shifrin I will need some time to understand the method you provided. I would like to try to implement it myself, but do you believe this method can be extended to Rational arguments? – Mr.Wizard Feb 16 '11 at 19:30
  • I don't see any problem with that. I don't use any specific information about the numbers, so you don't have to change the code. Try, for example, `getPartitions[19/3, Range[10]/3]`. As for the method, the key part is that whenever we anticipate a partition candidate, `p` is used, while for the final number in a complete partition, `pf` is used as wrappers. The tricky part then is to reconstruct back the positions of preceding numbers in a given partition, from position of the final number. This is done by nextPosition - it moves "up" the tree, getting each time the position of previous number – Leonid Shifrin Feb 16 '11 at 19:34
  • Odd, I tried something like that but got a null set, so I thought it didn't work. Maybe my inputs were sloppy. – Mr.Wizard Feb 16 '11 at 19:50
  • To clarify a bit more the `nextPosition` (which should have been named `prevPosition`): we only have, as positions, 1-s and 2-s, since this is a binary tree. An alternative definition for `nextPosition` would be: `nextPosition[pos_List] := pos /. {prev___, 1, 2 .., 1 ..} :> {prev, 1, 1}`. This is easier to understand, but several times slower (which matters a lot for large number of partitions) - thus using `Compile` with a C target. The meaning is: go up, ignoring consecutive "left" branches, then the same with "right" branches, until we meet the left branch, then take its first leaf. – Leonid Shifrin Feb 16 '11 at 19:59
8

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

Daniel Lichtblau
  • 6,854
  • 1
  • 23
  • 30
  • +1 Nice trick to avoid the combinatorial explosion on the mask – Dr. belisarius Feb 16 '11 at 16:53
  • @Daniel Lichtblau I got errors trying this on version 7. Can it be made to work? – Mr.Wizard Feb 16 '11 at 19:05
  • 1
    @Daniel Your solutions are certainly simpler and more concise than mine. Interestingly, your second solution is also about 7 times faster than the first one (according to my measurements), on the same example with `200` and `Prime[Range[150]]`. So, apparently, `Reduce` is much faster than `Solve`, for this problem at least. – Leonid Shifrin Feb 17 '11 at 14:29
  • 1
    @Leonid, O Leonid...The speed gain lies not in Reduce, but in reduction. Specifically, I first remove the primes larger than the target. Reduce alone does not see fit to do that. – Daniel Lichtblau Feb 17 '11 at 15:47
  • @Daniel Oh, I see. I should have paid closer attention to your second code. In my case, reduction was easy, although my solution would benefit from such reduction as well - I only reduce based on the sum, which has to be re-computed for each partition candidate. – Leonid Shifrin Feb 17 '11 at 18:49
6

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 Sows 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.

Sasha
  • 5,935
  • 1
  • 25
  • 33