Find number of solutions of a linear equation of N variables

We have to find the number of solutions to a linear equation of N variables. We have solved this using Dynamic Programming.

For the sake of simplicity, let us assume that we have to find the non-negative integral solutions for the equation, and all the coefficients on the equation are positive integers.

As an example, let us consider the linear equation 3x₁ + x₂ = 9. Note that this is a linear equation of 2 variables. So, here, n = 2. The non-negative integral solutions for this equation in the form (x₁, x₂) are: (0, 9), (1, 6), (2, 3) and (3, 0). Hence, there are 4 solutions for this linear equation of 2 variables.

Input format: The input by the user consists of 3 lines: first line contains n (the number of variables), second line contains the coefficients of variables in LHS, and the third line contains the constant term of the equation (the RHS).
Example input:

2
3 1
9

This input refers to the equation 3x₁ + x₂ = 9.

Another example input:

4
2 4 5 1
23

This input refers to the equation 2x₁ + 4x₂ + 5x₃ + x₄ = 23.

Brute force approach:

We will approach the problem first using the naive method. This method involves trying out all the possible solutions, recursively, starting from 0.

The idea is to subtract the first coefficient from the constant value as many times as we can, then subtract the second coefficient from the remainder, and so on...

Recursion can be done as follows: once we subtract a coefficient from the constant, we just have to apply the same function, albeit this time, the constant will be the new value, i.e. the value we got after subtracting the coefficient from the original constant. The pseudo-code for the function can be written as:

numberOfSolutions(coefficients, startingIndex, endingIndex, constant)
    if constant == 0:
        return 1
    result = 0
    for i = start to i = endingIndex:
        if coefficients[i] <= constant:
            result += numberOfSolutions(coefficients, i, endingIndex, constant - coefficients[i])
    return result

The idea behind this pseudo-code is the same as what has already been described.

Hence, the C++ code for this function can be written as:

#include <iostream>

using namespace std;

int numberOfSolutions(int coefficients[], int startingIndex, int endingIndex, int constant) {
    if (constant == 0) {
        return 1;
    }

    int result = 0;

    for (int i = startingIndex; i <= endingIndex; i++) {
        if (coefficients[i] <= constant) {
            result += numberOfSolutions(coefficients, i, endingIndex, constant - coefficients[i]);
        }
    }

    return result;
}

int main (void) {
    int n;
    cin >> n;
    int coefficients[n];
    for (int i = 0; i < n; i++) {
        cin >> coefficients[i];
    }
    int constant;
    cin >> constant;
    int noOfSolutions = numberOfSolutions(coefficients, 0, n - 1, constant);
    cout << noOfSolutions << "\n";
    return 0;
}

The above implementation works, but in exponential time, since we are looping through all the possible integral values until we get a solution.

Hence, this method is not much efficient, and turns out, we can solve this problem in pseudo-polynomial time (the time complexity depends on the numeric value of the input), using dynamic programming.

Using Dynamic Programming:

In Dynamic Programming, we solve problems by combining the solutions to subproblems ("Programming" in this context refers to a tabular method, not to writing computer code). Dynamic programming applies when the subproblems overlap, that is, when subproblems share subproblems. A dynamic programming algorithm solves each subproblem just once and then saves its answer in a table, thereby avoiding the work of recomputing the answer every time it solves each subproblem.

The idea behind this approach is the same as the brute force method, but here, instead of calling the function recursively every time we subtract a coefficient from the constant, we save the constant values after every step in a table(typically an array), and hence, the last value in the table will be the number of solutions of the equation.

The C++ code for this approach is as follows:

#include <iostream>

using namespace std;

int numberOfSolutionsDP(int coefficients[], int n, int constant) {
    int table[constant + 1];
    for (int i = 0; i < constant + 1; i++) {
        table[i] = 0;
    }
    table[0] = 1;

    for (int i = 0; i < n; i++) {
        for (int j = coefficients[i]; j <= constant; j++) {
            table[j] += table[j - coefficients[i]];
        }
    } 
    return table[constant];
}

int main (void) {
    int n;
    cin  >> n;
    int coefficients[n];
    for (int i = 0; i < n; i++) {
        cin >> coefficients[i];
    }
    int constant;
    cin >> constant;
    int noOfSolutionsDP = numberOfSolutionsDP(coefficients, n, constant);
    cout << noOfSolutionsDP << "\n";
    return 0;
}

Suppose we run the program with an input as follows:

2
3 1
9

Let us see how the program runs step-by-step:

1. First of all, the variables are: coefficients[] = {3, 1}, n = 2 and constant = 9.

2. An array, table[] is defined with size 10, and all its values are initialized to 0, except the first one, which is initialized to 1.

   So, the table, in the beginning, is:

   table[0] = 1
   table[1] = 0
   table[2] = 0
   table[3] = 0
   table[4] = 0
   table[5] = 0
   table[6] = 0
   table[7] = 0
   table[8] = 0
   table[9] = 0  
   

3. i = 0

   for j = coefficients[i] (=3) to j = constant (=9):
       table[j] += table[j - coefficients[i]]
   

   So, table[3] += table[3 - 3] = table[0] = 1
   Doing similar operations for table[4] to table[9], the table gets modified to this:

   table[0] = 1
   table[1] = 0
   table[2] = 0
   table[3] = 1
   table[4] = 0
   table[5] = 0
   table[6] = 1
   table[7] = 0
   table[8] = 0
   table[9] = 1
   

4. i = 1

   for j = coefficients[i] (=1) to j = constant (=9):
       table[j] += table[j - coefficients[i]]
   

   So, table[1] += table[1 - 1] = table[0] = 1
   Doing similar operations for table[2] to table[9], the table gets modified to this:

   table[0] = 1
   table[1] = 1
   table[2] = 1
   table[3] = 2
   table[4] = 2
   table[5] = 2
   table[6] = 3
   table[7] = 3
   table[8] = 3
   table[9] = 4
   

5. Both the loops end here, and then, the value of table[constant] is table[9] = 4. Hence, the returned value is 4.

Hence, this algorithm gives an answer of 4 for the given input, which is correct since the only non-negative integral solutions for the equation 3x₁ + x₂ = 9 in the form (x₁, x₂) are: (0, 9), (1, 6), (2, 3) and (3, 0).

The time complexity of this algorithm is n * constant (n * 9 in this case, since here constant = 9), since the for loop runs twice, first time looping through n, and second time looping through the value of constant. And this complexity is much better than the exponential time we were getting from the brute force algorithm.

You can hit me up on Twitter or LinkedIn!