Multiplying Roots of a Polynomial

edit

You can also get exact results by using the Algebraic Field QQbar, which avoids the need of manually constructing a splitting field (like rburing shows in his answer).

In Sage 8.4, I'd also suggest setting QQbar's display option to radical instead of decimal (the default). Try both and see the difference; it's the simplest way I know how to explain it!

sage: QQbar.options.display_format = 'radical'
sage: R.<x> = PolynomialRing(QQbar)
sage: f = x^2*((x+1/x)^2-1)
sage: f.factor()
(x - 1/2*I*sqrt(3) - 1/2) * (x + 1/2*I*sqrt(3) - 1/2) * (x - 1/2*I*sqrt(3) + 1/2) * (x + 1/2*I*sqrt(3) + 1/2)
sage: mul([p.constant_coefficient() for p,m in f.factor()])
1

Do you want exact results or approximations? For exact results, you can define

R.<x> = PolynomialRing(QQ)
f = R(x^2*((x+1/x)^2-1))
K.<a> = f.splitting_field()

so you can do:

sage: f.change_ring(K).factor()
(x - a) * (x + a - 1) * (x - a + 1) * (x + a)
sage: f.change_ring(K).roots()
[(a, 1), (-a + 1, 1), (a - 1, 1), (-a, 1)]
sage: beta = f.leading_coefficient()
sage: myprod = beta*prod(r[0]^r[1] for r in f.change_ring(K).roots()); myprod
1

By Vieta's formulas what you get is just $(-1)^n f_0$ where $n = \deg(f)$ and $f_0$ is the constant coefficient of $f$.

sage: myprod == (-1)^f.degree()*f.constant_coefficient()
True

If you are wondering what $a$ is in the discussion above: it is a root of

sage: a.minpoly()
x^2 - x + 1

As for numerics, you can do e.g.:

sage: f.change_ring(QQbar).roots()
[(-0.500000000000000? - 0.866025403784439?*I, 1),
(-0.500000000000000? + 0.866025403784439?*I, 1),
 (0.500000000000000? - 0.866025403784439?*I, 1),
 (0.500000000000000? + 0.866025403784439?*I, 1)]
sage: beta*prod(r[0]^r[1] for r in f.change_ring(QQbar).roots())
1.000000000000000? + 0.?e-18*I

Let's see...

sage: f=x^2*((x+1/x)^2-1)

The roots are :

sage: f.roots(multiplicities=False)
[-sqrt(1/2*I*sqrt(3) - 1/2),
 sqrt(1/2*I*sqrt(3) - 1/2),
 -sqrt(-1/2*I*sqrt(3) - 1/2),
 sqrt(-1/2*I*sqrt(3) - 1/2)]

Therefore, the factorized polynom is :

sage: prod([x-u for u in f.roots(multiplicities=False)])
(x + sqrt(1/2*I*sqrt(3) - 1/2))*(x - sqrt(1/2*I*sqrt(3) - 1/2))*(x + sqrt(-1/2*I*sqrt(3) - 1/2))*(x - sqrt(-1/2*I*sqrt(3) - 1/2))

And the roots' products is:

sage: prod(f.roots(multiplicities=False)).expand()
1

Whereas moving to a polynomial ring, more specialized, is often useful, it is not necessary here...