# How to efficiently find all sets of primes that add to a given number?

For large numbers the naive approach falls down:

Select[Subsets[Range[1, 4]], PrimeQ@Total@# &]


{{2},{3},{1,2},{1,4},{2,3},{3,4},{1,2,4}}

Subsets[Range[1, 50]]


Subsets::toomany: "The number of subsets (1125899906842624) indicated by... is too large; it must be a machine integer"

Any suggestions?

• IntegerPartitions[#, All, Prime[Range@PrimePi@#]] & - or do you mean sets with distinct elements? – ciao Sep 7 '15 at 6:06
• @ciao sets in the mathy sense, so yes here sets have distinct elements. – M.R. Sep 7 '15 at 6:10
• You have 6000+ rep and still pasting images instead of code. OMG. – Dr. belisarius Sep 7 '15 at 6:11
• @M.R. In that case, With[{ip = IntegerPartitions[#, All, Prime[Range@PrimePi@#]] &@#}, Select[ip, Length@# == Length@DeleteDuplicates@# &]] & ? – ciao Sep 7 '15 at 6:17
• BTW, I note that what your title says, and what your example does, do not mesh - the example shows finding sets that sum to a prime. So, what is it you want, exactly? – ciao Sep 7 '15 at 8:21

I'll also assume the question wants sets of distinct primes summing to an integer n. Following on from the very informative discussion of distinct partitions here (thank you @ciao!), I include the generating function vector update method below.

SumDistinctPrimes[n_] :=
Block[{n1 = n + 1, v, t},
v = ConstantArray[0, n1];
v[] = 1;
Do[
v[[p + 1 ;; n1]] += ToString[p]*v[[1 ;; n1 - p]],
{p, Prime[Range[PrimePi[n]]]}];
Map[ToExpression,
Apply[List, Apply[List, Expand[v[[n1]]]], 1], {2}]]


Each generating function binomial 1+"p"*z^p is tagged by the unique string "p" corresponding to prime p, then Apply is used to change sums and products into lists of primes. Almost all the time is spent on Expand[v[[n1]]], which perhaps someone could speed up...

Timings become faster than @ciao's distPrimePart[n] as n increases beyond about 250. Has @ciao got another brilliant QBinomial approach here?

I'll assume the title reflects the goal, as opposed to the example. If this is not the case, comment and I'll delete.

distPrimePart[n_] := Module[{l = 1, p},
Join @@ Reap[While[(p =Prime@Range@(PrimePi@(n - Tr@Prime@Range@(l - 1)))) != {},
Sow[Select[IntegerPartitions[n, {l}, p], Length@# == Length@DeleteDuplicates@# &]];
l++]][[2, 1]]];


much more efficient than my comment solution. Example use:

distPrimePart
(* {{23, 2}, {17, 5, 3}, {13, 7, 5}, {13, 7, 3, 2}, {11, 7, 5, 2}} *)


I assume you know this is a subset-sum problem, making it $\scriptsize\mathcal{N}\mathcal{P}$-Complete, you're not going to do this "efficiently"...