This particular example can be handled with horner:
In [1]: a, b, c, d, x = symbols('a, b, c, d, x')
In [2]: y = a * x ** 3 + b * x ** 2 + c * x + d
In [3]: print(y)
a*x**3 + b*x**2 + c*x + d
In [4]: print(horner(y))
d + x*(c + x*(a*x + b))
In the multivariate case the result depends on the ordering of the symbols e.g.:
In [9]: x, y = symbols("x, y")
In [10]: p = x ** 2 * y + y ** 2 * x + x ** 2 + y ** 2
In [11]: print(horner(p, x, y))
x*(x*(y + 1) + y**2) + y**2
In [12]: print(horner(p, y, x))
x**2 + y*(x**2 + y*(x + 1))
I have heard of algorithms that attempt to find an optimal ordering
for hornerisation but I don't think any is implemented in sympy.
--
Oscar