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
2  
@LeonidShifrin A solution in Mathematica proper would be preferable. –  Sasha Jul 1 '13 at 21:04

3 Answers 3

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

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

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] &].
share|improve this answer
    
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 Apr 18 at 4:20
    
@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 Apr 18 at 4:24
    
@Sasha I will add as update with attribution to you –  ubpdqn Apr 18 at 4:29
    
@Sasha FYI nextn(1230233,factors=c(2,3,5,7,11)) in R also yields...1232000 but I assume you already know that –  ubpdqn Apr 18 at 4:59

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.