Let me answer your question as best I can; but I'll start by correcting a misconception in your question.
It actually never happens that "all units are powers of just one unit" except when the unit group is finite. The reason this happens is because $-1$ gets in the way! Suppose there's a unit $u$ such that every other unit $v$ can be written as $u^n$ for some $n \in \mathbf{Z}$. But $v = -1$ is a perfectly good unit, so we must be able to write $u^n = -1$; so this equation forces $u$ to be a root of unity, and in particular it can't generate an infinite group.
When the unit group is infinite, the best we can hope for is that there is a unit $u$ such that every unit $v$ is $u^n z$ where $z$ is a root of unity. This is precisely what happens for real quadratic fields: the fundamental unit $u$ has the property that every unit is $\pm u^n$ for some $n \in \mathbf{Z}$ (and you can't get rid of the $\pm 1$ there).
When $K$ has bigger degree you need more than one unit, and Dirichlet's unit theorem gives a precise formula for how many units you need. If the number field $K$ is generated by a single algebraic number $\theta$, then you look at the minimal polynomial of $\theta$ and count how many real and non-real roots it has: it will have $r$ real roots and $s$ conjugate pairs of non-real roots, for some $r$ and $s$. Dirichlet's theorem says that the smallest collection of units you need to get within a root of unity of every unit in $K$ is of size $r + s - 1$.
For instance, I'm particularly fond of the field $K = \mathbf{Q}(\sqrt[3]{2})$. The minimal polynomial of $\sqrt[3]{2}$ is $X^3 - 2$ which has one real root and two conjugate non-real roots; so in this case one unit suffices and the field $K$ does have a fundamental unit, which turns out to be $\sqrt[3]{2} - 1$.
On the other hand, the field $L = \mathbf{Q}(\sqrt[4]{2})$ corresponds to $X^4 - 2$, which has two real roots and one pair of conjugate roots; so we'll need $2 + 1 - 1 = 2$ units to generate everything -- there's no "fundamental unit" for $L$. There are computer programs for calculating in number fields, and my computer took approximately 0.07 seconds to tell me that if we take $u_1 = \sqrt[4]{2} + 1$ and $u_2 = \sqrt{2} - 1$ then every unit of $L$ is of the form $\pm u_1^a u_2^b$ for some $a, b \in \mathbf{Z}$.
PS: I will try and answer your updated questions.
You ask whether it wouldn't be better to take $\{ \sqrt[4]{2} + 1, \sqrt[4]{2} - 1\}$ rather than $\{ \sqrt[4]{2} + 1, \sqrt{2} - 1\}$ as "fundamental units". This is entirely a matter of taste: there's no really 'best' set to take, and lots of sets will work equally well. When $r + s -1$ is 1, then a fundamental $u$ will be unique up to replacing $u$ with $\omega \cdot u^{\pm 1}$ for $\omega$ one of the finite set of roots of unity in the field (usually just $\pm 1$); but as soon as you step to larger degrees there'll be an infinite amount of choice and it is a rather fruitless exercise to try and single out any one choice which is 'best'.
As for actually computing units in practice: the bible for such computations is Henri Cohen's book "A Course in Computational Algebraic Number Theory".
Let $\sigma_1,\sigma_2$ be the real and comlex embeddings, and let $f = (\sigma_1,\sigma_2) : K \to \Bbb R \times \Bbb C$.
Then $f$ preserves addition and multiplication, and $N(x) = \sigma_1(x)|\sigma_2(x)|$.
So units are all on the surface $S = \{(y,z) \in \Bbb R \times \Bbb C \mid y|z|=1 \}$.
Suppose that you have found a unit (of norm $1$) $u$ with $\sigma_1(u) > 1$. (If you got $\sigma_1(u)<1$, just pick $1/u$ instead).
Then any $n$th root of $u$ is sent by $f$ onto the compact piece of the surface $S$ given by the condition $1 \le y \le \sigma_1(u)$.
Then you only have to look for possible lattice points on that compact surface.
More precisely this gives you some bounds on the coordinates of the possible roots of $u$ in whatever basis for $\mathfrak O_K$ you have, then it is a finite computation to check all of them.
The elements of norm $-1$ are $(-1)$ times an element of norm $1$ so you don't need to check for a possible square root of $u$ of norm $-1$ once you have found the $u$ with the $y$-component closest to $1$.
Here's how to check @Lubin's comment :
Let $u = 2^{1/3}$ and $v = u-1$. We want to be sure that $v$ is a fundamental unit.
we have $f(a+bu+cu^2) = (a+bu+cu^2, a-\frac12(bu+cu^2) + i \frac{\sqrt3}2(bu-cu^2))$
Then $f^{-1}(x,y+iz) = (\frac 13(x+2y), \frac 1{3u}(x-y+\sqrt3 z), \frac 1{3u^2}(x-y-\sqrt3 z))$
Very roughly, the surface $x|y+iz|=1$, restricted to $\sqrt v \le x \le 1$ is
bounded by $|y|,|z| \le \frac 1{|x|} \le \frac 1 {\sqrt v}$.
Going back to $a,b,c$ we get the bounds
$-1.137698 = \frac 13 (\sqrt v -2/\sqrt v) \le a \le \frac13 (1 + 2/\sqrt v) = 1.640973$
$-1.282880 = \frac 1{3u}(\sqrt v-(1+\sqrt3)/\sqrt v) \le b \le \frac 1{3u}(1+(1+\sqrt 3)/\sqrt v) = 1.682329$
$-1.018222 = \frac 1{3u^2}(\sqrt v-(1+\sqrt3)/\sqrt v) \le c \le \frac 1{3u^2}(1+(1+\sqrt 3)/\sqrt v) = 1.335266$
Hence if $v$ has a norm $1$ root $a+bu+cu^2$, we have $a,b,c \in \{-1;0;1\}$ which leaves $27$ candidates. The elements of norm $1$ among those are $1,-1+u=v,1+u+u^2=1/v$.
As a reality check we can compute the first few roots in $\Bbb R \times \Bbb C$ then look at their coefficients in $\Bbb R \otimes_\Bbb Q K$
The square roots of norm $1$ are
$(-0.101474-0.371471u+0.679930u^2),(-0.441357+0.641236u-0.465818u^2)$
The cube roots are
$(0.480750-0.480750u+0.480750u^2),
(-0.605707+0.605707u+0.302853u^2),
(0.763143+0.381571u-0.381571u^2)$
Best Answer
Here’s the story, though I leave it to others to furnish a proof or give a proper reference.
For $d\ge2$ and not a square, you get all solutions to the Pell equation $m^2-dn^2=\pm1$ by looking at the continued fraction expansion of $\sqrt d$. If $k=\lfloor\sqrt d\rfloor$, then your expansion looks like this: $$ \sqrt d=k+\bigl[\frac1{\delta_1+}\>\frac1{\delta_2+}\cdots\frac1{\delta_{r-2}+}\>\frac1{\delta_{r-1}+}\>\frac1{2k+}\>\bigr]\,, $$ where the part in brackets repeats infinitely. Then, every time you evaluate a convergent to the continued fraction just before the appearance of the $2k$, you’ll get a solution of Pell from the numerator $m$ and the denominator $n$. For instance, $\sqrt7=2+\frac1{1+}\,\frac1{1+}\,\frac1{1+}\,\frac1{4+}\,\cdots$, repeating with a period of length four. You evaluate $2+\frac1{1+}\,\frac1{1+}\,\frac1{1}=8/3$, and lo and behold, the first solution of $m^2-7n^2=\pm1$ is $m=8$, $n=3$.
For $d$ squarefree and incongruent to $1$ modulo $4$, this gives you all the units in $\Bbb Q(\sqrt{d}\,)$, but for $d\equiv1\pmod4$, it may happen that what you get is the cube of a unit. But as I recall, in that case, the primitive unit’s coordinates always will show up as the numerator and denominator of an earlier convergent.