Yes, you can do this without any trigonometric functions, and even without any square roots. Take a look at Rodrigues' formula for a rotation matrix given an angle $\theta$ and an axis along a unit vector $k$:
$R=I\cos\theta+\sin\theta[k]_\times+(1-\cos\theta)kk^{\text T}\;,$
where $I$ is the identity matrix and
$[k]_\times=\begin{pmatrix}0&-k_z&k_y\\k_z&0&-k_x\\-k_y&k_x&0\end{pmatrix}\;.$
You have not an axis and an angle, but a vector $N$ that you want to rotate the $z$ axis into. The axis should be perpendicular to this vector and the $z$ axis, and the angle should be the angle between this vector and the $z$ axis. Since the $z$ component of your vector is the cosine of the angle it forms with the $z$ axis, you already have $\cos\theta=N_z$. Now consider
$N_\perp=\begin{pmatrix}0\\0\\1\end{pmatrix}\times N=\begin{pmatrix}-N_y\\N_x\\0\end{pmatrix}\;.$
This vector is perpendicular to both $e_z$ and $N$, so its direction is along the desired rotation axis. Its magnitude is $\sin\theta$. So
$ \begin{eqnarray} N_\perp&=&\sin\theta k\;,\\ [N_\perp]_\times&=&\sin\theta[k]_\times\;,\\ \end{eqnarray} $
where $k$ is a unit vector along the rotation axis. Now we have most of the ingredients of Rodrigues' formula. The only problem left is that we don't have $kk^{\text T}$ but $\sin^2\theta kk^{\text T}$. But $\sin^2\theta=1-\cos^2\theta$, and we have $\cos\theta$, so we can correct for that using only elementary arithmetic. Putting it all together, we have
$ \begin{eqnarray} R &=& I\cos\theta+[N_\perp]_\times+(1-\cos\theta)N_\perp N_\perp^{\text T}/(1-\cos^2\theta) \\ &=& I\cos\theta+[N_\perp]_\times+N_\perp N_\perp^{\text T}/(1+\cos\theta) \\ &=& IN_z+ \begin{pmatrix}0&0&N_x\\0&0&N_y\\-N_x&-N_y&0\end{pmatrix} +\frac1{1+N_z} \begin{pmatrix}N_y^2&-N_xN_y&0\\-N_xN_y&N_x^2&0\\0&0&0\end{pmatrix} \\ &=& \begin{pmatrix}N_z&0&N_x\\0&N_z&N_y\\-N_x&-N_y&N_z\end{pmatrix} +\frac1{1+N_z} \begin{pmatrix}N_y^2&-N_xN_y&0\\-N_xN_y&N_x^2&0\\0&0&0\end{pmatrix} \;. \end{eqnarray} $
You can check that this is an orthogonal matrix that rotates a unit vector along the $z$ axis into $N$. If you organize the operations efficiently, you only need two negations, three additions, one division and five multiplications to calculate the rotation matrix.
Note that the result becomes undefined at $N_z=-1$, corresponding to the fact that there's no preferred choice of axis in this case. You can either arbitrarily choose an axis specifically for that case, e.g. the $x$ axis, or, since this might also cause numerical problems near $N_z=-1$, you can instead solve the problem for $-N_z$ whenever $N_z<0$, and then invert the resulting vectors. This will produce a mirror image of your distribution, but I'd assume that this is rotationally invariant around the $z$ axis anyway, and in that case the inversion wouldn't make a difference.