Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
1  
So, using RLink is not an option, I guess? –  Leonid Shifrin Jul 1 '13 at 20:27
1  
@LeonidShifrin A solution in Mathematica proper would be preferable. –  Sasha Jul 1 '13 at 21:04
add comment

2 Answers

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 *)
share|improve this answer
    
Yes, this is using ILP. –  Brett Champion Jul 1 '13 at 21:25
    
@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 Jul 2 '13 at 5:09
add comment

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.

share|improve this answer
    
For all the practical applications this is probably the way to go, but it does not take the second argument into account. –  Sasha Jul 1 '13 at 21:05
    
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 Jul 1 '13 at 21:29
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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