Next highly composite number?

Question

R language has this function 'nextn' (link) which computes the next highly composite number greater than a given one, which is used to find the optimal padding size for the subsequent FFT operation.

So I am looking for similar Mathematica function

NextHighlyCompositeNumber[ n_Integer?Positive, primes:{__Integer?PrimeQ} ]

which will output the smallest number $m \geqslant n$ such that $m = p_1^{e_1} \cdots p_k^{e_k}$ for some non-negative exponents $e_k$.

This seems to be a combinatorial search kind of problem, which could be solved using integer linear programming (if available).

Answer 1

Minimization can be expressed with plain NMinimize in a bit awkward manner:

NextHighlyCompositeNumber[n_Integer?Positive, primes:{__Integer?PrimeQ}] :=
  With[{parms = Unique[] & /@ primes}, 
    With[{eqn = Inner[Power, primes, parms, Times]}, 
      Floor@First@NMinimize[{eqn, eqn >= n &&
        parms \[Element] Integers && And @@ (# >= 0 & /@ parms)}, parms]]]

AbsoluteTiming@Table[NextHighlyCompositeNumber[n, {2, 3, 5}], {n, 1, 10}]

(* {3.063287, {1, 2, 3, 4, 5, 6, 8, 8, 9, 10}} *)

Alternatively, a shorter but more explicit exposition is possible using LinearProgramming:

NextHighlyCompositeNumber[n_Integer?Positive, primes:{__Integer?PrimeQ}] := 
  Quiet@Inner[Power, primes, 
    LinearProgramming[Log[primes], {Log[primes]}, {Log[n]}, Automatic, Integers],
    Times]

AbsoluteTiming@Table[NextHighlyCompositeNumber[n, {2, 3, 5}], {n, 1, 50}]

(* {0.140417, {1, 2, 3, 4, 5, 6, 8, 8, 9, 10, 12, 12, 15, 15, 15, 16, 18,
    18, 20, 20, 24, 24, 24, 24, 25, 27, 27, 30, 30, 30, 32, 32, 36, 36,
    36, 36, 40, 40, 40, 40, 45, 45, 45, 45, 45, 48, 48, 48, 50, 50}} *)

Idea in the latter case is to turn product-of-powers representing all possible values representable by these primes into sum-of-products by applying logarithm on all constraints. This makes the problem a dot product which is, naturally, amenable to linear programming.

Explicit LinearProgramming implementation (or variant of NMinimize that would have been written to the same form) is dramatically faster and numerically easier to evaluate.

EDIT:

Not so surprisingly, generating consecutive highly composite numbers is much faster than integer linear programming approach. This tail-recursive code achieves it (sorry, it's bit of a mess, nowhere near elegant):

HighlyCompositesTo[n_Integer?Positive, primes:{__Integer?PrimeQ}] :=
  HighlyCompositesTo[n, primes, {1}, 1 & /@ primes, primes, Min[primes]]

HighlyCompositesTo[n_, primes_, list_, pos_, next_, min_] := list /; Last[list] >= n

HighlyCompositesTo[n_, primes_, list_, pos_, next_, min_] := 
  With[{nlist = Append[list, min],
    npos = MapThread[If[#1 == min, #2 + 1, #2] &, {next, pos}]}, 
    With[{nnext = MapThread[If[#1 == min, #3 nlist[[#2]], #1] &, {next, npos, primes}]},
      HighlyCompositesTo[n, primes, nlist, npos, nnext, Min[nnext]]]]

This performs much better than NextHighlyCompositeNumber in extracting values up to 10000:

HighlyCompositesTo[10000, {2, 3, 5, 7}]; // AbsoluteTiming

(* {0.007496, Null} *)

NestWhileList[NextHighlyCompositeNumber[# + 1, {2, 3, 5, 7}] &,
  1, # < 10000 &]; // AbsoluteTiming

(* {5.247371, Null} *)

Both produce the same result:

NestWhileList[NextHighlyCompositeNumber[# + 1, {2, 3, 5, 7}] &, 1, # < 10000 &] ==
  HighlyCompositesTo[10000, {2, 3, 5, 7}]

(* True *)

Answer 2

One way to solve this problem is to make a list and then just choose from the list. The first 1200 highly composite numbers are given here. The largest one is order 10^87 so it is unlikely that you would need to go higher than this when doing a FFT. Here is a stripped down version:

highlyComposite[n_] := Module[{comps = {2, 4, 8, 16, 32, 64, 120, 128, 180, 240, 
  256, 360, 512, 720, 840, 1024, 1260, 1680, 2048, 2520, 4096, 5040, 7560, 
  8192, 10080, 15120, 16384, 20160, 25200, 27720, 32768, 45360, 
  50400, 55440, 65536, 83160, 110880, 131072, 166320, 221760, 
  262144, 277200, 332640, 498960, 524288, 554400, 665280, 720720, 
  1048576, 1081080, 1441440, 2097152, 2162160, 2882880, 3603600, 
  4194304, 4324320, 6486480, 7207200, 8388608, 8648640, 10810800, 
  14414400, 16777216, 17297280, 21621600, 32432400, 33554432, 
  36756720, 43243200, 61261200, 67108864, 73513440, 110270160, 
  122522400}}, 
  First[Select[comps, # > n &]]]

You can use this to find the next larger (highly composite number) than n. For instance, highlyComposite[121] returns 128 and highlyComposite[10000] gives 10080.

Answer 3

The numbers you seek are a generalized form of smooth numbers, composed of an arbitrary list of primes rather than simply a list of the first few primes. Dijkstra's algorithm, given here, generates numbers, in increasing order up from 1, which contain only non-negative powers of primes less than a maximum. Translated and tweeked for your problem, the code becomes the following.

NextHighlyComposite[n_Integer, p_List] :=
   Block[{i = Length[p], r, c, v, s, m = {1}},
         r = Range[i];
         c = ConstantArray[1, i];
         v = c;
         While[
               (s = Min[v = p*m[[c]]]) < n,
               c[[Pick[r, v, s]]] += 1;
               m = Flatten[{m, s}]];
         s = Min[v = p*m[[c]]];
         c[[Pick[r, v, s]]] += 1;
         s]

This code is very fast, not slowed by FactorInteger or NMinimize, and generalizes the approach by @bills to include an arbitrary list of primes.

Map[NextHighlyComposite[#, {2, 3, 5}] &, {47, 100, 101, 1001, 10001}]
(* {48,100,108,1024,10125} *)

NextHighlyComposite[1017, {7, 11, 13}]
(* 1183 *)

Answer 4

From the definition in R and testing:

fac[n_, p_] := 
 Or @@ (FactorInteger[n][[All, 1]] == # & /@ (Rest@Subsets[p]))
nc[n_, p_] := NestWhile[# + 1 &, n, ! fac[#, p] &]

Some R examples:

and

nc[#, {2, 3, 5}] & /@ {47, 100, 101, 1001, 10001}

yields:

{48, 100, 108, 1024, 10125}

Similarly,

nc[103, {3, 5, 7}]
nc[1017, {7, 11, 13}]

yield 105 and 1183 respectively.

UPDATE

See comments. Sasha has instructively pointed out code inefficiency. Sasha's code follows:

fac2[n_, p_List] := With[{pr = FactorInteger[n][[All, 1]]}, Intersection[pr, p] =!= pr]; 
nc2[n_, p_] := NestWhile[# + 1 &, n, fac2[#, p] &].