**Short answer**: the vector $(s_z\,(z + s_z) - x^2, -x y, -x\,(z + s_z))$ with $s_z := \text{sign}(z) \, \|(x,y,z)\|$ is orthogonal to the vector $(x,y,z)$.

Note that we assume that $\text{sign}(x)$ is defined as $+1$ for $x \ge 0$ and as $-1$ otherwise.

Let $(x,y,z)$ be a vector with norm s and z > -s then the following matrix is an orthogonal basis where every basis vector has norm s:

$\left(
\begin{array}{ccc}
s - \frac{x^2}{z+s} & -\frac{x y}{z+s} & x \\
-\frac{x y}{z+s} & s - \frac{y^2}{z+s} & y \\
-x & -y & z \\
\end{array}
\right)$

There are two notable cases if z = -s:

- The vector is of form $(0,0,z)$ with z < 0 and we can simply invert it before applying the formula above. As shown below this can be exploited to get a branch-free implementation.
- The vector is the zero vector $(0,0,0)$. "perpendicular" doesn't make much sense in case of the null vector. If you interpret it as "dot product is zero" than you can just return the zero vector.

We can deal with these two problems as follows:

Let's look at the first vector: $(s - \frac{x^2}{z+s}, -\frac{x y}{z+s}, -x)$. The singularity at $(0,0,-1)$ can be avoided by inverting the input vector and then inverting the result which gives: $(-s - \frac{x^2}{z-s}, -\frac{x y}{z-s}, -x)$.

Following this idea we can set $s_z := \text{sign}(z) \, s$ and compute an orthogonal basis vector for *any* non-null vector $(x,y,z)$ as:

$(s_z - \frac{x^2}{z + s_z}, -\frac{x y}{z + s_z}, -x)$

This leads to a nice branch-free C++ implementation for a normalized vector:

```
Vector3 OrthoNormalVector(double x, double y, double z) {
const double g = std::copysign(1., z);
const double h = z + g;
return Vector3(g - x*x/h, -x*y/h, -x);
}
```

Check the implementation of copysign on your platform to make sure that copysign(1., 0.) returns 1 and not 0.

For an arbitrary vector, not necessarily normalized, we can use a little trick to get an orthogonal vector: we scale the vector by the factor $z+s_z$ to get:

$(s_z\,(z + s_z) - x^2, -x y, -x\,(z + s_z))$

This vector is still orthogonal to the original vector $(x,y,z)$ as it was just scaled by a factor. It also has zero norm if and only if the norm of the original vector is 0.

This leads again to a branch-free implementation:

```
Vector3 OrthogonalVector(double x, double y, double z) {
const double s = std::sqrt(x*x + y*y + z*z);
const double g = std::copysign(s, z); // note s instead of 1
const double h = z + g;
return Vector3(g*h - x*x, -x*y, -x*h);
}
```