All possible solutions of an equations in R

Question

I want to write an R code to generate all distinct solution of the equation x1+x2+...+xm=n, where n is a positive integer and xi>=0, is a non-negative interger.

I have written an R code for this purpose, but I don't like its run time. Furthermore, it dose not work for relatively large values of m and n, like n=m=10. My code is as follows.

 A<-matrix(ncol=m,nrow=choose((n+m-1),(m-1)))
  a<-1
  for (h in 1:((n+1)^m)) {
        l <- rep(NA, m)
        for (il in 1:m) {
            l[il] <- (((h-1) %/% ((n+1)^(il-1))) %% (n+1)) + 1
            l[il]<-l[il]-1
        }
         if(sum(l)==n) { A[a,]<-l;a<-a+1}
    }
  A

Dose anybody know a more efficient way for finding the solutions of that equation? Thank you.

Answer 1

It seems like you are going through many more combinations of numbers than needed. As you put in your code, n + m - 1 choose m - 1 represents the number of solutions that will satisfy this sort of problem. Think about why this works as follows: we want to split up n into m groups. If we lined up our n elements side by side, we need to make m-1 chops to create m groups. We can think about each chop as it's own element. Let's use n=8, m=5 to think about it concretely. It's as if we need to divide up 8 a's into 5 groups (4 chops, represented by "/"). This is the same as if we had a row of 12 empty spaces, and we had to put 8 a's and 4 /'s in the row. A sample output could be:

a a / / a / a a a a / a

If we find the difference between the indices of the slashes and subtract one (and pretend there are 2 extra slashes on the outside indices of 0 and n+m), we can find each solution. For the example above, the 5 numbers would be 2 (3 - 0 - 1), 0 (4 - 3 - 1), 1 (6 - 4 - 1), 4 (11 - 6 - 1), and 1 (13 - 11 - 1). Therefore, if we generate all of the (n + m - 1) choose (m - 1) results, we can do simple subtraction, and don't need to worry about any extra combinations.

We can emulate this with the code too, using the combinations() function in the gtools library.

library(gtools)
nums <- rep(0, m)
combs <- combinations(n + m - 1, m - 1)
for(i in 1:choose(n + m - 1,m - 1)){
  for(j in 1:m){
    #edge case for beginning
    if(j == 1){
      nums[j] <- combs[i, j] - 1
    }
    #edge case for end
    else if(j == m){
      nums[j] <- n + m - combs[i, j-1] - 1
    }
    else {
      nums[j] <- combs[i, j] - combs[i, j-1] - 1
    }
  }
  A[i,] <- nums
}

If you want to further speed things up, you can put the inner for loop in a function that takes in a row of the combinations matrix, n, and m, and returns a single row of the solution output you want. Then, you can use an apply function on that function:

rowSummer <- function(spaces, size, max){
  nums <- rep(0, size)
  for(j in 1:size){
    if(j == 1){
      nums[j] <- spaces[j] - 1
    }
    else if(j == size){
      nums[j] <- max - spaces[j-1] - 1
    }
    else {
      nums[j] <- spaces[j] - spaces[j-1] - 1
    }
  }
  nums
}

combs <- combinations(n+m-1, m-1)
A <- apply(combs, 1, rowSummer, size=m, max=n+m)
t(A)

However, the original solution above already speeds things up a great deal by only going through the necessary combinations.