Algorithm vs. model

From [1]:

We are given a plane defined by Ax+By+Cz-D = 0 where D is significantly larger than A,B,C and GCD(A,B,C) = 1. How would I find all points (x, y, z), where x,y,z are integers and >= 0, that lie on the plane in an efficient manner?

So the poster asks for an algorithm to find $\color{darkred}x,\color{darkred}y,\color{darkred}z \in \{0,1,2,\dots\}$ such that $$\color{darkblue}A \cdot \color{darkred}x + \color{darkblue}B \cdot \color{darkred}y + \color{darkblue}C \cdot \color{darkred}z = \color{darkblue}D$$ Besides the assumptions stated in the question, I'll further assume $\color{darkblue}A,\color{darkblue}B,\color{darkblue}C,\color{darkblue}D\gt 0$.

Programmers are inclined to think about algorithms. Here I want to do something different. Instead of writing an algorithm, I think it is a good idea to try a Constraint Programming model. For this experiment, I am going to use Python Constraint [2], which is a simple, somewhat minimalistic CP solver. We need to provide a domain for the variables, so my initial approach is just to use $\color{darkred}x,\color{darkred}y,\color{darkred}z \in \{0,1,2,\dots,\color{darkblue}D\}$. The code can be:

from constraint import *
A = 3; B = 4; C = 7; D = 1000
problem = Problem()
problem.addVariables(['x','y','z'],range(0,D+1))
problem.addConstraint(ExactSumConstraint(D,[A,B,C]))
solutions = problem.getSolutions()
print(f"number of solutions: {len(solutions)}")
display(solutions)

The output looks like:

Not bad! The code for the model is just a few lines, and the solution time is quite good (2.9 seconds).

I was a bit sloppy with the domain. We can do a bit better: $$\begin{align} & \color{darkred}x \le \left\lfloor \color{darkblue}D/\color{darkblue}A\right\rfloor\\ & \color{darkred}y \le \left\lfloor \color{darkblue}D/\color{darkblue}B\right\rfloor \\ & \color{darkred}z \le \left\lfloor \color{darkblue}D/\color{darkblue}C\right\rfloor\end{align}$$ where $\left\lfloor . \right\rfloor$ is the floor function. In Python, this can look like:

problem.addVariable('x',range(0,D//A+1))
problem.addVariable('y',range(0,D//B+1))
problem.addVariable('z',range(0,D//C+1))

Interestingly, this did not make any difference in performance (also 2.9 seconds). Looks like the CP solver is outsmarting me here. (That is to be expected, but I wanted to try this anyway.)

Conclusion

Python constraint is a great little CP solver for tasks like this. 

References

1. Calculate all points of non negative integer values that lie on a plane, https://stackoverflow.com/questions/75740138/calculate-all-points-of-non-negative-integer-values-that-lie-on-a-plane
2. Python-constraint 1.4.0, https://pypi.org/project/python-constraint/