# Factoring polynomials with symbolic expressions

asked 2017-11-24 11:17:56 -0600

This post is a wiki. Anyone with karma >750 is welcome to improve it.

This question is about writing some code that factors a desired polynomial within $\mathbb{Q}$, $\mathbb{R}$ and $\mathbb{C}$ for educational purposes.

So far, I have the following code:

@interact
def _(p = input_box(default=x^5 + x^4 - 8*x^3 + 11*x^2 - 15*x + 2, label = 'P(x) = ')):
q.<x> = PolynomialRing(QQ, 'x')
r.<x> = PolynomialRing(RR, 'x')
c.<x> = PolynomialRing(CC, 'x')
a.<x> = PolynomialRing(AA, 'x')
qp = PolynomialRing(QQ, 'x')(p)
print(factor(qp))
rp = r(p)
print(factor(rp))
cp = c(p)
print(factor(cp))
ap = a(p)
print(factor(ap))


whose output is

(x - 2) * (x^4 + 3*x^3 - 2*x^2 + 7*x - 1)
(x - 2.00000000000000) * (x - 0.147637797932293) * (x + 3.96552349222940) * (x^2 - 0.817885694297105*x + 1.70805524786870)
(x - 2.00000000000000) * (x - 0.408942847148552 - 1.24129810909174*I) * (x - 0.408942847148552 + 1.24129810909174*I) * (x - 0.147637797932293) * (x + 3.96552349222940)
(x - 2.000000000000000?) * (x - 0.1476377979322930?) * (x + 3.965523492229398?) * (x^2 - 0.8178856942971048?*x + 1.708055247868697?)


The factorization is perfect, but I would like the output to be symbolic whenever possible (ie. written with radicals instead of decimals). I know about the use of the ring 'AA' and that the numbers ending with '?' may be translated into symbolic, but I do not know the smartest way to do this.

Which would be the smartest way to achieve this goal? In the same way, ideas to improve the previous code are also welcome

edit retag close merge delete

The key is to use radical_expression for algebraic numbers. See the answers to similar questions: https://ask.sagemath.org/questions/query:radical_expression/ .

( 2017-11-24 19:14:33 -0600 )edit

Sort by » oldest newest most voted

Here are some hints:

sage: f = factor(ap)
sage: f
(x - 2.000000000000000?) * (x - 0.1476377979322930?) * (x + 3.965523492229398?) * (x^2 - 0.8178856942971048?*x + 1.708055247868697?)

sage: list(f)
[(x - 2.000000000000000?, 1),
(x - 0.1476377979322930?, 1),
(x + 3.965523492229398?, 1),
(x^2 - 0.8178856942971048?*x + 1.708055247868697?, 1)]

sage: c = list(f)[1][0].coefficients()[0]
sage: c
-0.1476377979322930?

