The question is not really precise at the points that matter for an algorithm. Please let me complete with some words...
So is there an algorithm that returns
the representation as sum of two (?) polynomials with coefficients (which) are (in) some
"simple" form - like (integers or)
roots of some integer coefficient
polynomial with degree not too
large?
In this situation, i'd like to be more precise, thus fixing the frame:
We give a polynomial $P\in\mathbb Z[x]$ and know it has a simple representation as a sum of two polynomials $A,B\in\mathbb Z[x]$. (Any real polynomial in one variable takes non-negative values on $\mathbb R$ if and only if it is a sum of two squares of real polynomials.)
A first question is now if this representation is unique, else we cannot ask for "the" representation. Simple examples show that we may have more representations for more reasons, e.g. $65 = 1^2+8^2 = 7^2 + 4^2$, or $65x^2 = x^2+(8x)^2 = (7x)^2 + (4x)^2$, and it is hard to make sentences, that eliminate this "trivial case", but even then one may come with the not so trivial, but in an other sense trivial example $(5x)^2 + 5^2 = 25x^2 + 25 =(3x-4)^2+(4x+3)^2. So i think it is clear that we cannot assume uniqueness.
Then the relation
$$ P = A^2 + B^2 $$
is a relation in $\mathbb Q[x]$, which transferred to the ring
$R=\mathbb Q[i][x]$ becomes:
$$ P = (A+iB)(A-iB) \ ,$$
so an algorithm would be to factor over this ring!
Sage can do this.
From this decomposition and from each choice of an element $u+iv$ of norm one, $u,v\in \mathbb Q$, $u^2 + v^2 = 1$, we can then use the "contorsion" $P=(u+iv)(A+iB)\cdot (u-iv)(A-iB)$ to get some other solution $A'+iB'=(u+iv)(A+iB)$ instead of the $A+iB$.
In the given case we can proceed as follows:
K.<j> = QuadraticField( -1 )
R.<x> = PolynomialRing( K, sparse=False )
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 ) / 58
# P.factor()
print "In our case we have exactly two simple factors, so we take f to be the first one..."
print "The the other one, g, is its Galois conjugate..."
f = P.factor()[0][0] # the first factor, taken without multiplicity
fcoeffs = f.coefficients()
gcoeffs = [ c.galois_conjugate() for c in fcoeffs ]
g = R( gcoeffs )
print "The factors are:"
print "f = %s" % f
print "g = %s" % g
print "\nWe associate A = (f+g)/2 and B = (f-g)/2/j and get"
A = (f+g)/2
B = (f-g)/2/j
print "A = %s" % A
print "B = %s" % B
print "Final check: Is P == A^2 + B^2 ? %s" % bool( P == A^2+B^2 )
I took the monic version of $P$. Results:
In our case we have exactly two simple factors, so we take f to be the first one...
The the other one, g, is its Galois conjugate...
The factors are:
f = x^5 + (-49/58*j - 21/58)*x^4 + (33/58*j - 19/58)*x^3 + (-9/29*j + 21/29)*x^2 + (-3/29*j + 7/29)*x - 17/29*j - 28/29
g = x^5 + (49/58*j - 21/58)*x^4 + (-33/58*j - 19/58)*x^3 + (9/29*j + 21/29)*x^2 + (3/29*j + 7/29)*x + 17/29*j - 28/29
We associate A = (f+g)/2 and B = (f-g)/2/j and get
A = x^5 - 21/58*x^4 - 19/58*x^3 + 21/29*x^2 + 7/29*x - 28/29
B = -49/58*x^4 + 33/58*x^3 - 9/29*x^2 - 3/29*x - 17/29
Final check: Is P == A^2 + B^2 ? True
To find the "other solution", we have to combine in an other way...
In this case we knew there is a "beautiful solution", but in general...
I will pick an other example, where the positivity is clear, but where we cannot expect to factorize over the Gaussian integers.
$$ P = x^6 + 4x^4 + 4x + 2= (x^2)^3 +(2x^2-1)^2+(2x+1)^2\ . $$
Then
sage: P = R( x^6 + 4*x^4 + 4*x + 2 )
sage: P.factor()
x^6 + 4*x^4 + 4*x + 2
So we need a bigger ring. The following still works:
R.<x> = PolynomialRing( QQbar, sparse=False )
P = R( x^6 + 4*x^4 + 4*x + 2 )
P.factor()
f1, f2, f3, f4, f5, f6 = [ f[0] for f in P.factor() ]
Then we have the polynomials:
sage: f1
x - 0.7197390712152117? - 0.8347870198157105?*I
sage: f2
x - 0.7197390712152117? + 0.8347870198157105?*I
sage: f3
x + 0.10620918972850244? - 2.047873213713062?*I
sage: f4
x + 0.10620918972850244? + 2.047873213713062?*I
sage: f5
x + 0.6135298814867092? - 0.12275817772089546?*I
sage: f6
x + 0.6135298814867092? + 0.12275817772089546?*I
see how the complex conjugation "builds pairs", consider then:
sage: f = f1*f3*f5
sage: g = f2*f4*f6
sage: f
x^3 + (0.?e-18 - 3.005418411249669?*I)*x^2 + (-2.516269913339240? - 0.3080100693808800?*I)*x - 0.9256991333675843? + 1.069149715653755?*I
sage: g
x^3 + (0.?e-18 + 3.005418411249669?*I)*x^2 + (-2.516269913339240? + 0.3080100693808800?*I)*x - 0.9256991333675843? - 1.069149715653755?*I
and obtain in a similar way two polynomials $A,B$. But:
sage: A = f+g
sage: A
already takes a looong time. Of course, understanding such coefficients is possibly / in general as hard as understanding the Galois closure of the splitting field of $P$.
Which is the question exactly?
I saw problems where was given polynomials with real coefficients and one variable and was asked to prove them to be positive. I was wondering a method to show them positive by completing them to a sum of squares. But I was not sure in what cases that gives a full proof as the coefficients are not always in closed form. I managed to do such proofs using Samuelson's inequality and then Sturm's theorem but I was wondering in which cases the sum of squares method provides a shorter proof.
None of the three questions above is answered, moreover, some new actors come into play, making the situation more unclear. Let us go back to the question, as is stated. To see why, just take a look at the <-1> i have got because somebody considered my answer to be misleading, or unrelated, or offensive, or mathematically wrong, or possibly abusively using bold face words, or all together. Not my feeling, so this is the reason why i still politely insist to make things clear. The points again, starting a new try to make things clear: