# Revision history [back]

### Positive polynomial as a sum of squares

I found a question to prove the polynomial $58x^{10}-42x^9+11x^8+42x^7+53x^6-160x^5+118x^4+22x^3-56x^2-20x+74$ to be positive on $\mathbb{R}$. The solution was $(7x^5-4x^3+6x^2+2x-5)^2+(-3x^5+7x^4-3x^3+7)^2$. Is there a method in Sage to find such a representations? I tried the code modified from https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/

R = AA[x]
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)
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 % 2 != 0:
return False
else:
Q = Q * (R(x)-i)^i
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())^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)
R = AA[x]
P = R(58*x^10-42*x^9+11*x^8+42*x^7+53*x^6-160*x^5+118*x^4+22*x^3-56*x^2-20*x+74)
A,B = sum_of_two_squares(P)
A
7.615773105863908?*x^5 - 2.75743509005418?*x^4 - 44.2356981337271?*x^3 + 11.9854478203666?*x^2 + 25.5052589522359?*x - 3.05770635146248?


It looks like the code does not always find the simplest form of solution. So is there an algorithm that return the representation with coefficients are some "simple" form like integers or roots of some integer coefficient polynomial with degree not too large?

### Positive polynomial as a sum of squares

I found a question to prove the polynomial $58x^{10}-42x^9+11x^8+42x^7+53x^6-160x^5+118x^4+22x^3-56x^2-20x+74$ to be positive on $\mathbb{R}$. The solution was $(7x^5-4x^3+6x^2+2x-5)^2+(-3x^5+7x^4-3x^3+7)^2$. Is there a method in Sage to find such a representations? I tried the code modified from https://ask.sagemath.org/question/10062/polynomials-as-a-sum-of-squares/

R = AA[x]
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)
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 % 2 != 0:
return False
else:
Q = Q * (R(x)-i)^i
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())^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)
R = AA[x]
P = R(58*x^10-42*x^9+11*x^8+42*x^7+53*x^6-160*x^5+118*x^4+22*x^3-56*x^2-20*x+74)
A,B = sum_of_two_squares(P)
A
7.615773105863908?*x^5 - 2.75743509005418?*x^4 - 44.2356981337271?*x^3 + 11.9854478203666?*x^2 + 25.5052589522359?*x - 3.05770635146248?


It looks like the code does not always find the simplest form of solution. So is there an algorithm that return the tries to find a representation with such that it first tries to represent the polynomial as a sum of two constants, then sum of two polynomials of degree at most 1, then sum of two polynomials of degree at most 2 and so on and in every case the one can see the minimal polynomial of the coefficients are some "simple" form like integers or roots of some integer coefficient polynomial with degree not too large?over $\mathbb Q[x]$.