1/12*sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) - 1/2*sqrt(-(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 71/9/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 321/2/sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 43/6) + 3/4

Symbolic Ring


Then you should be able to loop over the factorization and the coefficients (do not forget multiplicities!) and reconstruct the symbolic polynomial

Note however, that some algebraic numbers can not be expressed by radicals (and note that even if it does, sometimes Sage is not able to do it...):

sage: d = list(f)[3][0].coefficients()[0]
sage: d
1.708055247868697?
1.708055247868697?
sage: d.degree()
6

more

@tmonteil Thank you very much! This is exactly what I needed

( 2017-11-25 04:57:44 -0600 )edit

Let me please also join the party.

"Short" answer, which gives sage code for the roots, manually computed assisted by sage:

R.<x> = PolynomialRing( QQ )
f = x^5 + x^4 - 8*x^3 + 11*x^2 - 15*x + 2
f.factor()

# the factor (x-2) is simple, the other one...
g = x^4 + 3*x^3 - 2*x^2 + 7*x - 1
e,d,c,b,a = g.coefficients()

p = ( 8*a*c - 3*b^2 ) / ( 8*a^2 )
q = ( b^3 - 4*a*b*c + 8*a^2*d ) / ( 8*a^3 )

U = 649/27 + 2*sqrt(7214/27)
V = 649/27 - 2*sqrt(7214/27)

U3 =    U ^(1/3)
V3 = -(-V)^(1/3)    # since V^(1/3) is not real

S = (1/2) * sqrt( 43/12 + U3 + V3 )
signs = [ -1, 1 ]

myroots = [ -b/4/a + s1*S + s2*sqrt( - S^2 - p/2 - s1*q/4/S )
for s1 in signs
for s2 in signs ]

for myroot in myroots:
print myroot


The printed roots are somehow parallely digestible in this form. The first printed one is after manual rearrangement numerically

sage: (
....: -1/2*sqrt(  (2/3*sqrt(7214/3) + 649/27)^(1/3)
....:           - (2/3*sqrt(7214/3) - 649/27)^(1/3) + 43/12)
....: - sqrt( - 1/4*(2/3*sqrt(7214/3) + 649/27)^(1/3)
....:         + 1/4*(2/3*sqrt(7214/3) - 649/27)^(1/3)
....:         + 107/16/sqrt(  (2/3*sqrt(7214/3) + 649/27)^(1/3)
....:                       - (2/3*sqrt(7214/3) - 649/27)^(1/3) + 43/12)
....:         + 43/24)
....: - 3/4
....: ).n()
....:
-3.96552349222940

sage: g.roots( ring=QQbar, multiplicities=False )[0]
-3.965523492229398?


and it coincides with the first root, computed by sage in QQbar.

I appologize in advance for this long answer, did the work for myself, to make it once for all times clear how to figure out the four roots of g = x^4 + 3*x^3 - 2*x^2 + 7*x - 1 in a structural way.

Here, very quickly, the problem translates in finding the roots of the factor g of degree four in

sage: R.<x> = PolynomialRing( QQ )
sage: f = x^5 + x^4 - 8*x^3 + 11*x^2 - 15*x + 2
sage: f.factor()
(x - 2) * (x^4 + 3*x^3 - 2*x^2 + 7*x - 1)
sage: g = x^4 + 3*x^3 - 2*x^2 + 7*x - 1


and to give their description using only radicals. We can build

sage: K.<a> = NumberField( g )
sage: L.<b> = K.galois_closure()
sage: L.degree()
24

sage: L
Number Field in b with defining polynomial x^24 + 24*x^23 + 334*x^22 + 8432*x^21 + ...


and so on, result was truncated before coefficients of L.polynomial() get to big.

There are many examples in many answers on this site, where we compute in one of the "big" fields AA, or QQbar, or SR. Let us work here in K and L.

(I appologize again, there will be a dry exposition till the point with the direct connection to the question, which is bold face marked in the sequel. The reader in hurry should please scroll down.)

The answer applies only for (cubic, and) quartic equations. This particular polynom g is a good example to understand the Galois group of the Galois closure of the roots of the quartic. So let us use it.

Note that we have:

sage: g.discriminant().factor()
-1 * 2^5 * 3607


and also, looking into some subfields:

for M, mapM, _ in L.subfields( degree=2 ):
print "DISC = %s for %s" % ( M.discriminant().factor(), M )


This gives the one field of degree two,

DISC = -1 * 2^3 * 3607 for Number Field in b0 with defining polynomial x^2 + 15072*x + 73412352


and the same, using degree=3 as parameter above gives:

DISC = -1 * 2^3 * 3607 for Number Field in b0 with defining polynomial x^3 + 92*x^2 - 24064*x - 6486016
DISC = -1 * 2^3 * 3607 for Number Field in b1 with defining polynomial x^3 + 92*x^2 + 3896*x + 223904
DISC = -1 * 2^3 * 3607 for Number Field in b2 with defining polynomial x^3 + 92*x^2 - 18200*x + 1414048


The solution of the quartic is given for instance here: wiki/Quartic_function#General_formula_for_roots The main idea is to find the roots in the form: $$A\pm_1 B\pm_2\sqrt{C\mp_1D}$$ where the signs $\pm_1$ and $\mp_1$ do correspond, their product is $-1$, and the sign $\pm_2$ can be chosen independently. A simple computation

var( 'A,B,C,D,X' );
signs = [ -1, 1 ]
prod( [ X - ( A+s1*B + s2*sqrt( C+s1*D ) ) for s1 in signs for s2 in signs ] ).expand()


shows how to search for A,B,C,D, when the equation is given. We usually clear first the factor in $x^3$ by a linear substitution, so that after this step $A=0$. The value of $A$ is simple. The main point is to find the value of $B$, which becomes an $S$ in the sequel.

More precisely now. The wiki article gives for the equation $$ax^4 + bx^3 + cx^2 + dx + e = 0$$ the solution $$-\frac b{4a}\pm_1 S\pm_2 \sqrt{ -S^2 -\frac p2 \mp_1 \frac q{4S}} \ ,$$ and the main complication is getting $S$. Above, $p,q$ are simple rational expressions involving $a,b,c,d$. If we give weights $0,1,2,3,4$ to the coefficients $a,b,c,d,e$, then $S$ is of weight $1$, as $b/(4a)$, $p$ of weight $2$, as $-S^2$, and $q$ of weight three.

So the particular case extracted from the question, is asking for $S$, if we have $S$, we have basicly solved and understood the situation.

Let us ask sage for the algebraic equation satisfied by the sum of two roots $x_0+x_1$, for some choice of the order of the roots. Note that considering all choices...

g_roots = g.roots( ring=L, multiplicities=False )
for k in [0..3]:
for kk in [k+1..3]:
xk  = g_roots[ k]
xkk = g_roots[kk]
print "x%s + x%s has minpoly = %s" % ( k, kk, (xk+xkk).minpoly() )


is given:

x0 + x1 has minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
x0 + x2 has minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
x0 + x3 has minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
x1 + x2 has minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
x1 + x3 has minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
x2 + x3 has minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82


and the choice of the two roots is rather immaterial. Alternatively, if we work in L:

sage: K.<a> = NumberField( g )
sage: L.<b> = K.galois_closure()
sage: maps = K.embeddings(L)
sage: len( maps )
4

sage: for m1, m2 in Combinations( maps, 2 ):
....:     print "minpoly =", ( m1(a) + m2(a) ).minpoly()
....:
minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82
minpoly = x^6 + 9*x^5 + 23*x^4 + 3*x^3 - 7*x^2 + 87*x - 82


Now we can try:

x0, x1, x2, x3 = g_roots
( x0 + x1 + 9/6 ).minpoly()    # this is (2S).minpoly()


(Above, $9/6=3/2=2\cdot b/(4a)$ is clearing the monomial in $x^5$, and it corresponds to the term $-b/(4a)$ from the wiki formula.) This gives:

x^6 - 43/4*x^4 + 995/16*x^2 - 11449/64


Which is a cubic equation in $x^2$. We can further ask for...

sage: ( ( x0 + x1 + 9/6 )^2 - 43/12 ).minpoly()
x^3 + 71/3*x - 1298/27
sage: _.discriminant().factor()
-1 * 2^5 * 3607

sage: ( ( x0 + x1 + 9/6 )^2 - 43/12 ).minpoly().roots( ring=QQbar, multiplicities=False )
[1.789260758493839?,
-0.8946303792469193? - 5.105659331867054?*I,
-0.8946303792469193? + 5.105659331867054?*I]

sage: for r in _:
....:     print r, '=', r.radical_expression()

1.789260758493839? = (2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 71/9/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)


and two further solutions. All three are of the form $$\epsilon \left( \frac{649}{27} +\frac 29\sqrt{7214\cdot 3} \right)^{1/3} + \epsilon^2\left( \frac{649}{27} -\frac 29\sqrt{7214\cdot 3} \right)^{1/3}\ ,$$ where $\epsilon$ is one of the three roots of unity, either $1$, or one of the two primitive roots, and the two radicals are taken so that their products is the real number $71/9$:

( ( 649/27 + 2/9*sqrt(7214)*sqrt(3) ) * \
( 649/27 - 2/9*sqrt(7214)*sqrt(3) )   ).expand()^(1/3)


which gives:

71/9*(-1)^(1/3)


We can now write simple sage code, using radicals, which generates the four roots of g.

coeffs = g.coefficients()
coeffs . reverse()
a,b,c,d,e = coeffs

p = ( 8*a*c - 3*b^2 ) / ( 8*a^2 )
q = ( b^3 - 4*a*b*c + 8*a^2*d ) / ( 8*a^3 )

U = 649/27 + 2*sqrt(7214/27)
V = 649/27 - 2*sqrt(7214/27)

U3 =    U ^(1/3)
V3 = -(-V)^(1/3)    # since V^(1/3) is not real

S = (1/2) * sqrt( 43/12 + U3 + V3 )
signs = [ -1, 1 ]

myroots = [ -b/4/a + s1*S + s2*sqrt( - S^2 - p/2 - s1*q/4/S )
for s1 in signs
for s2 in signs ]

print "The roots computed with bare hands starting from the sage value\nS = %s" % S
print "and setting root = %s (+/-) S + [+/-] sqrt( -S^2 - %s (-/+) %s/S )" % ( -b/4/a, -p/2, q/4 )
for r in myroots:
print r.n()

print
print "The roots computed with sage in QQbar"
for rr in g.roots( ring=QQbar, multiplicities=False ):
print rr


(The code in the above block is basicly the short answer.) Results:

The roots computed with bare hands starting from the sage value
S = 1/2*sqrt((2/3*sqrt(7214/3) + 649/27)^(1/3) - (2/3*sqrt(7214/3) - 649/27)^(1/3) + 43/12)
and setting root = -3/4 (+/-) S + [+/-] sqrt( -S^2 - 43/16 (-/+) 107/32/S )
-3.96552349222940
0.147637797932293
0.408942847148552 - 1.24129810909174*I
0.408942847148552 + 1.24129810909174*I

The roots computed with sage in QQbar
-3.965523492229398?
0.1476377979322930?
0.4089428471485524? - 1.241298109091741?*I
0.4089428471485524? + 1.241298109091741?*I


The explicit computation shows also how the Galois group of L acts on the roots (and on "similar" expressions).

sage: G.<v> = K.galois_group()
sage: G.order()
24
sage: G.gens()
[(1,12)(2,19)(3,17)(4,22)(5,11)(6,8)(7,16)(9,10)(13,14)(15,23)(18,21)(20,24),
(1,18)(2,21)(3,11)(4,16)(5,20)(6,7)(8,23)(9,17)(10,24)(12,14)(13,19)(15,22),
(1,19)(2,12)(3,24)(4,8)(5,9)(6,22)(7,15)(10,11)(13,18)(14,21)(16,23)(17,20),
(1,20,7)(2,4,11)(3,14,8)(5,22,19)(6,13,17)(9,15,18)(10,21,23)(12,16,24)]

sage: G.is_solvable()
True
sage: G.is_polycyclic()
True
sage: G.is_supersolvable()
False

sage: P = PermutationGroup( G.gens() )
sage: for N in P.normal_subgroups():
....:     print "Normal subgroup of order %s" % N.order()
....:
Normal subgroup of order 1
Normal subgroup of order 4
Normal subgroup of order 12
Normal subgroup of order 24

sage: N = [ N for N in P.normal_subgroups() if N.order() == 4 ][0]
sage: PQ = P.quotient(N)
sage: N.gens()
[(1,18)(2,21)(3,11)(4,16)(5,20)(6,7)(8,23)(9,17)(10,24)(12,14)(13,19)(15,22),
(1,19)(2,12)(3,24)(4,8)(5,9)(6,22)(7,15)(10,11)(13,18)(14,21)(16,23)(17,20)]
sage: PQ.order()
6
sage: PQ.is_commutative()
False


So $G=Gal(L:\mathbb Q)$ has a normal subgroup $\cong (Z/2)^2$ and it is a twisted product of it with the symmetric group with $3!$ elements.

The one element of order 3 acts by replacing U3, V3 above by an other choice for the roots of order three of U, V, constrained to have a fixed product.

The elements of order two act by changing signs in the places where we used $\pm$...

more

It turns out that the "brute force" approach (i.e. asking for the roots in the (default) SR ring and building the product of the resulting monomials) gives (a not especially appetizing) radical exression :

sage: prod([(x-t[0])^t[1] for t in (x^5 + x^4 - 8*x^3 + 11*x^2 - 15*x + 2).roots()])
1/20736*(12*x + sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 6*sqrt(-(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 71/9/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 321/2/sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 43/6) + 9)*(12*x + sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) - 6*sqrt(-(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 71/9/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 321/2/sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 43/6) + 9)*(12*x - sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 6*sqrt(-(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 71/9/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 321/2/sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 43/6) + 9)*(12*x - sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) - 6*sqrt(-(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) + 71/9/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 321/2/sqrt((36*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(2/3) + 129*(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3) - 284)/(2/9*sqrt(7214)*sqrt(3) + 649/27)^(1/3)) + 43/6) + 9)*(x - 2)


viewing this monster may help you to convince yourself that is indeed the product of 5 first-degree polynomials.

But note that this is pure happenstance : your polynomial being of degree 5, its pure luck that its roots can be expressed by radicals.

more