All possible combinations of x1 + ... + xk = s

Questions
,
1

I have a vague memory of having done this some time ago but I don’t remember the formula. What I need is all possible combinations of positive integers xi (i=1,..,k) that yield s when they are summed. For instance, for k=2 and s=2 there are the 3 combinations {0,2}, {1,1} and {2,0} the first number in a bracket representing x1 and the second x2. This can be done with the code

f[k_, s_] := Module[{c, x},
  c = 0;
  Do[
   If[Sum[x[i], {i, 1, k}] == s, c++],
   Evaluate[Apply[Sequence, Table[{x[i], 0, s}, {i, 1, k}]]]];
  Return[c]]

But I need a fast function as I will calculate it for very big k and s. So this code is too slow. Do you remember the formula involving factorials? Or is there a Mathematica function that does this?

2

Very big k and s worries me and positive versus zero puzzles me a little, but is IntegerPartitions the function that you might vaguely remember?

3

Hey Bill, IntegerPartitions is almost what I am looking for but I need only the number of partitions not the partitions themselves. For the example I cited above the number is 3 as there are the partitions {0,2}, {1,1} and {2,0} for k=2 and s=2. For k=2 and s=3 it is 4 and for k=3 and s=2 it’s 6 and so on. Of course, I could count the elements with Length[IntegerPartitions[s,{k}]] but I think there is a direct formula if I remember well. Unfortunately I don’t find it… And apparently IntegerPartitions does not include the number 0.

4

Ah, I didn’t understand that. Is it perhaps PartitionsP that you were thinking of? And A000041 - OEIS for a lot of background information on the subject.

5

No this cannot be PartitionsP because this function only takes one argument. It must take k and s.

6

You say that you want positive integers but your example has zeros. Would you clarify?

7

Well I assumed that 0 is a positive integer. Sorry if I am wrong… To be more precise, k>1 and s>1 but
xk>=0 for every k. By the way, it seems that the code above is not correct for some low values. I will be back tomorrow. Here in Europe it’s time to go sleep now…

8

Yes, formula exist and is very simple:

 myF[k_, s_]:=Pochhammer[k, s]/s! (*for : k >= s  *)

Regards M.I.

9

Wow this works. Thanks Mariusz!
Meanwhile, I also found the recursive formula

m[2, s_] := 1 + s; m[k_, 2] := m[k - 1, 2] + k;
m[k_, s_] := m[k, s] = m[k, s - 1] + m[k - 1, s];

But your function is better. By the way, in my code above for f[k,s] there is an error: the xi do not go from 0 to k but from 0 to s… I will edit it if I can.

10

Perhaps people might be more familiar with the (equivalent!) expression in terms of binomial coefficients:

myF[k_, s_] := Binomial[s + k - 1, s]
11

Hey Mariusz, I still have a question. In fact I realized that what I need is the values of this function, which is almost the same as the one in my first post expect that the xi go from 0 to d (where d>1 is a positive integer) instead of going from 0 to s:

f[k_, d_, s_] := Module[{c, x},
   c = 0;
   Do[
    If[Sum[x[i], {i, 1, k}] == s, c++],
    Evaluate[Apply[Sequence, Table[{x[i], 0, d}, {i, 1, k}]]]];
   Return[c]];

I guess this is somewhat more complicated. Do you know the formula for this? Or someone else?

Merry Christmas

12

Thanks J.M. Do you also know a book where is explained why this is so?

13

In WL both Binomial and Pochhammer are based on the Gamma function. Can be seen by evaluating:

Binomial[n, r] // FunctionExpand
Pochhammer[n, r] // FunctionExpand

For positive integers \Gamma (x)=(x-1)!
Use this in a replacement rule:

gRule = Gamma[x_] -> Factorial[x - 1];
FunctionExpand[Binomial[k + s - 1, s]] /. gRule
FunctionExpand[Pochhammer[k, s]/s!] /. gRule

Both cases will evaluate to \frac{(k+s-1)!}{(k-1)! s!}

14

In case anyone is interested: the analytical function for the numerical code above is

f[k_, d_, s_] := 
 Sum[(-1)^i Binomial[k, i] Binomial[k + s - 1 - (d + 1) i, k - 1], {i,
    0, Floor[s/(d + 1)]}]

And a faster algorithm is this one:

f[k_, d_, s_] := Module[{c, x},
   c = 0; x = 0;
   Do[
    If[Total[IntegerDigits[x, d + 1]] == s, c++];
    x++,
    {(d + 1)^k}];
   Return[c]];

And happy new year 2020.