# Count the number of positive solutions for a linear diophantine equation

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.

• 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. Apr 3 '11 at 17:18

## 4 Answers

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

• Sorry, you lost me at the last step. Why are there 6 cases and what are these 6 cases? Apr 3 '11 at 20:31
• 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. Apr 3 '11 at 20:45
• 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? Apr 4 '11 at 5:27
• Ehrhart polynomials are defined for inequalities rather than equations. Is there an easy way to use them to count solutions to equations? Nov 28 '18 at 16:02
• @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)$ Nov 29 '18 at 17:09

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 $$\sum_{k=1}^N{a_kx_k} = n, \space x\in\mathbb{N}^N$$

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{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}$$

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,

• There is no need of passing j around since it's just the length of list a
– fons
Jul 3 '16 at 1:40
• Can you please explain a bit the two-line code in the 'else' clause. Dec 24 '17 at 1:18
• p(a[1:], m - k * a[0], j - 1) is the number solutions once you fixed the first variable x[0] to some $k$ ... Dec 24 '17 at 9:27
• I feel only one answer must be correct. I mean one yielding $2282$. Jan 25 at 10:12
• it is stated in the answer that $P_3(3,2,1,1,40)=2282$ right, his code consider nonnegative solution if I understand correctly. Jan 26 at 2:35

For the number of solutions of:

$$\sum_{k=1}^{m}p_k a_k = M, p_k,a_k \in \mathbb{N}$$

Make a generating function:

$$P(x) = \prod_{k=1}^{m}\frac{x^{a_k}}{1-x^{a_k}}$$

Notice that the coefficient of the expansion of $$P(x)$$ standard notation $$[x^n]P(x)$$ are giving the number of solutions of:

$$\sum_{k=1}^{m}p_ka_k = n$$

Therefore the number of solutions of

$$\sum_{k=1}^{m}p_ka_k=M, p_k,a_k \in \mathbb{N}$$

is given by $$M$$-th derivative at $$0$$ divided by $$M!$$

$$\frac{1}{M!}\frac{\mathrm{d^M}P(x) }{\mathrm{d} x^M} (0)$$

Now you can use one of the methods for calculating higher derivative:

$$f^{(n)}=\lim_{h \to 0} \frac{1}{h^n}\sum_{k=0}^{n}(-1)^{k+n}\binom{n}{k}f(x+kh)$$ $$f^{(n)}=\frac{n!}{2\pi i}\oint\limits_\gamma \frac{f(z)}{z^{n+1}} \mathrm{d}z$$

Where $$\gamma$$ is a circle around the origin.

Notice that the second gives a great way to directly calculate the number of the solutions:

$$S(M)=\frac{1}{2\pi i}\oint\limits_\gamma \frac{P(z)}{z^{M+1}} \mathrm{d}z$$

In your case that is then:

$$S(M)=\frac{1}{2\pi i}\oint\limits_\gamma \frac{1}{z^{41}(1-z^3)(1-z^2)(1-z)^2} \mathrm{d}z$$

Replacing (we have a pole at $$1$$):

$$z=\frac{1}{2}e^{i\theta}$$

we have:

$$S(47)=\frac{2^{39}}{\pi}\int_{0}^{2\pi} \frac{1}{e^{40i\theta}(1-\frac{1}{8}e^{3i\theta})(1-\frac{1}{4}e^{2i\theta})(1-e^{i\theta})^2} \mathrm{d}\theta$$

(The integral is solved by partial fractions if we are to do it manually or numerically, or for larger values only asymptotically.)

Notice that out of all partial fractions the form $$\frac{1}{z^n}$$ or $$\frac{1}{z^{2n}+z^{n}+1}$$ does not count in (the integral being equal to $$0$$). So we are left with:

$$z=\frac{1}{2}e^{i\theta}$$

$$-\frac{306425}{144(z-1)}+\frac{10577}{72(z-1)^2}-\frac{83}{12(z-1)^3}+\frac{1}{6(z-1)^4}+\frac{1}{16(z+1)}$$

which gives the integral evaluation (counting that odd exponents are changing the sign)

$$\frac{306425}{144}+\frac{10577}{72}+\frac{83}{12}+\frac{1}{6}+\frac{1}{16}=2282$$