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.

You want to prove that a certain polynomial is positive or to decompose it as a sum of squares ? For the first problem you may compute the number of real roots (which is algorithmically simple). For the second one is there an algorithm to do it ?

Sorrry. I meant decompose the polynomial as sum of squares.