Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

In the following paper, D. Kane describes an algorithm for subset sum that runs in logspace:

I am trying to implement it

subsetsum[sum0_, set0_] := Module[{sum = sum0, set = set0},
c = 0;
csum = Abs[sum] + Total[set] + 1;
p = NextPrime[csum];
n = Length[set];
check = False;
While[c <= n,
Print["c=" <> ToString[c] <> "p=" <> ToString[p]];
If[IntegerQ[Sum[x^(-sum)*Product[(1 + x^(set[[i]])), {i, 1, n}], {x, 1, p - 1}]/p],
check = True; Break[];];
c = c + Floor[Log[2, p]]; p = NextPrime[p];

This initial value $C$ which seems to be the $C$th prime or the next prime after $C$ is given as:

$ C = \lvert B \rvert + \sum_{i=1}^n \lvert m_i \rvert +1 $

$B$ is the sum to which a subset of some set $\{m_1,m_2,m_3,\ldots,m_n\}$ sum to.

Furthermore, on a simple instance such as:

    set = {2, 3, 5, 7, 8, 9};
    sum = 9;
    subsetsum[sum, set] 

the program seems to stop after 2 loops returning False... which is wrong. It should be returning true...

Anyone has any ideas on this?

share|improve this question
I don't see anything unary in your inputs, which IIRC is the target of DK's algorithm. Ping him at Stanford if you need clarification on paper, nice guy and responsive. – ciao May 20 '14 at 8:52
I think you want the test to be If[!IntegerQ[...]] since that is when you have a "witness" to existence of one or more solutions. Also, strictly speaking, that algorithm calls for more space efficient arithmetic than is being done above. – Daniel Lichtblau May 20 '14 at 15:29

[Okay, I realize this isn't an answer but it's a sort'a neat way to get those coefficients.]

One can alternatively use the Cauchy Integral Formula to get the appropriate coefficient of the product polynomial that enumerates the possible subset sums. Below is a numerical implementation. A symbolic version can of course be coded, but it would in effect be computing the coefficients explicitly.

residue[ll_, n_] := Module[
  {x, logsum, pts, m = 1101, vals},
  logsum = Total[Log[1 + x^ll]] - n*Log[x];
  pts = N[1/2*Exp[2*Pi*I*Range[0, m - 1]/m], 200];
  vals = Exp[logsum /. Transpose[{Thread[x -> pts]}]];

Here is an example of usage.

vals = Union[NextPrime[RandomInteger[{1, 100}, 20]]]

(* Out[184]= {7, 11, 17, 29, 37, 41, 43, 53, 59, 67, 79, 89, 101} *)

This next function computes the coefficients explicitly so we can compare.

coeffpoly[ll_] := 
 Rest[Module[{x}, CoefficientList[Times @@ (1 + x^ll), x]]]


(* Out[198]= 38 *)

SO there are 38 ways to form the value 358 as a subset sum of vals. Now we do this using residue.

residue[vals, 358]

(* Out[203]= \ 38.0000000000000000000000000000000000000000000000000000000000000000000\ 0000000000000000000000 + 0.*10^-90 I *)

Almost seems like magic. The weaknesses are

(1) One cannot get away with m=1101 in general. The integral approximation could need arbitrarily many terms.

(2) One cannot rely on numerical precision of 100 in general. This too might have to be ratcheted arbitrarily high.

As subset sums are NP complete clearly one or the other will have to go beyond polynomially large in the size of the input in order to make this into a fully reliable method. I do not know if it is possible to get "usually" acceptable approximations in a way that stays in polynomial complexity. I suspect that, at least for the class of subset sums for which lattice reduction methods usually work, the approach above would not require exponential resources in the problem size.

share|improve this answer
Neat indeed! The numerical stability of complex-valued functions is always tricky... – user13675 May 27 '14 at 4:11

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.