Poly should choose radical as generator if possible

[smichr]

The following was written by @oscarbenjamin in response to a comment here about trying to recognize expressions that can be written as factors of independent generators, e.g. (cos(x) - 1)*(x - 2) which has generator cos(x) and x but they appear independently in different factors. One could also imagine the cos(x) being replaced with sqrt(x):

Poly doesn't handle square roots of symbols very well:

In [53]: Poly(sqrt(x) + x)
Out[53]: Poly(x + (sqrt(x)), x, sqrt(x), domain='ZZ')

In [54]: Poly(sqrt(x) + x, sqrt(x))
Out[54]: Poly((sqrt(x)) + x, sqrt(x), domain='ZZ[x]')

There are two possible approaches to handling that in Poly:

Identify sqrt(y) as a generator and recognise y as sqrt(y)**2. This is a relatively simple approach that reduces such expressions to multivariate polynomials that have well used routines. This works even if the coefficients are EX domain (in so far as anything does).

Construct an extension field so the domain is QQ[y,t]/(t**2-y). This can be done but there is less code for implementing such a domain and it doesn't necessarily combine well with other more complicated expressions.
Your last example can be factored by manually implementing 1.:

In [49]: a, b, c, x = symbols('a, b, c, x')

In [50]: eq = expand((sqrt(x)-a)*(sqrt(x)-b)*(x-c))

In [51]: eq
Out[51]: 
                             3/2               3/2          2
-a⋅b⋅c + a⋅b⋅x + a⋅c⋅√x - a⋅x    + b⋅c⋅√x - b⋅x    - c⋅x + x 

In [52]: eq.subs(sqrt(x), t).factor().subs(t, sqrt(x))
Out[52]: (-a + √x)⋅(-b + √x)⋅(-c + x)

The extension approach is needed for more complicated examples where it isn't possible to just substitute out one symbol of interest e.g. if you have sqrt(x-1) and sqrt(x-2) etc.

[eagleoflqj]

I'm working on this with the first approach.

[smichr]

A tool that you might find useful is one which will, given a list of generators, identify those that are factors of one another: e.g.
[x, sqrt(x), exp(y), exp(2*y), exp(x+2*y), exp(x)] -> {x: [(sqrt(x), 2)], exp(x + 2*y): [(exp(x), 1), (exp(y), 2)], exp(2*y): [(exp(y), 2)]} . As a first step you might sort these by base (expr.as_base_exp()[0]) and then for a given base, by free symbols. So that list of generators could be reduced to [a**2, a, b, b**2, c*b**2, c] where {a:sqrt(x), b:exp(y), c:exp(x)}.