In the following, we adjust $\delta$ in the range specified by $(2)$, $\sigma_B=\pm1$, and $\sigma_C=\pm1$ until we get the required rotation.
Angles $\alpha$, $\beta$, and $\gamma$
Let $\alpha$, $\beta$, and $\gamma$ be the angles from the $X$, $Y$, and $Z$ axes to the axis of rotation, $A$. These three angles are related by the Pythagorean Theorem:
$$
\cos^2(\alpha)+\cos^2(\beta)+\cos^2(\gamma)=1\tag{1}
$$
These angles, along with the angle $\delta$, will be assumed to be in the range $[0,\pi]$.
Angles $\rho_X$, $\rho_Y$, and $\rho_Z$
Consider the diagram on the sphere
$\hspace{3.2cm}$
For each value of $\delta$ so that
$$
\sin^2(\delta)\ge\max\left(\cos^2(\alpha),\cos^2(\gamma)\right)\tag{2}
$$
the circle of radius $\delta$ about $Y$ intersects both the circle of radius $\alpha$ about $X$ and the circle of radius $\gamma$ about $Z$ at one or two points.
Using the Spherical Law of Cosines, we can compute
$$
\angle AXY=2\tan^{-1}\left(\frac{\cos(\gamma)}{\sin(\alpha)+\cos(\beta)}\right)\tag{3}
$$
For each $\alpha$ and $\delta$, there are one or two solutions for $B$. To account for this, let $\sigma_B=\pm1$. Then
$$
\begin{align}
\rho_X&=\angle AXY-\sigma_B\cos^{-1}\left(\frac{\cos(\delta)}{\sin(\alpha)}\right)\tag{4}\\
\angle YBX&=\pi-\sigma_B\cos^{-1}(\cot(\alpha)\cot(\delta))\tag{5}\\
\angle XYB&=\sigma_B\cos^{-1}\left(\frac{\cos(\alpha)}{\sin(\delta)}\right)\tag{6}
\end{align}
$$
Using the Spherical Law of Cosines, we can compute
$$
\angle YZA=2\tan^{-1}\left(\frac{\cos(\alpha)}{\sin(\gamma)+\cos(\beta)}\right)\tag{7}
$$
For each $\gamma$ and $\delta$, there are one or two solutions for $C$. To account for this, let $\sigma_C=\pm1$. Then
$$
\begin{align}
\rho_Z&=\angle YZA-\sigma_C\cos^{-1}\left(\frac{\cos(\delta)}{\sin(\gamma)}\right)\tag{8}\\
\angle ZCY&=\pi-\sigma_C\cos^{-1}(\cot(\gamma)\cot(\delta))\tag{9}\\
\angle CYZ&=\sigma_C\cos^{-1}\left(\frac{\cos(\gamma)}{\sin(\delta)}\right)\tag{10}
\end{align}
$$
For consistency, define $\sigma_A=\mathrm{sgn}(\cos(\beta))$. Since $\angle XYZ=\frac\pi2$,
$$
\begin{align}
\rho_Y&=\frac\pi2-\angle XYB-\angle CYZ\tag{11}\\
\angle XAZ&=\pi-\sigma_A\cos^{-1}(\cot(\alpha)\cot(\gamma))\tag{12}
\end{align}
$$
Angle of Rotation
The rotations $\rho_X$, $\rho_Y$, and $\rho_Z$ computed above will fix the axis $A$. Accounting for parallel transport, the rotation about the axis $A$ is equal to the total geodesic curvature
$$
\rho_X\cos(\alpha)+\rho_Y\cos(\delta)+\rho_Z\cos(\gamma)\tag{13}
$$
minus the area of the light red deltoid region (since it is traversed clockwise).
Using Girard's Theorem, the areas of the purple triangles are
$$
\begin{align}
|\triangle A|&=\angle ZXA+\angle AZX+\angle XAZ-\pi\\
|\triangle B|&=\angle BXY+\angle XYB+\angle YBX-\pi\\
|\triangle C|&=\angle YZC+\angle CYZ+\angle ZCY-\pi
\end{align}\tag{14}
$$
The areas of the white sectors are
$$
\begin{align}
|\unicode{x2AA6} X|&=\rho_X\,(1-\cos(\alpha))\\
|\unicode{x2AA6} Y|&=\rho_Y\,(1-\cos(\delta))\\
|\unicode{x2AA6} Z|&=\rho_Z\,(1-\cos(\gamma))
\end{align}\tag{15}
$$
The area of the light red deltoid is
$$
\frac\pi2-|\triangle A|-|\triangle B|-|\triangle C|-|\unicode{x2AA6} X|-|\unicode{x2AA6} Y|-|\unicode{x2AA6} Z|\tag{16}
$$
Note that depending on the signs of $\sigma_B$ and $\sigma_C$, the light red deltoid may have one or two lunes attached.
Subtracting $(16)$ from $(13)$ yields a rotation of
$$
|\triangle A|+|\triangle B|+|\triangle C|+\rho_X+\rho_Y+\rho_Z-\frac\pi2\tag{17}
$$
Four applications of Girard's Theorem reduces $(17)$ to
$$
\angle XBY+\angle YCZ+\angle ZAX-2\pi\tag{18}
$$
Example:
Suppose $\alpha=\gamma=\frac\pi3$ and $\beta=\frac\pi4$. Note that $\cos^2(\alpha)+\cos^2(\beta)+\cos^2(\gamma)=1$.
By $(2)$, we can choose any $\delta$ so that $\sin^2(\delta)\ge\frac14$. Let $\delta=\frac\pi6$.
$(3)$: $\angle AXY=0.615479708670388$
$(4)$: $\rho_X=0.615479708670388$
$(5)$: $\angle YBX=\pi$
$(6)$: $\angle XYB=0$
$(7)$: $\angle YZA=0.615479708670388$
$(8)$: $\rho_Z=0.615479708670388$
$(9)$: $\angle ZCY=\pi$
$(10)$: $\angle CYZ=0$
$(11)$: $\rho_Y=\frac\pi2$
$(12)$: $\angle XAZ=1.91063323624902$
$(18)$ says that the rotation is $1.91063323624902$.
Note that the angles specified in your question do not satisfy $(1)$:
$$
\cos^2(\alpha)+\cos^2(\beta)+\cos^2(\gamma)=2.25134287511099\ne1
$$
There is no point that is $-17^\circ$, $+40^\circ$, and $-30^\circ$ from the coordinate axes.
The simplest approach to extract correctly Euler angles from a rotation matrix for any sequence of angles is using the $\mathrm{atan2}$ function. In the end, it is done in the same way (and maybe also explained why) in the text you linked. Note that compared to other inverse trigonometric functions, $\mathrm{atan2}$ has the range $(-\pi,\pi]$ (the full circle).
In your case of the $z$-$y$-$z$ rotation, if $\sin\theta\neq 0$, then
$$
\phi=\mathrm{atan2}(R_{32},R_{31}), \quad \psi=\mathrm{atan2}(R_{23},-R_{13}).
$$
We can obtain the $\theta$ angle from the last row or the last column of $R$. For example, considering the last row of $R$, we have
$$
R_{31}\cos\phi+R_{32}\sin\phi=\sin\theta(\cos^2\phi+\sin^2\phi)=\sin\theta
$$
so
$$
\theta=\mathrm{atan2}(R_{31}\cos\phi+R_{32}\sin\phi,R_{33}).
$$
For these Euler angles to be well-defined, the condition $\sin\theta\neq 0$ is required (that is, $\theta\neq k\pi$, $k\in\mathbb{Z}$). Otherwise, e.g., if $\theta=0$, the two $z$-angles are not uniquely defined.
Best Answer
One of the advantages 'Z-Y-Z' convention is that the inverse operation is defined also with this convention.
Generally for matrices we have
$(ABC)^{-1}=C^{-1}B^{-1}A^{-1}$
For rotations 'ZYZ' we have
$(\ Rot_z(\alpha)\ Rot_y(\beta)\ Rot_z(\gamma)\ )\ ^{-1}=Rot_z(-\gamma)\ Rot_y(-\beta)\ Rot_z(-\alpha) $.
So the convention is preserved, only angles are taken with minuses and in reverse order.
If we had chosen 'ZYX' the reverse order, based on the general form for the inverse of matrices product, would have been 'XYZ' and to have the operation coded in 'ZYX' order it would require more additional calculations.