After a lot of searching, I have finally found a reference that has what I needed. Although it doesn't appear in a 'Quaternion to Euler' or 'Quaternion to Tait-Bryan' google search, when I gave up and started looking for 'Quaternion to Axis-Angle' with the intention of going through that representation as an intermediate step, I came across the following wonderful document:
Technical Concepts
Orientation, Rotation, Velocity and Acceleration,
and the SRM
Version 2.0, 20 June 2008
Author: Paul Berner, PhD
Contributors: Ralph Toms, PhD, Kevin Trott, Farid Mamaghani, David Shen,
Craig Rollins, Edward Powell, PhD
http://www.sedris.org/wg8home/Documents/WG80485.pdf
It covers a lot of the formalisms, but most importantly, shows derivations and solutions for 3-1-3 and 3-2-1 Euler angle representation. It also seems to cover inter-conversion between pretty much every other rotation representation I'm aware of, and so I would also recommend it as a good general reference.
Oh, and the actual solution for a 3-2-1 ($z-y-x$) Tait-Bryan rotation convention from that reference:
$$
\phi = \operatorname{arctan2}\left(q_2 q_3 + q_0 q_1,\frac{1}{2}-(q_1^2 + q_2^2)\right) \\
\theta = \arcsin(-2(q_1 q_3 - q_0 q_2)) \\
\psi = \operatorname{arctan2}\left(q_1 q_2 + q_0 q_3,\frac{1}{2}-(q_2^2 + q_3^2)\right)
$$
Note that the gimbal-lock situation occurs when $2(q_1 q_3 + q_0 q_2) = \pm1$ (which gives a $\theta$ of $\pm \frac{\pi}{2} $), so it can be clearly identified before you attempt to evaluate $\phi$ and $\psi$.
(Convention for arctan2 is $\operatorname{arctan2}(y, x)$, as hinted on page 3.)
Any nonzero quaternion you like will give you a rotation; it doesn't matter how big or small its coefficients are. This is because quaternions act on vectors by conjugation: that is, the rotation of a vector $v$ by a quaternion $q$ is a the vector $qvq^{-1}$ (where we think of vectors as being imaginary quaternions). So if $r$ is any scalar (real number), then the rotation of $v$ by the quaternion $rq$ will be
$$(rq)v(rq)^{-1}=rqv\frac{q^{-1}}{r}=qvq^{-1}$$
That is, scaling $q$ is irrelevant because $q^{-1}$ will scale inversely.
So if all you want to do is give some quaternion and have it specify a rotation, then any invertible quaternion will do. In your case, however, you want to find a random rotation using quaternions. If you do this by allowing all your coefficients to vary over $[-1,1]$, you will find that your random rotation is biased. Rotations near $\pm 1\pm i\pm j\pm k$ will occur more frequently than rotations near $\pm 1$, $\pm i$, $\pm j$, or $\pm k$, because there are more possible scalings of them in your range.
If you want to generate a uniformly random set of quaternions, you can get around this by rejection sampling: allow your coefficients to vary over $[-1,1]$, but then throw out any quaternion which has norm greater than $1$ and start over. This will cancel out the bias from the previous paragraph. Probably whoever told you that the $L_2$ norm shouldn't exceed $1$ was talking about doing something like this.
If you additionally want to ensure that the Euler angles don't exceed some specified value, that probably means more rejection sampling. The set of quaternions whose Euler angles are all small doesn't have a particularly nice description, but you can generate a random quaternion, find the corresponding Euler angles via these formulae, and throw out the quaternion if they're not small enough.
(If your Euler angle threshold is really small, this becomes inefficient and you should look for an approximate solution instead, but if it's 30˚ or 45˚, the rejection sampling approach should work fine...)
Best Answer
By definition of the regular representation, the degree of the representation is the group order i.e. $8$. If you want to explicitely find the matrices (over any field $k$), consider the action of the element on each of the other elements in the group.
Concisely, this means (as reuns noted in the comment) $$\varrho(g)_{ij} = \begin{cases} 1 & g.g_j = g_i \\ 0 & ~ \text{otherwise} \end{cases} $$
(Again, as reuns noted), it suffices to explicitely compute the $8 \times 8$ matrices for the two generators $i,j$ where $Q_8 := \langle i,j \vert i^4 = 1, i^2 = j^2, jij^{-1} = i^{-1} \rangle$. Then the representing matrices for the other elements of the group can be computed via matrix multiplication.
If you are learning representation theory of finite groups and want some explicit examples, have a look at James & Liebeck.
Also, if you are aquainted with representation theory via $\mathbb{F}G$-modules, remember that the regular representation is defined as the $\mathbb{F}$-vector space with basis $g_1, \ldots, g_n$ where $\vert G \vert = n$. Then it becomes clear why the coefficients of the representing matrix are determined that way.