12
$\begingroup$

Given a linear Diophantine equation, how can I count the number of positive solutions?

More specifically, I am interested in the number of positive solutions for the following linear Diophantine equation:

$3w + 2x + y + z = 47$

Update: I am only interested in non-zero solutions.

  • 3
    Regarding your update, as hinted at by chandok, you can replace your equation with an equivalent equation $3a+2b+c+d=40$ by setting $a=w-1$, $b=x-1$, $c=y-1$, $d=z-1$. Then every non-negative integer solution of this equation corresponds to a positive integer solution of the original equation.2011-04-03
  • 0
    Thanks for the clarification. It wasn't clear to me what he meant by "shifting the variables".2011-04-03

3 Answers 3

14

The number of solutions of (some linear expression) = n, is always something that looks like a polynomial (see http://en.wikipedia.org/wiki/Ehrhart_polynomial) You can compute them recursively or compute enough small values for it so that it determines the polynomial.

Let me do the recursive method here. First, it's simpler for me to allow $0$ as a value, I am not sure if you allow them, if you don't, then by shifting the variables, it is the same as looking a sum that gives $40$ instead of $47$, allowing things like $x=0$ and then adding $1$ to everyone.

For any $n \ge 0$, the number of solutions of $z=n$ is $P_1(n) = 1$.

For any $n \ge 0$, the number of solutions of $y+z=n$ is $P_2(n) = P_1(n) + P_1(n-1) + \ldots P_1(0) = (n+1)*1 = n+1$.

For any $n \ge 0$, the number of solutions of $2x+y+z=n$ is $P_3(n) = P_2(n) + P_2(n-2) + \ldots $. It starts to get tricky : I have to separate two cases accorging to the parity of $n$ : $$ P_3(2m) = (2m+1) + (2m-1) + \ldots + (1) = (m+1)^2 = (n^2+4n+4)/4 $$ $$ P_3(2m+1) = (2m+2) + (2m) + \ldots + (2) = (m+1)(m+2) = (n^2+4n+3)/4$$

For the last step, you have to split into 6 cases so I will only do the one you need : $$P_4(6m+4) = P_3(6m+4) + P_3(6m+1) + \ldots + P_3(1) = \Sigma_{k=0}^m P_3(6k+1) + P_3(6k+4) $$ $$ = \Sigma_{k=0}^m ((6k+1)^2 + 4(6k+1)+3 + (6k+4)^2 + 4(6k+4)+4)/4 = \Sigma_{k=0}^m 18k^2+27k+11 $$ $$ = 3m(m+1)(2m+1)+27m(m+1)/2+11(m+1) = (12m^3 + 45m^2 + 55m +22)/2 $$

So $P_4(40) = P_4(6*6+4) = 2282$.

  • 0
    Sorry, you lost me at the last step. Why are there 6 cases and what are these 6 cases?2011-04-03
  • 0
    We have 3 tiny variations according to n mod 3, which are wether the sum starts at P3(0), P3(1), or P3(2), because whenever we add 1 to w, we decrease n by 3 because of the coefficient 3 in front of w. But we also have to carry the 2 tiny variations from earlier. According to n mod 2, we start with an "even" case or an "odd" case for P3, and this will make a tiny difference somewhere. So what happens in the small details depends on n mod 6.2011-04-03
  • 0
    Ok, thanks. So I suppose that the 6 cases would be to split $n$ into $6m$, $6m + 1$, $6m + 2$, $6m + 3$, $6m + 4$, and $6m + 5$. Correct?2011-04-04
  • 0
    Ehrhart polynomials are defined for inequalities rather than equations. Is there an easy way to use them to count solutions to equations?2018-11-28
  • 1
    @AndrewMacFie well in this case, if $f(n)$ is the number of 4-uples having $3w+2x+y+z \le n$, then you get the number of solutions for the equality by computing $f(n)-f(n-1)$2018-11-29
4

Let $\mathbb{N} := \{0, 1, 2, ...\}$ be the set of natural numbers (i.e nonnegative integers). Given $N\in\mathbb{N}\backslash\{0\}$, $n\in\mathbb{N}$, and $(a_1, a_2, ..., a_N) \in\mathbb{N}^N$, let $P_{N}(a_1, a_2, ..., a_N, n)$ be the number of solutions to the linear Diophatine equation \begin{equation} \sum_{k=1}^N{a_kx_k} = n, \space x\in\mathbb{N}^N \end{equation}

Note that $P_{N}$ defines a function from $\mathbb{N}^{N+1}$ to $\mathbb{N}\backslash\{0\}$.

Then, the following recursive formula is easily established by induction on $N$ (for $N>2$, fix $x_1\in[0, n/a_1]\cap \mathbb{N}$, and solve for $(x_2$, ..., $x_N)\in\mathbb{N}^{N-1}$) \begin{equation} \begin{split} P_1(a_1, n) = 1, \text{if $n$ is a multiple of $a_1$, $0$ otherwise;}\\ P_N(a_1, a_2, ..., a_N, n) = \sum_{x_1\in [0, n/a_1]\cap \mathbb{N}}P_{N-1}(a_2, ..., a_N, n - x_1a_1)\text{, if $N>1$} \end{split} \end{equation}

The following (naive) code snippet in Python will solve any instance of the problem in finite time.

def diophantine_count(a, n):
    """Computes the number of nonnegative solutions (x) of the linear
    Diophantine equation
        a[0] * x[0] + ... a[N-1] * x[N-1] = n

    Theory: For natural numbers a[0], a[2], ..., a[N - 1], n, and j,
    let p(a, n, j) be the number of nonnegative solutions.

    Then one has:
        p(a, m, j) = sum p(a[1:], m - k * a[0], j - 1), where the sum is taken
        over 0 <= k <= floor(m // a[0])

    Examples
    --------
    >>> diophantine_count([3, 2, 1, 1], 47)
    3572
    >>> diophantine_count([3, 2, 1, 1], 40)
    2282
    """

    def p(a, m, j):
        if j == 0:
            return int(m == 0)
        else:
            return sum([p(a[1:], m - k * a[0], j - 1)
                        for k in xrange(1 + m // a[0])])

    return p(a, n, len(a))

Now that we have the hammer, let's find the nails...

OK, returning to your problem, we have $P_3(3, 2, 1, 1, 47) = 3572$.

N.B: $P_3(3, 2, 1, 1, 40) = 2282$, in accordance with mercio's brute-force computation :)

Cheers,

  • 0
    There is no need of passing j around since it's just the length of list a2016-07-03
  • 0
    Can you please explain a bit the two-line code in the 'else' clause.2017-12-24
  • 0
    p(a[1:], m - k * a[0], j - 1) is the number solutions once you fixed the first variable x[0] to some $k$ ...2017-12-24
3

Building on dohmatob's answer, here is a Pseudo Polynomial algorithm in Python, using Dynamic Programming.

Its time (and memory) complexity is O(U*N) where N is the target sum (47 in your example) and U is the number of terms/unknowns.

def diophantine_count(a, n):
    # Dynamic Programming storage.
    # dp[i][j] is the number of ways to get a sum of exactly
    # "i", only using the "j" first terms of the equation
    dp = [[0] * len(a) for _ in range(n+1)]

    # Base case.
    # There is a single way of obtaining a sum of 0
    # (with a value of zero for all the variables)
    dp[0] = [1]*len(a)

    for i in xrange(1, n+1):
        for j, c in enumerate(a):
            # dp[j][i] can be obtained adding two sub-cases:
            # (1). using one term less (j-1), and still target the same sum (i)
            # (2). using the same terms (j), but target a (i-c) sum,
            #      where c is the coefficient of the term eliminated in (1).
            #      This is because all the solutions for (i-c) are also solutions
            #      for i by adding a c to each of them.
            s1, s2 = 0, 0
            if j > 0:
                s1 = dp[i][j-1]
            if i >= c:
                s2 = dp[i-c][j]
            dp[i][j] = s1 + s2

    return dp[n][len(a)-1]

Here's an example on how to use it, adjusting the equation to obtain the count of non-zero solutions as suggested by mercio.

>>> diophantine_count([3, 2, 1, 1], 40)
2282