## Thursday, March 16, 2023

### 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()
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))


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/

#### 1 comment:

1. Actually, for solving a diophantine equation in three variables 2.9 sec sounds rather large :-). You could just do an LLL transformation to find the lattice of integer points on the plane, and then you have a generator for all the integer feasible solutions.