If there are only unnested radicals and not more than 4 we will always be able to "do it right" because they will collect properly and not introduce further radicals after squaring. But if (as in the inline comment) we have an expression which started well but now has more than 4 radicals it could be that we made the wrong choice in our first step so we abandon the process. In fact, if we do smart inspection of the radicals (something that smells a bit ellusively recursive) we might find that if we select them in a certain way we can get more than 4 radicals to reduce to less as a result of fortuitous cancellation of terms:

```
the original expression

>>> eq=sqrt(5)*r+2*sqrt(5)*sqrt(-sqrt(5)+5)+5*r
>>> eq
sqrt(5)*sqrt(sqrt(5) + 5) + 2*sqrt(5)*sqrt(-sqrt(5) + 5) + 5*sqrt(sqrt(5) + 5)
>>> powsimp(_)
sqrt(5*sqrt(5) + 25) + 2*sqrt(-5*sqrt(5) + 25) + 5*sqrt(sqrt(5) + 5)
>>> _.args
(sqrt(5*sqrt(5) + 25), 5*sqrt(sqrt(5) + 5), 2*sqrt(-5*sqrt(5) + 25))
>>> a,b,c=_

multiply by (c - (a + b))

>>> factor_terms(powsimp(powsimp(expand(eq*(c-a-b)))), radical=True)
5*(-30 - 10*sqrt(5) - 5*sqrt(2)*sqrt(sqrt(5) + 3) - sqrt(10)*sqrt(sqrt(5) + 3) + 4*sqrt(10)*sqrt(-sqrt(5) + 3))
>>> _/5
-30 - 10*sqrt(5) - 5*sqrt(2)*sqrt(sqrt(5) + 3) - sqrt(10)*sqrt(sqrt(5) + 3) + 4*sqrt(10)*sqrt(-sqrt(5) + 3)
>>> args=_.args

run through all combinations, multiplying 
>>> X = Add(*args)
>>> for c in combinations(args,3):
...  print powsimp(powsimp(expand(X*(Add(*c) - (X - Add(*c))))))

...
-240*sqrt(-10*sqrt(5) + 30) - 200*sqrt(5) - 40*sqrt(10*sqrt(5) + 30) + 600
100*sqrt(2*sqrt(5) + 6) + 60*sqrt(10*sqrt(5) + 30) + 800 + 880*sqrt(5)
200 + 60*sqrt(10*sqrt(5) + 30) + 400*sqrt(-2*sqrt(5) + 6) + 280*sqrt(5) + 300*sqrt(2*sqrt(5) + 6)
-240*sqrt(-10*sqrt(5) + 30) - 400*sqrt(-2*sqrt(5) + 6) + 320*sqrt(5) + 1600
-240*sqrt(-10*sqrt(5) + 30) - 280*sqrt(5) + 200*sqrt(2*sqrt(5) + 6) + 1000
100*sqrt(10*sqrt(5) + 30) + 300*sqrt(2*sqrt(5) + 6) + 1200 + 800*sqrt(5)
-200*sqrt(2*sqrt(5) + 6) - 400*sqrt(-2*sqrt(5) + 6) - 200*sqrt(5) - 200
-800*sqrt(5) - 800  # <=========================================== the best
-600 + 100*sqrt(2*sqrt(5) + 6) + 280*sqrt(5) + 240*sqrt(-10*sqrt(5) + 30) + 100*
sqrt(10*sqrt(5) + 30)
-280*sqrt(5) - 400*sqrt(-2*sqrt(5) + 6) + 200 + 40*sqrt(10*sqrt(5) + 30)
```
So if we try all combinations it *might* work...but it might fail. We only know that things will get better if we have no more than 4 terms. What I am going to do is use minpoly as a fallback so that if we exceed 4 terms we see if the solution to the minpoly gives an easier expression to deal with. For example, for this equation we get the following equivalent expression:

```
>>> eq
sqrt(5)*sqrt(sqrt(5) + 5) + 2*sqrt(5)*sqrt(-sqrt(5) + 5) + 5*sqrt(sqrt(5) + 5)
>>> [w for w in solve(minpoly(eq,x)) if test_numerically(eq,w)]
[10*sqrt(sqrt(5) + 5)]

easy to solve

>>> radsimp(1/_[0])
(-sqrt(5)*sqrt(sqrt(5) + 5) + 5*sqrt(sqrt(5) + 5))/200
```