Next highly composite number?

16

1

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).

Sasha

Posted 2013-07-01T20:01:10.637

Reputation: 7 203

1So, using RLink is not an option, I guess? – Leonid Shifrin – 2013-07-01T20:27:41.810

3@LeonidShifrin A solution in Mathematica proper would be preferable. – Sasha – 2013-07-01T21:04:10.080

Answers

9

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 *)

kirma

Posted 2013-07-01T20:01:10.637

Reputation: 13 550

Yes, this is using ILP. – Brett Champion – 2013-07-01T21:25:04.310

@BrettChampion Yes, my first attempt was already ILP. In case someone wonders how your comment relates to my answer, the answer was extensively edited... – kirma – 2013-07-02T05:09:14.323

3

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.

bill s

Posted 2013-07-01T20:01:10.637

Reputation: 62 963

For all the practical applications this is probably the way to go, but it does not take the second argument into account. – Sasha – 2013-07-01T21:05:19.977

If you are using it for FFT length, why does it matter how the current number factors? It would appear that the only thing of importance is to get a highly composite number for the FFT calculation. – bill s – 2013-07-01T21:29:23.663

2

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 *)

KennyColnago

Posted 2013-07-01T20:01:10.637

Reputation: 14 269

1

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:

enter image description here

and

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

yields:

{48, 100, 108, 1024, 10125}

Similarly,

enter image description here

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] &].

ubpdqn

Posted 2013-07-01T20:01:10.637

Reputation: 53 491

Thank you. I should point few inefficiencies in your code. Note that FactorInteger[n] is called within the body of the pure function, i.e. your code factors n for 2^Length[p]-1 times. This can be inefficient. Perhaps more efficient would be fac2[n_, p_List] := With[{pr = FactorInteger[n][[All, 1]]}, Intersection[pr, p] =!= pr]; nc2[n_, p_] := NestWhile[# + 1 &, n, fac2[#, p] &]. Then your version nc[1230233, {2, 3, 5, 7, 11}] takes 0.4 seconds, while nc2[1230233, {2, 3, 5, 7, 11}] takes only 0.03, while giving the same answer of 1232000. – Sasha – 2015-04-18T04:20:35.020

@Sasha thank you for the instructive advice...my time limited aim was to reproduce the R implementation and agree entirely subset for large list of primes would be unwieldy but for specified purpose, most of time 3 primes, not too expensive...but accept entirely your more efficient code. – ubpdqn – 2015-04-18T04:24:28.833

@Sasha I will add as update with attribution to you – ubpdqn – 2015-04-18T04:29:52.677

@Sasha FYI nextn(1230233,factors=c(2,3,5,7,11)) in R also yields...1232000 but I assume you already know that – ubpdqn – 2015-04-18T04:59:48.613