3
$\begingroup$

As an example I want to find all Integers up to a given limit which only have certain prime factors. I can do it efficiently for the case of 2, 3, and 5 by

g = 100;
t1 = Table[{2^x*3^y*5^z}, {z, 0, Log[g]/Log[5]}, {y, 0, Log[g/(5^z)]/Log[3]}, {x, 0, Log[g/(5^z*3^y)]/Log[2]}]
t2 = Sort[Flatten[t1]]

But I am looking for a more general solution where I might have a set of k different factors. I tried Nest but didn' t get far. And I don't want to use a function like FactorInteger and then Select the good ones. I am looking for a solution combining depth of Table and a (maybe) recursive function of the limits.

$\endgroup$

2 Answers 2

3
$\begingroup$

Perhaps you are looking for smooth numbers. The function pSmoothOuter is quite fast, but requires memory as the upper bound m grows larger. Input pmax is the maximum prime allowed.

pSmoothOuter[pmax_Integer, m_] :=
   Block[{s},
      s = 2^Range[0, Log[2, m]];
      Do[
         s = Pick[s = Flatten[Outer[Times, s, p^Range[0, Log[p, m]]]], UnitStep[m - s], 1],
         {p, Prime[Reverse[Range[2, PrimePi[pmax]]]]}];
      Sort[s]]

For example,

pSmoothOuter[5,10^2]

{1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32,
36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100}

AbsoluteTiming[{Total[#], Length[#]} &[pSmoothOuter[30, 10^8]]]

{0.015348, {2364148327261, 88415}}

The following version allows an input list of primes p, not necessarily contiguous.

pSmoothOuter[p_List, m_] :=
   Block[{s},
      s = Min[p]^Range[0, Log[Min[p], m]];
      Do[
         s = Pick[s=Flatten[Outer[Times, s, q^Range[0, Log[q, m]]]], UnitStep[m - s], 1],
         {q, Most@Reverse[Sort[p]]}];
      Sort[s]]

For example,

pSmoothOuter[{5, 13, 11}, 800]

{1, 5, 11, 13, 25, 55, 65, 121, 125, 143, 169, 275, 325, 605, 625, 715}

$\endgroup$
1
  • $\begingroup$ Wow! Runs like crazy! $\endgroup$
    – user57467
    Commented Aug 22, 2020 at 8:53
2
$\begingroup$

This should work. The construction is unwieldy, and it can be stream-lined, but at least this works:

makePrimes[numPrimes_, bound_] := Module[
  {powers = (Prime@Range[numPrimes])^Array[x, numPrimes] // Reverse},
  Sort@Flatten@Table[
    Times @@ powers // Evaluate,
    Sequence @@ MapThread[{#1, 0, #2} &,
      {Reverse@Array[x, numPrimes], 
       Log[bound/Prepend[Most@Exp@Accumulate@Log@powers, 1]]/Reverse@Log@Prime@Range[numPrimes]}
    ] // Evaluate
  ]
]

Then

makePrimes[3, 100]
(* {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100} *)

is the same list as in the OP's example.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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