Where can I find the derivation of the following formula for a geodesic in 3D Euclidean space expressed as:
Calling $\vec u$ the tangent unit vector $\mathrm dr/\mathrm d s$ along the geodesic, and $\vec n$ the unit vector normal to the surface, the resulting equation of a geodesic reads
$$\frac{\mathrm d \vec u}{\mathrm d s} = -\left ( \vec u \cdot \frac{\mathrm d \vec n}{\mathrm ds} \right ) \vec n=\vec u\times \left(\frac{\mathrm d\vec n}{\mathrm ds} \times \vec n\right)$$
?
I see a dot product between the tangent of the geodesic curve $\vec u$, and the tangent to the surface $\vec u \cdot \frac{\mathrm d \vec n}{\mathrm ds},$ that is to say, the projection of the vector representing the change in the normal vector to the surface (tangent to the surface) onto the unit tangent to the curve times the unit normal must be equal and opposite to the change in the tangent to the curve $\frac{\mathrm d \vec u}{\mathrm d s}$ in the middle expression. In the RHS there are two cross products.

There's nothing too interesting going on here. For ease of typing, I will use $'$ for $d/ds$ and let $\vec v = \vec n\times\vec u$.
The first equation just follows from $0 = (\vec u\cdot\vec n)' = \vec u'\cdot\vec n + \vec u\cdot\vec n{}'$.
The second equation is just vector algebra: Write $\vec n{}' = \lambda\vec u + \mu\vec v$, where, of course, $\lambda = \vec n{}'\cdot\vec u$. Then \begin{align*} \vec u\times(\vec n'\times\vec n) &= \vec u\times\big((\lambda\vec u + \mu\vec v)\times\vec n\big) \\ &= -\lambda\vec u\times\vec v + \mu\vec u\times\vec u = -\lambda\vec n, \end{align*} as desired.