Yes, I finally understand the logic that _ is returned by the solve function.

I guess it wouldn't matter what root I take also.

    >>> gh = solve(x**3*y-t , x)
    >>> (x**9*y**3).subs(x, _[0])  # Output is t^3 
    >>> simplify(x**9 * y**3).subs(x, _[1])  # Output is t^3

Also in the first piece of code, is finding the smallest mul, is just to identify if the factor is a function of x*y?
Wouldn't just substituting x*y as t do and save the computation?