1 | initial version |

Hi,

here is a possible solution, where the algorithm comes from those slides, page 12.

```
# the ring
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == Q:
return (sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of f
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: is_positive(P)
True
sage: A,B = sum_of_square(P)
```

I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.

2 | No.2 Revision |

Hi,

here is a possible solution, where the algorithm comes from those slides, page 12.

```
# the ring
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return
```~~(sqrt(Q),R(0))
~~(sqrt(LC) * sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of ~~f
~~T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: is_positive(P)
True
sage: A,B =
```~~sum_of_square(P)
~~sum_of_two_squares(P)

I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.

3 | No.3 Revision |

Hi,

here is a possible solution, where the algorithm comes from those slides, page 12.

`# the `~~ring
~~ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
```

I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.

4 | No.4 Revision |

Hi,

here is a possible solution, where the algorithm comes from the example given in those slides, page 12.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
```

5 | No.5 Revision |

Hi,

here is a possible solution, where the algorithm comes from the example given in those slides, page 12.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
```

6 | No.6 Revision |

Hi,

here is a possible solution, where the algorithm comes from the example given in those slides, page 12.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field
```

7 | No.7 Revision |

Hi,

here is a possible solution, where the algorithm comes from the example given in those slides, page ~~12.~~12. Actually, when a polynomial is positive, it is possible to write it as a sum of two squares only.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
#IM = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field
```

8 | No.8 Revision |

Hi,

here is a possible solution, where the algorithm comes from the example given in those slides, page 12. Actually, when a polynomial is positive, it is possible to write it as a sum of two squares only.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
```~~#IM = R(1)
~~for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field
```

9 | No.9 Revision |

Hi,

here is a possible solution, ~~where the ~~which is exact. The algorithm comes from the example given in those slides, page 12. Actually, when a polynomial is positive, it is possible to write it as a sum of two squares only.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field
```

10 | No.10 Revision |

Hi,

here is a possible solution, which is ~~exact. ~~exact (the coefficients are assumed to be algebraic real numbers). The algorithm comes from the example given in those slides, page 12. Actually, when a polynomial is positive, it is possible to write it as a sum of two squares only.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
# IM = (R(x)-f.roots()[0][0].imag())^mult
#extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)
```

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field
```

11 | No.11 Revision |

Hi,

here is a possible solution, which is exact (the coefficients are assumed to be algebraic real numbers). The algorithm comes from the example given in those slides, page 12. Actually, when a polynomial is positive, it is possible to write it as a sum of two squares only.

```
# the ring of Univariate Polynomials in x over Algebraic Real Field
R = AA[x]
def is_positive(P):
r"""
Tells whether the Polynomial P defined on the Algebraic Real Field is
positive.
"""
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
for i in P.roots():
if i[1] % 2 != 0:
return False
return True
def sum_of_two_squares(P):
r"""
P is assumed to be a polynomial defined on the Algebraic Real Field.
Returns False is P is not positive.
Returns a pair of polynomial (A,B) such that A^2 + B^2 = P otherwise
"""
# try to convert P if it is defined on a subfield of AA, say QQ.
if P.parent() != R:
P = P.change_ring(AA)
LC = P.leading_coefficient()
if LC < 0:
return False
# Q will be the part of P with real roots.
Q = R(1)
for i in P.roots():
if i[1] % 2 != 0:
return False
else:
Q = Q * (R(x)-i[0])^i[1]
if P == LC * Q:
return (sqrt(LC) * sqrt(Q),R(0))
T = R(1)
for fact,mult in R(P/Q).factor():
f = fact.change_ring(QQbar)
T = T * (R(x)-f.roots()[0][0])^mult
```~~ ~~# ~~ IM = (R(x)-f.roots()[0][0].imag())^mult
#extract ~~extract real and imaginary part of T
RE = R(0)
IM = R(0)
for i in range(T.degree()+1):
RE += T[i].real()*R(x)^i
IM += T[i].imag()*R(x)^i
SLC = sqrt(LC)
SQ = sqrt(Q)
return (SLC*SQ*RE, SLC*SQ*IM)

To use it:

```
sage: R = AA[x]
sage: P = R(x^6+x^5+x^4-x^3+x^2-x+2/5)*R((x-1)^2*(x+4)^2)
sage: P
x^10 + 7*x^9 + 8*x^8 - 18*x^7 - 12*x^6 - 4*x^5 + 177/5*x^4 - 193/5*x^3 + 202/5*x^2 - 128/5*x + 32/5
sage: is_positive(P)
True
sage: A,B = sum_of_two_squares(P)
sage: A^2+B^2
x^10 + 7.00000000000000?*x^9 + 8.00000000000000?*x^8 - 18.000000000000?*x^7 - 12.000000000000?*x^6 - 4.000000000000?*x^5 + 35.400000000000?*x^4 - 38.600000000000?*x^3 + 40.400000000000?*x^2 - 25.6000000000000?*x + 6.40000000000000?
sage: (A^2+B^2).parent()
Univariate Polynomial Ring in x over Algebraic Real Field
```

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.