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\$.