# Product of set elements not exceeding a bound

I have a set of distinct positive integers l (if you wish, you may assume it is a list ordered by size), and a bound lim. In general, l has thousands of elements, but all are at most equal to lim. I would like to retrieve the set of all products of elements of l that do not exceed lim. That is, I'd like an efficient way of executing the following:

Select[ Times @@@ Subsets[lst] , #<=lim & ]


The only thing I have approaching a solution (other than the above, which is clearly not workable) is a rather complicated recursive implementation that is not very efficient.

• Are you able to give rough values for lim and the max integers in l? (Also I assume we assume the integers are all positive?) Jan 4, 2017 at 17:14
• Depending on l and lim, you might first Select[l, #<=lim], and then use your Subsets approach (so I second @Quantum_Oli 's question). Jan 4, 2017 at 18:25
• Yes, the integers are all positive and at most equal to lim. Jan 4, 2017 at 18:28
• For the use-cases you are interested in, is Length[l] typically only a bit smaller than lim (l is "dense") or is lim much, much larger (l is "sparse")? I.e. what is a typical value for Length[l]/lim? Jan 4, 2017 at 22:32
• l is definitely sparse. For example, when lim=10^8, Length[l]/lim is about 0.003. (A word: clearly this is part of another problem. I've chosen not to post the real problem because I'm trying to solve a Project Euler problem, and my solution path leads me in this direction. If it's wrong, it's wrong, but I want to find that out for myself.) Jan 4, 2017 at 22:57

Fun problem! Here is another recursive implementation. Note that I'm assuming that the input l is sorted and duplicate-free, like OP describes, i.e. Union[l] == l is True.

prods[l_, lim_] :=
Block[{list = l, len = Length@l, storage = ConstantArray[0, Length@l]},
Reap[
Do[storage[[1]] = list[[i]];
prodsRec[1, i, Quotient[lim, list[[i]]]]
, {i, len}
]
][[2, 1]]
]

prodsRec[factorInd_, listPos_, lim_] := Block[{j = listPos + 1},
Sow[storage[[1 ;; factorInd]]];
While[j <= len && list[[j]] <= lim,
storage[[factorInd + 1]] = list[[j]];
prodsRec[factorInd + 1, j, Quotient[lim, list[[j]]]];
j++
]
]


Looks more procedural than perhaps is usual around here, but should be pretty fast. Very little list manipulation means very little copying of lists due to Mathematica lists being immutable.

Testing one dense and one slightly sparse situation:

AbsoluteTiming[Length[prods[Range[100], 100]]]


{0.005645, 539}

rand = Union[RandomInteger[10000, 300]];
Length[rand]


299

AbsoluteTiming[Length[prods[rand, 10000]]]


{0.003750, 399}

An example closer to what OP describes in comments:

rand2 = Union@RandomInteger[10^8, 10^8*3/1000];
Length[rand2]


299 563

AbsoluteTiming[Length[bigExample = prods[rand2, 10^8]]]


{3.114989, 311 435}

Since the numbers in rand2 are sampled uniformly from Range[10^8], it's no surprise almost all products with two or more factors is waaaay larger than 10^8, so most "products" have only one factor:

Tally[Length /@ bigExample]


{{1, 299 563}, {2, 11 828}, {3, 44}}

It's much harder when the allowed factors are all "small":

AbsoluteTiming[Length[smallFactors = prods[Range[300], 10^5]]]


{14.771860, 1 470 607}

Tally[Length /@ smallFactors]


{{1, 300}, {2, 44 850}, {3, 344 016}, {4, 598 930}, {5, 384 643}, {6, 91 472}, {7, 6 345}, {8, 51}}

• Nice solution. It is indeed fast. Unfortunately, I am interested in a limit of 10^12, and even this solution isn't fast enough. It appears that there are no non-recursive clever solutions, though I'll wait another day to see. If that's true, I'll have to find another approach. Thanks (+1). Jan 4, 2017 at 23:32
• I guess one would need to know more about the elements of l in order to proceed, if there is a pattern etc. I would be surprised if the number of valid products isn't huge for 10^12, haha. Jan 4, 2017 at 23:48

a recursive implementation.

test[ list_, rem_ /; (Length@rem > 0), lim_] := Module[{next, r},
next = Select[ {Append[list, #], Complement[rem, {#}]} & /@ rem ,
(Times @@ #[[1]] <= lim) &];
(Sow[#[[1]]]; test[#[[1]], #[[2]], lim]) & /@ next];
Union[Sort /@ Reap[test[{}, Range[10], 10]][[2, 1]]]


{{1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2}, {1, 3}, {1, 4}, {1, 5}, {1, 6}, {1, 7}, {1, 8}, {1, 9}, {1, 10}, {2, 3}, {2, 4}, {2, 5}, {1, 2, 3}, {1, 2, 4}, {1, 2, 5}}

Select[Subsets[Range[10]], Times @@ # <= 10 &][[2 ;;]] == %


True

A bigger example, (This will break the Subsets approach )

 Union[Sort /@ Reap[test[{}, Range[100], 100]][[2, 1]]] //
Length // AbsoluteTiming


{2.31678, 539}