1 | initial version |

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

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.