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