I seem to be giving lots of answers that depend on very special properties of the supplied example. Here’s an argument, tailored to your polynomial $f(x)=x^3-x-1$. Not the general method you were hoping for at all.
First set $\alpha$ to ba a root of your polynomial, which we all know is irreducible over $\Bbb Q$. It’s not too hard to calculate the discriminant of the ring $\Bbb Z[\alpha]$ as the norm down to $\Bbb Q$ of $f'(\alpha)=3\alpha^2-1$; this number is $23$, surprisingly small for a cubic extension. The fact that it’s square-free implies that $\Bbb Z[\alpha]$ is the ring of integers in the field $k=\Bbb Q(\alpha)$.
Our field $k$ clearly is not totally real, since $f$ has only one real root. So in the jargon of algebraic number theory, $r_1=r_2=1$, one real and one (pair of) complex embedding(s). We can apply the Minkowski Bound
$$
M_k=\sqrt{|\Delta_k|}\left(\frac4\pi\right)^{r_2}\frac{n!}{n^n}\,,
$$
which for $n=3$ gives a bound less than $2$, so that $\Bbb Z[\alpha]$ is automatically a principal ideal domain.
Let’s factor the number $23$ there: we certainly know that it’s not prime, since $23$ is ramified.
Now, we already know a number of norm $23$, necessarily a prime divisor of the integer $23$, it’s $3\alpha^2-1$, and indeed $23/(3\alpha^2-1)=4 + 9\alpha - 6\alpha^2$. But better than that, $23/(3\alpha^2-1)^2=3\alpha^2-4$. This number has norm $23$ (because the norm of $23$ itself is $23^3$). So we’ve found the complete factorization of $23$.
Now let’s look more closely at $f(x)=(x-\alpha)g(x)$, for a polynomial $g$ that we can discover by Euclidean Division to be $g(x)=x^2+\alpha x+\alpha^2-1$. And the roots of $g$ are the other roots of $f$; the Quadratic Formula tells you what they are, and the discriminant of $g$ is $\alpha^2-4(\alpha^2-1)=4-3\alpha^2$. which we already know as $-23/(3\alpha^2-1)^2$. Going back to the Quadratic Formula, our other roots are
$$
\rho,\rho'=\frac{\alpha\pm\sqrt\delta}2\>,\>\delta=\frac{-23}{(3\alpha^2-1)^2}\>,\>\sqrt\delta=\frac{\sqrt{-23}}{3\alpha^2-1}\>.
$$
And that gives you your roots of this one very special cubic polynomial in terms of one root $\alpha$ and $\sqrt{-23}$.
Assume $\mbox{char}(F)=2$ and let $\alpha_1,\dots,\alpha_n$ be the distinct roots of $f(X)$. Define
$$
\delta_f = \sum_{i<j} \frac{\alpha_i}{\alpha_i + \alpha_j}
$$
and
$$
D(f) = \sum_{i<j} \frac{\alpha_i\alpha_j}{\alpha_i^2 + \alpha_j^2}.
$$
The element $D(f)$ is called the Berlekamp discriminant of $f(X)$. It is straighforward to verify that
$$
\delta_f^2 + \delta_f + D(f) = 0.
$$
If $\tau\in G_f$ is the transposition $(k \; \; k+1)$, we obtain
$$
\sigma\delta_f = \delta_f + \frac{\alpha_{k}}{\alpha_k + \alpha_{k+1}} + \frac{\alpha_{k+1}}{\alpha_{k+1} + \alpha_{k}} = \delta_f + 1.
$$
Thus $\sigma\delta_f = \delta_f$ if and only if $\sigma\in A_n\cap G_f=SG_f$. By the Galois correspondence you obtain that
$$
F(\alpha_1,\dots,\alpha_n)^{SG_f} = F(\delta_f)
$$
as desired (it is necessary to use the fact that $[S_n:A_n]=2$).
Best Answer
The answer to your question just lies in a discriminant computation. Given an odd prime $p$, you know that the cyclomic field $K=\mathbf Q(\zeta_p)$ is a normal extension of $\mathbf Q$, with cyclic Galois $G\cong \mathbf F_p^*$. By cyclicity, $G$ contains a unique subgroup $H$ of index $2$, and the fixed field $F$ fixed by $H$ is the unique quadratic subfield of $K$, which you want to describe explicitly.
By definition, the dicriminant $D=D(1,\zeta_p,..., {\zeta_p}^{p-2})$ is a square in $K$. The classical formula $D=\pm N(f'(\zeta_p))$, where $f$ is the $p$-th cyclotomic polynomial, shows that $D=\pm p^{p-2}$. It follows that $F=\mathbf Q(\sqrt {\pm p})$. Actually, taking more care of signs, one gets $F=\mathbf Q(\sqrt {p^*})$, where $p^*=(-1)^{\frac {p-1}2}p$ . Short proof: introduce the subfield $K^+$ of $K$ fixed by complex conjugation, which is also the maximal totally real subfield of $K$, and equally the unique (cyclic) subextension of degree $\frac {p-1}2$ over $\mathbf Q$, so $F \subset K^+$ iff $p\equiv 1$ (mod $4$).