Revision history [back]

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