Polynomials as a sum of squares

edit

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
    # 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.