Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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.

like image 875
Mr.Wizard Avatar asked Feb 16 '11 15:02

Mr.Wizard


People also ask

What is a partition of a set of integers?

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 form for the integer 3?

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.

What is Ramanujan partition theory?

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.


2 Answers

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}
like image 135
Leonid Shifrin Avatar answered Sep 29 '22 10:09

Leonid Shifrin


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

like image 45
Daniel Lichtblau Avatar answered Sep 29 '22 12:09

Daniel Lichtblau