We want to find a perspective camera lookat
direction $\hat{d}$ so that a given direction $\hat{a}$ in world coordinates becomes $\hat{b}$ in camera coordinates. The camera rotation matrix $R$ is defined by:
World local up
vector:
$$
\hat{y} = \begin{bmatrix} 0 & 1 & 0 \end{bmatrix}^{\mkern-1.5mu\mathsf{T}}
$$
Camera lookat
vector:
$$
\hat{d} = \begin{bmatrix} d_x & d_y & d_z \end{bmatrix}^{\mkern-1.5mu\mathsf{T}}
$$
Camera right
vector:
$$
\begin{align*}
\vec{r} &= \hat{y} \times \vec{d} \quad \text{normalized to:} \\
\hat{r} &= \frac{\begin{bmatrix} d_z & 0 & -d_x \end{bmatrix}^{\mkern-1.5mu\mathsf{T}}}{\sqrt{d_x^2 + d_z^2}} \\
&= \frac{\begin{bmatrix} d_z & 0 & -d_x \end{bmatrix}^{\mkern-1.5mu\mathsf{T}}}{\sqrt{1 - d_y^2}}
\end{align*}
$$
Camera up
vector:
$$
\begin{align*}
\vec{u}
&= \hat{d} \times \hat{r} \\
&= \frac{\begin{bmatrix} -d_x d_y & d_x^2 + d_z^2 & -d_y d_z \end{bmatrix}^{\mkern-1.5mu\mathsf{T}}}{\sqrt{d_x^2 + d_z^2}}
\end{align*}
$$
So the camera rotation matrix and its inverse are:
$$
\begin{align*}
R &=
\begin{bmatrix}
\vdots & \vdots & \vdots \\
\hat{r} & \hat{u} & \hat{d} \\
\vdots & \vdots & \vdots
\end{bmatrix}
\\
R^{-1} = R^\mathsf{T} &=
\begin{bmatrix}
\dots \hat{r}^\mathsf{T} \dots \\
\dots \hat{u}^\mathsf{T} \dots \\
\dots \hat{d}^\mathsf{T} \dots
\end{bmatrix}
\end{align*}
$$
To find $\hat{d}$ we have an equation:
$$
\begin{align*}
R \cdot \hat{a} &= \hat{b} \\
\hat{a} &= R^{-1} \hat{b} \quad \text{|| easier to solve} \\
\begin{bmatrix}
a_x \\
a_y \\
a_z
\end{bmatrix}
&=
\begin{bmatrix}
\frac{d_z}{\sqrt{d_x^2 + d_z^2}} &
0 &
\frac{-d_x}{\sqrt{d_x^2 + d_z^2}} \\
\frac{-d_x d_y}{\sqrt{d_x^2 + d_z^2}} &
\sqrt{d_x^2 + d_z^2} &
\frac{-d_y d_z}{\sqrt{d_x^2 + d_z^2}} \\
d_x & d_y & d_z
\end{bmatrix}
\begin{bmatrix}
b_x \\
b_y \\
b_z
\end{bmatrix} \\
\end{align*}
$$
That's a system of 3 equations:
$$
\left{\begin{align*}
a_x &= \frac{d_z b_x - d_x b_z}{\sqrt{d_x^2 + d_z^2}} \\
a_y &= \sqrt{d_x^2 + d_z^2} \cdot b_y - \frac{d_x d_y b_x + d_y d_z b_z}{\sqrt{d_x^2 + d_z^2}} \\
a_z &= d_x b_x + d_y b_y + d_z b_z
\end{align*}\right.
$$
Let's start with:
$$
\begin{align*}
a_y &= \sqrt{d_x^2 + d_z^2} \cdot b_y - \frac{d_x d_y b_x + d_y d_z b_z}{\sqrt{d_x^2 + d_z^2}} \\
&= \sqrt{1 - d_y^2} \cdot b_y - \frac{d_x b_x + d_z b_z}{\sqrt{1 - d_y^2}} d_y \\
&= \frac{(1 - d_y^2) b_y - (d_x b_x + d_z b_z) d_y}{\sqrt{1 - d_y^2}} \\
&= \frac{b_y - (d_x b_x + d_y b_y + d_z b_z) d_y}{\sqrt{1 - d_y^2}} \\
&= \frac{(b_y - a_z d_y)}{\sqrt{(1 - d_y^2)}} \quad \text{|| substituted $a_z$} \\
0 &= (a_y (1 - d_y^2))^2 - (b_y - a_z d_y)^2 (1 - d_y^2) \\
d_y &= \pm \frac{a_z b_y \pm a_y \sqrt{ a_y^2 + a_z^2 - b_y^2 }}{a_y^2 + a_z^2}
\end{align*}
$$
We found $d_y$ so next let's investigate:
$$
\begin{align*}
a_x &= \frac{d_z b_x - d_x b_z}{\sqrt{d_x^2 + d_z^2}} \\
(d_x^2 + d_z^2) a_x^2 &= (d_z b_x - d_x b_z)^2 \\
d_x^2 a_x^2 + d_z^2 a_x^2 &= d_z^2 b_x^2 - 2 d_z b_x d_x b_z + d_x^2 b_z^2 \\
0 &= (a_x^2 - b_x^2) d_z^2 + 2 b_x d_x b_z d_z + (a_x^2-b_z^ 2)d_x^2 \\
d_z &= \frac{-2 b_x d_x b_z \pm \sqrt{4 b_x^2 d_x^2 b_z^2 - 4(a_x^2 - b_x^2)(a_x^2-b_z^ 2)d_x^2}}{2(a_x^2-b_x^2)} \\
&= d_x \frac{-b_x b_z \pm \sqrt{b_x^2 b_z^2 - (a_x^2 - b_x^2)(a_x^2-b_z^ 2)}}{a_x^2-b_x^2} \\
&= d_x \frac{-b_x b_z \pm a_x \sqrt{-a_x^2 + b_x^2 + b_z^2}}{a_x^2-b_x^2} \\
&= d_x \cdot r \\
\end{align*}
$$
We found the ratio $r$ between $d_x$ and $d_z$ which allows solving:
$$
\begin{align*}
d_x^2 + d_y^2 + d_z^2 &= 1 \\
d_x^2 + d_x^2 r^2 &= 1 - d_y^2 \\
d_x &= \pm \sqrt{\frac{1 - d_y^2}{1 + r^2}}
\end{align*}
$$
And the lookat
vector $\hat{d}$ is solved, except that there are multiple solutions due to ambiguous signs in the equations and usually only one solution makes sense.
We could compute them all and choose the one that correctly projects $\hat{a}$ to $\hat{b}$ and is closest to the camera's previous lookat
but let's try to do better.
There are 4 ambiguous signs. Let's represent each possible combination using a 4-bit number. We'll then try all combinations for different inputs and pick the best one,
but also compare the correct signs with 4 empirically chosen signs taken from the input:
$$
\left{\begin{align*}
a_x &> 0 \\
a_z &> 0 \\
b_x &> 0 \\
b_x &> a_x
\end{align*}\right.
$$