Count the number of positive solutions for a linear diophantine equation

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$.

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

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,