Revision history [back]
 | 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)
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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)
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
| | 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
I had to ascillate between AA and QQbar because of some bug appearing in Sage (to be reported soon). Computations are exact.
