19
$\begingroup$

I need to answer the following for a number of parameters: How many ways can the integer $k$ be written as a sum of $n$ different integers ranging from $1$ to $m$?

My initial attempt was the following function:

NumberOfWays[k_, n_, m_] := 
  Count[Map[Length, 
    Map[DeleteDuplicates, 
     IntegerPartitions[k, {n}, Range[m]]]], 
   n];

This works, but becomes very slow as the parameters get big. I then thought I might do it using a generating function and attempted the following:

GenFuncy[m_] := Product[1 + y*x^j, {j, 1, m}];
NumberOfWays2[k_, n_, m_] := Coefficient[GenFuncy[m], x^k*y^n];

Again this works, but surprisingly (to me) it is even slower.

Is there any way I can speed these functions up, or maybe another faster way to do the calculation altogether?

$\endgroup$
  • 2
    $\begingroup$ Capitalized I is reserved for the imaginary unit, and N is a function meant for making exact numbers inexact; consider using different variables. $\endgroup$ – J. M. will be back soon Aug 14 '15 at 12:28
  • 1
    $\begingroup$ Anyway: look up QPochhammer[]. $\endgroup$ – J. M. will be back soon Aug 14 '15 at 12:40
  • 2
    $\begingroup$ Have you read this?: dlmf.nist.gov/26.10 $\endgroup$ – Mr.Wizard Aug 14 '15 at 13:10
  • 1
    $\begingroup$ @J. M. I must admit I'd like to see your solution. jorgen, this is not a good solution but at least it is faster than your NumberOfWays2: NumberOfWays3[k_, n_, m_] := {k, n} /. CoefficientRules[GenFuncy[m], {x, y}] $\endgroup$ – Mr.Wizard Aug 14 '15 at 13:32
  • 2
    $\begingroup$ Give SeriesCoefficient[ QPochhammer[-x y, x, #3], {x , 0, #1}, {y, 0, #2}] &[k, n, m] a whirl... a one-shot test using {k = 107, n = 11, m = 100} was nearly order of mag. faster than NumberOfWays3 and/or straight Coefficient on the expanded product (the latter two ~ in speed in q&d tests.) $\endgroup$ – ciao Aug 14 '15 at 20:41
18
$\begingroup$

This seems pretty quick, particularly on larger cases / larger k, e.g. 451, 29, 101 finishes in a few seconds on the loungebook.

N.B. - I have not tested this exhaustively, just thrown together from ideas...

If[Min[#3, #1 - Tr@Range@(#2 - 1)] < 0, 0, 
   SeriesCoefficient[QPochhammer[-x y, x, Min[#3, #1 - Tr@Range@(#2 - 1)]], 
                     {x , 0, #1}, {y, 0, #2}]] &[n, k, m]

UPDATE:

This seems to be very fast, particularly on larger cases. n.b.: posted with testing in progress, I'd like to prove correctness, but so far empirical testing matches prior methods, and appears faster than answers prior on large cases...

myDP[n_, k_, m_] := If[n < Binomial[k + 1, 2] || m < k, 0, 
                       SeriesCoefficient[QBinomial[m, k, q], {q, 0, n - Binomial[k + 1, 2]}]]

For a huge case of {n, k, m} = {5050, 100, 5050} this took a fraction of a second on the loungebook to return the result of 1 (for this case, there would be ~$2.74235\times 10^{68}$ partitions generated for any of the partition massaging methods like the OP's NumberOfWays, making use of these absurd for anything other than minimal cases.) The neat follow-up solution from KennyColnago took (unsure - aborted it after 5 minutes, monitoring progress indicated over an hour would be needed, figure 10X faster or so for both on a workstation...) for the same case - but I'd prefer to perhaps have his benchmark post extended with results on his hardware for a fair comparison.

Update 2: A further optimization, taking advantage of the symmetry of the gaussian polynomial:

myDPc[n_, k_, m_] := 
 Module[{mn = Binomial[k + 1, 2], mx = (k - k^2 + 2 k m)/2},
  If[mn <= n <= mx && m >= k,
     SeriesCoefficient[QBinomial[Min[n - Binomial[k, 2], m], k, q],
                       {q, 0, If[n > (mn + mx)/2, mx - n, n - mn]}],0]];

On an exhaustive search for all valid n for {k,m}={45,60} this was over 4X faster than myDP, and for large cases (e.g., {n,k,m}={18775, 50, 400} it was over 20,000X faster than myDP.

There's an additional optimization possible that might be advantageous when searching ranges of {n,k,m}: for any given {n,k,m}, by symmetry of the Q-Binomial, there's a dual of {n', k',m} where n' and k' are simple transformations of n and k that has precisely the same polynomial. Memoization on that can about double the performance for such searches.

Update 3 2015/08/20: Added an optimization (in edited myDPc above) for larger k, resulting in over 2 orders of magnitude performance boost to e.g. {n,k,m}={5100,100,5100} and about three orders of magnitude boost to {n,k,m}={12000,154,12000}.

I think I've run out of ideas...

$\endgroup$
  • $\begingroup$ Looks like I don't have to write an answer anymore... this nails the deal. :) $\endgroup$ – J. M. will be back soon Aug 15 '15 at 7:36
  • $\begingroup$ @J. M. You flatter me. Thanks! $\endgroup$ – ciao Aug 15 '15 at 10:42
  • $\begingroup$ @jorgen: Appreciate the accept, perhaps you did not see KennyColnago's follow-up answer? Unless you're specifically wanting to use GF (perhaps to generalize to some other cases), or need to do this for cases where the memory usage of his follow-up becomes prohibitive, it's a fast and clever answer. If this were my question, and the two mentioned items did not apply, that's the answer I'd accept... $\endgroup$ – ciao Aug 16 '15 at 22:44
  • 1
    $\begingroup$ @J. M.: Enough to make it painfully obvious my combinatorics has become very rusty, but, I've not seen this derivation for this partition case before. Might need to swing by Persi's office for a sanity check... $\endgroup$ – ciao Aug 17 '15 at 7:57
  • 1
    $\begingroup$ I wish I understood this method but it has gone way over my head. $\endgroup$ – Mr.Wizard Aug 20 '15 at 10:12
18
$\begingroup$

Here is a summary of comments (before @ciao's best answer above), with a change in notation. These functions calculate the number of partitions of n into exactly k distinct parts of size at most m.

NumberOfWays000[n_, k_, m_] := 
   Count[Map[Length,Map[DeleteDuplicates, IntegerPartitions[n,{k},Range[m]]]], k]

NumberOfWays001[n_, k_, m_] := 
   Total[Boole[Apply[Unequal, IntegerPartitions[n, {k}, Range[m]], 1]]]

NumberOfWays002[n_, k_, m_] := 
   Coefficient[Expand[Product[1 + t*z^j, {j, 1, m}]], t^k*z^n]

NumberOfWays003[n_, k_, m_] := 
   Coefficient[SeriesCoefficient[QPochhammer[-t*z, z, m], {z, 0, n}], t^k]

NumberOfWays004[n_, k_, m_] := 
   If[# =!= {n, k}, #, 0] &[{n, k} /. 
      CoefficientRules[Product[1 + t*z^j, {j, 1, m}], {z, t}]]

And some timings:

Block[{n = 50, kmax},
   kmax = Quotient[Sqrt[8 n + 1] - 1, 2];
   {AbsoluteTiming[Sum[NumberOfWays000[n, k, m], {k, 1, kmax}, {m, 1, n}]],
    AbsoluteTiming[Sum[NumberOfWays001[n, k, m], {k, 1, kmax}, {m, 1, n}]],
    AbsoluteTiming[Sum[NumberOfWays002[n, k, m], {k, 1, kmax}, {m, 1, n}]],
    AbsoluteTiming[Sum[NumberOfWays003[n, k, m], {k, 1, kmax}, {m, 1, n}]],
    AbsoluteTiming[Sum[NumberOfWays004[n, k, m], {k, 1, kmax}, {m, 1, n}]]}
]

(* {{2.462624, 109279}, {1.644974, 109279}, {19.891302, 109279},
    {9.002684, 109279}, {23.642686, 109279}} *)

A larger calculation with $n=107, k=11, m=100$ shows that methods based on IntegerPartitions slow dramatically to over 10 seconds. Methods 002, 003, 004 have timings of 1.5, 0.5, 1.8 s, respectively, subject to warnings about caching, etc. by ciao. Code in the answer by ciao returns in just 0.17 s.

$\endgroup$
  • $\begingroup$ Might want to add some test with large cases, since that's what op seems to be after, and any of the "counting" of massaged generated partitions that s/b fastest on small cases will utterly fail... $\endgroup$ – ciao Aug 14 '15 at 21:05
  • $\begingroup$ Also might want to add the expansion to method 001 as stated in comments... $\endgroup$ – ciao Aug 14 '15 at 21:13
  • $\begingroup$ @ciao Is method 004 the rewrite of method 001 to which your refer? A large example added and, as you suggest, the methods dependent on IntegerPartitions approach an utter fail... $\endgroup$ – KennyColnago Aug 14 '15 at 22:14
  • 1
    $\begingroup$ No, I meant taking the OP and going from Coefficient[GenFuncy[m] , x^k*y^n] to Coefficient[GenFuncy[m] // Expand, x^k*y^n] - that can speed things greatly. As an aside, use care in timings - MMA caches stuff from various coefficient & series functions, and I think some of it is opaque, that is, not cleared by ClearSystemCache etc. - I've found when benchmarking such things that running each case in a fresh kernel can lead to wildly differing results vs testing sequentially with ClearSystemCache interspersed.... $\endgroup$ – ciao Aug 14 '15 at 23:08
  • $\begingroup$ Ok, added method 001a with Expand, which makes a huge increase in speed, as you were aware. $\endgroup$ – KennyColnago Aug 14 '15 at 23:35
16
+100
$\begingroup$

Here is a totally different approach based on the fact that successive products forming the generating function are due to multiplication by a binomial $1+t*z^j$. Form a matrix $v$ of zeros with $n+1$ rows and $k+1$ columns. Initialize the top left corner to 1. Iterate $v=v+w$ where $w$ is the matrix $v$ shifted down by $j$ rows and to the right by 1. The element $v[[n+1,k+1]]$ after $m$ iterations is the number of partitions of $n$ into exactly $k$ parts of size at most $m$. It takes more memory, but is very fast.

DistinctPartitions2D[n_, k_, m_] :=
   Block[{n1 = n + 1, k1 = k + 1, v},
      v = ConstantArray[0, {n1, k1}];
      v[[1, 1]] = 1;
      Do[
         v = v + ArrayPad[v[[1 ;; n1 - j, 1 ;; k]], {{j,0}, {1,0}}],
         {j, 1, m}];
      v[[n1, k1]]
   ]

Some timings

AbsoluteTiming[DistinctPartitions2D[451, 29, 101]]
(* {0.010868, 231} *)

AbsoluteTiming[DistinctPartitions2D[3000, 30, 1729]]
(* {46.823481, 2704951781342880353088665158660429357} *)

Update

Clearly(!) @ciao's brilliant new code using QBinomial deserves his own bounty and the Accept. For large cases, it is orders of magnitude faster than DistinctPartitions2D. Nolo contendere. For example,

AbsoluteTiming[myDP[3000, 30, 1729]]

took about 1.5 s compared to 47 s for DistinctPartitions2D shown above. My method is slow because approximately 100000 array elements are added 1729 times, for about 161 million operations. ciao's case of

AbsoluteTiming[myDP[5050, 100, 5050]]

returned 1 in a fraction of a second, but DistinctPartitions2D would require about 2.6 billion operations. No need to benchmark a bunch of cases, math beats brute force every time. Well done.

$\endgroup$
  • $\begingroup$ Clever idea, +1 o/c. You can ~ double speed, mind if I add a derivation here or on mine? $\endgroup$ – ciao Aug 16 '15 at 3:51
  • $\begingroup$ @ciao I was secretly hoping you would have a go. Looking forward to it. $\endgroup$ – KennyColnago Aug 16 '15 at 15:43
  • 1
    $\begingroup$ @J. M.: End result that coeff. is taken from is not sparse, and you'd pay to overhead of updating the zeroes along the way...e.g., for {1000,10,1000} case, only ~11% zeroes in end result. We think alike, though.... $\endgroup$ – ciao Aug 17 '15 at 8:21
  • 1
    $\begingroup$ @ciao Updated this post with a couple of comparisons of your amazing QBinomial code. Never heard of QBinomial, don't understand it, I've got some homework to do. Thanks! $\endgroup$ – KennyColnago Aug 17 '15 at 17:47
  • 1
    $\begingroup$ 100 rep. bounty awarded, thanks again for prodding further thought on this interesting question. $\endgroup$ – ciao Aug 20 '15 at 3:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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