I was able to derive an algorithm to extract all solutions that satisfy the constraints in the general case. However I'm still not able to count that solutions.
For simplicity let consider the case with three variables: $$ax + by + cz = z$$ with $x,y,z \in [\theta_L, \theta_H]$.
If the g.c.d. condition it is satisfied, it is possible to derive a parametric solution, using various algorithms, see for example Morito (1979) or Fox (2000), of this type:
$$x = x_0 + x_{F_1}k_1 +x_{F_2}k_2 $$
$$y = y_0 + y_{F_1}k_1 +y_{F_2}k_2 $$
$$x = z_0 + z_{F_1}k_1 +z_{F_2}k_2 $$
where $(x_o, y_0, z_0)$ is a solution of the diophantine equation, $(x_{F_i}, y_{F_i}, z_{F_i})$ are called generator vector and $k_1, k_2$ are any integer
Now I can replace that solution in the initial equation, and I obtain a diophantine equation in two variables, that is: $$(ax_{F₁}+by_{F₁}+cz_{F₁})k₁+(ax_{F₂}+by_{F₂}+cz_{F₂})k₂=k-ax₀-by₀-cz₀$$
For simplicity just change name to the parameters: $$a_1k_1 + a_2k_2 = a_3$$ The solution of this equation will have this form: $$k_1 = a_{01} + a_{11}k_{11}$$ $$k_2 = a_{02} +a_{12}k_{11} $$
Replace in the second diophantine equation to obtain: $$a_1(a_{01} + a_{11}k_{11}) + a_2(a_{02} +a_{12}k_{11}) = a_3$$ It is possible now to solve for $k_{11}$: $$k_{11}=(a_3 -a_1a_{01} -a_2a_{02} )/(a_1a_{11}-a_2a_{12})$$ Now replace $k_{11}$ within the equation of $k_1$ and obtain: $$k_1 = a_{01}+a_{11}*((a_3 -a_1a_{01} -a_2a_{02} )/(a_1a_{11}-a_2a_{12}))$$. In theory you could also replace $k_2$ but don't do it. Instead now we have $x, y ,z$ as a function of only one variable, that is $k_2$
The following algorithm is finite and will extract all the solutions that satisfies the constraint in $x,y,z$. However I'm still not able to count the number of solutions in the general case (the one I need, even at the cost of some restrictions in the parameters) while this algorithm make it possible this number in every practical example with numbers (and not only to count but also actually to extract all solutions).
ALGORITHM
For every $i$, $i$ from $0$ to $n-1$ such that $i=1$ then $x=\theta_L$, $i=2$ then $x=\theta_L+1...$ $i=n-1$ then $x = \theta_H$
STEP 1: calculate the value of $k_2$ using $x = x_0 + x_{F_1}k_1 +x_{F_2}k_2 $
STEP 2: use $k_1$ and $k_2$ to obtain $y, z$
STEP 3: accept the solution if both $y,z \in [\theta_L, \theta_H]$, else reject.
Moreover if $$\theta_L\sum a_i = b$$ or $$\theta_H\sum a_i=b$$ I have only one solution. Out of this interval there are no solutions. I also know the maximum number of solutions that is somewhere within the interval and is less or equal than $r^n$ where $r$ is the number of possible values $x_i$ can take in $[\theta_L,\theta_H]$.
– Riccardo Jan 24 '14 at 00:28