Revision history [back]
 | 1 | initial version |
There is nothing "out of the box" within Sage but it is realtively straightforward to design an algorithm. First of all, given a polynomial you need to compute its roots which can be done as follows
sage: R.<x> = PolynomialRing(QQ,'x')
sage: P = x^6 - x^5 + x^4 - x^3 + x^2 - x + 2/5
sage: r = P.roots(CDF)
sage: print r
[(-0.54474977598 - 0.816253990699*I, 1),
(-0.54474977598 + 0.816253990699*I, 1),
(0.355140182797 - 0.844015166249*I, 1),
(0.355140182797 + 0.844015166249*I, 1),
(0.689609593183 - 0.140734077916*I, 1),
(0.689609593183 + 0.140734077916*I, 1)]
Ideally, you would use QQbar (algebraic numbers: exact computation) instead of CDF (complex double field: approximate computation) but I found that it was just impossible to do the computation within QQbar... As you see in the result above, the roots are all complex (which proves that your polynomial is positive as its dominant term is positive). Moreover, the list is well ordered by conjugates and you can rewrite your polynomial as a product of polynomials of degree 2 as follows
sage: (r0,m0),(r1,m1),(r2,m2),(r3,m3),(r4,m4),(r5,m5) = r
sage: P1 = x^2 - (r0 + r1) * x + (r0*r1)
sage: P2 = x^2 - (r2 + r3) * x + (r2*r3)
sage: P3 = x^2 - (r4 + r5) * x + (r4*r5)
sage: print P1 * P2 * P3
x^6 - x^5 + x^4 - x^3 + x^2 - x + 0.4
The above form is not that useful, but each of the polynomial P1, P2, P3 may be rewritten as a sum of two squares :
sage: a1 = (r0+r1)/2
sage: b1 = (r0*r1 - (r0+r1)**2/4).real().sqrt()
sage: a2 = (r2+r3)/2
sage: b2 = (r2*r3 - (r2+r3)**2/4).real().sqrt()
sage: a3 = (r4+r5)/2
sage: b3 = (r4*r5 - (r4+r5)**2/4).real().sqrt()
sage: ((x-a1)**2 + b1**2) * ((x-a2)**2 + b2**2) * ((x-a3)**2 + b3**2)
x^6 - x^5 + x^4 - x^3 + x^2 - x + 0.4
The form $((x-a_1)^2 + b_1^2) ((x-a_2)^2 + b_2^2) ((x-a_3)^2 + b_3^2)$ is not exactly the one you asked for but you can obtain it using the identities $(a^2 + b^2)(c^2 + d^2) = (ac - bd)^2 + (ad + bc)^2$ several times.
