Find extra arbitrary two points for a plane, given the normal and a point that lies on the plane For a plane, I have the normal $n$, and also a point $P$ that lies on the plane.
Now, how am I going to find extra arbitrary two points ($P_1$ and $P_2$) for the plane so that these three points $P$, $P_1$ and $P_2$ completely define the plane?
The solution here suggests that one assumes a certain $x$ and $y$ to substitute into the plane equation and find the remaining  $z$. But this method is only suitable for hand calculation; it breaks down for plane $z=0$. As such, it is not suitable for computer implementation.
I would need an algorithm that is robust and can handle all the cases, any idea how to construct it? 
There is a similar question here, and the answer suggests me to use  Gram-Schmidt, but I don't see how it can be applied in my case here. 

The fundamental problem here is to find two vectors $u, v$ that are orthogonal to the 
vector $n = (n_x, n_y, n_z)^T$. Morevover, the procedure should be numerically stable. 
I assume that $|n| = 1$.
Approach 1. Find the component of $n$ that has the smallest absolute value, say $n_z = \epsilon_z$.
The vector $n$ is hence "almost" in the x/y-plane (z = 0). Take the vector $w = (0, 0, 1)^T$,
with all components zero except that one on the place of the smallest component. 
Compute $u = w \times n$. $u$ is orthogonal to $n$ and is computed numerically stable 
because it is "almost" orthogonal to $n$. Finalize the procedure by computing
$v = n \times u$.
Approach 2. Often Householder orthogonalization is recommended for ill-conditioned
orthogonalization. It could work better than the Approach 1.


*

*Build that "mirror" vector $w = (n_x + 1, n_y, n_z)$ or $w = (n_x - 1, n_y, n_z)$
whose first component has bigger absolute value.

*Compute the 3x3 Householder matrix $H = I - 2 w w^T / (w^T w)$

*The first row of H will be a unit vector parallel to $n$, and the
other two rows will be unit vectors orthogonal to $n$ and to each other, i.e. these are the vectors $u^T, v^T$.


After $u, v$ are computed, we set $P_1 = P + u, P_2 = P + v$.
