Radical expressions for roots of polynomials using Galois theory

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

Given a polynomial from ℚ[X], we need a better way to express its roots using radical expressions if such an expression is possible.

The current approach, as used by e.g. NumberFieldElement._symbolic_ (and after #14239 gets merged probably AlgebraicNumber_base.radical_expression instead), delegates this task to the solve method for expressions from the symbolic ring. That in turn will delegate to Maxima. But Maxima is not able to find a radical expression in all cases where they do exist:

sage: p = x^6-300*x^5+30361*x^4-1061610*x^3+1141893*x^2-915320*x+101724
sage: p.solve(x, explicit_solutions=True)
[]
sage: r = 1/8*((sqrt(4*(1/9*sqrt(109)*sqrt(3) + 2)^(1/3) - 4/3/(1/9*sqrt(109)*sqrt(3) + 2)^(1/3) + 17) + 5)^2 + 4)*(sqrt(4*(1/9*sqrt(109)*sqrt(3) + 2)^(1/3) - 4/3/(1/9*sqrt(109)*sqrt(3) + 2)^(1/3) + 17) + 5)
sage: r.minpoly() == QQ[x](p)
True

In #14239 comment:77 Jeroen Demeyer stated that a proper solution here would use Galois Theory, and that we might be able to leverage PARI for this. So the goal of this ticket here is a function or method which constructs radical expressions for the roots of all polynomials where PARI can compute the Galois group of the splitting field (perhaps with a bound on the degree just to limit the running time of the algorithm).

Depends on #14239

CC: @pjbruin @sagetrac-tmonteil @videlec @seblabbe @slel

Component: number fields

Keywords: radical, galois, symbolic

Issue created by migration from https://trac.sagemath.org/ticket/17516

[jdemeyer]

Description changed:

--- 
+++ 
@@ -11,4 +11,4 @@
 True

-In #14239 comment:77 Jeroen Demeyer stated that a proper solution here would use Galois Theory, and that we might be able to leverage PARI for this. So the goal of this ticket here is a function or method which constructs radical expressions for the roots of all polynomials where doing so is possible, perhaps with some explicitely stated bound on the degree. +In #14239 comment:77 Jeroen Demeyer stated that a proper solution here would use Galois Theory, and that we might be able to leverage PARI for this. So the goal of this ticket here is a function or method which constructs radical expressions for the roots of all polynomials where PARI can compute the Galois group of the splitting field (perhaps with a bound on the degree just to limit the running time of the algorithm).

[jdemeyer]

Dependencies: #14239

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

comment:3

Computing the Galois group for a polynomial works pretty fast, but I'm not sure how much use that really is. I fear we might need the map between the original polynomial and the Galois closure. Converting the Galois closure for the example above takes some time, but not so much as converting that closure to pari. The time is apparently spent somewhere inside _pari_integral_basis. Now I wonder, do we actually have to call galoisinit on the galois closure as a number field? Or could we do that call on its defining polynomial instead? Something along these lines:

sage: Qx.<x> = QQ[]
sage: p = x^6-300*x^5+30361*x^4-1061610*x^3+1141893*x^2-915320*x+101724
sage: K = NumberField(p, names="a")
sage: GC, GCm = K.galois_closure(names="b", map=True)
sage: q = GC.defining_polynomial()
sage: gal = pari(q).galoisinit()
sage: G = PermutationGroup(sorted(gal[6], cmp=cmp))
sage: G.is_solvable()
True
sage: ds = G.derived_series()

I still have to figure out how to turn that derived series into a radical expression. I hope you don't mind me posting thoughts along the way.

[jdemeyer]

comment:4

Replying to @gagern:

Computing the Galois group for a polynomial works pretty fast, but I'm not sure how much use that really is. I fear we might need the map between the original polynomial and the Galois closure. Converting the Galois closure for the example above takes some time, but not so much as converting that closure to pari. The time is apparently spent somewhere inside _pari_integral_basis. Now I wonder, do we actually have to call galoisinit on the galois closure as a number field? Or could we do that call on its defining polynomial instead? Something along these lines:

sage: Qx.<x> = QQ[]
sage: p = x^6-300*x^5+30361*x^4-1061610*x^3+1141893*x^2-915320*x+101724
sage: K = NumberField(p, names="a")
sage: GC, GCm = K.galois_closure(names="b", map=True)
sage: q = GC.defining_polynomial()
sage: gal = pari(q).galoisinit()

Yes indeed

sage: G = PermutationGroup(sorted(gal[6], cmp=cmp))
sage: G.is_solvable()
True
sage: ds = G.derived_series()

I don't think need you need to convert the group to Sage, I would use galoisfixfield() from PARI. The tricky part will be adding the roots of unity, you need to add them manually.

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

comment:5

Replying to @jdemeyer:

I don't think need you need to convert the group to Sage, I would use galoisfixfield() from PARI.

Passing which subgroup as an argument? I guess I should iterate over all normal subgroups, using galoissubgroups and galoisisnormal, which aren't available in Sage at the moment. #17519 should fix that.

I just started reading the chapter “How to solve a solvable equation” in the book Classical Galois Theory by Lisl Gaal. That sounds very promising, since it apparently concentrates on the algorithm, not the theorem.

[jdemeyer]

comment:6

Replying to @gagern:

Replying to @jdemeyer:

I don't think need you need to convert the group to Sage, I would use galoisfixfield() from PARI.

Passing which subgroup as an argument?

I think it suffices to take the information from galoisinit().gen and galoisinit().orders, but I haven't checked the details.

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

comment:7

Replying to @jdemeyer:

I think it suffices to take the information from galoisinit().gen and galoisinit().orders, but I haven't checked the details.

Current experiments seem to indicate that gal.galoisfixfield(gal[6][:-1]) might be enough. At least for the given example. But I'm not sure whether that's a general rule. What one can do is look for a sequence of normal subgroups in such a way that the index of each one in its supergroup is a prime number. In the given example, there is exactly one such sequence which ends up at the trivial group. Not sure whether that's a general rule either, and if not, whether it would be worthwhile to try different sequences. If I understand the GAP code in local/gap/latest/lib/grpperm.gi correctly, its DerivedSubgroup method makes that choice greedily, without a graph search. So if it really is the derived series we want, then we might have a closer look at that as well, or simply call GAP as I did in my first attempt.

Here is what I currently use in my experiments, to be on the safe side:

def findRadicalSeries(t, collect=None):
    g = t[0]
    n = t[1]
    for s in g.galoissubgroups():
        if g.galoisisnormal(s):
            m = ZZ(prod(s[1]))
            k = n // m
            if is_prime(k):
                t2 = (s, m, t)
                if m != 1:
                    r = findRadicalSeries(t2, collect)
                    if r is not None and collect is None:
                        return r
                else:
                    r = []
                    while t2 is not None:
                        r.append(t2[:2])
                        t2 = t2[2]
                if collect is None:
                    return r
                else:
                    collect.append(r)
    return collect
rs = findRadicalSeries((gal, ZZ(prod(gal[7])), None))

I'm currently experimenting with the sequence that found. And just caused yet another crash, this time somewhere in the symbolic expressions engine. Investigating…

[edd8e884-f507-429a-b577-5d554626c0fe]

comment:8

Just cc-ing myself, this feature was in my plans for years.

[jdemeyer]

comment:9

Replying to @gagern:

Current experiments seem to indicate that gal.galoisfixfield(gal[6][:-1]) might be enough.

If you apply this recursively (going from the given field to a subfield and then recursing), sure.

[jdemeyer]

comment:10

I doubt we need to use GAP for this. I think the subfield sequence given by initial segments of galoisinit().gen is sufficient.

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

comment:11

I still have trouble figuring out what element we should adjoin. One option (Lagrange resolvent?) would apparently be sum(x_izⁱ for in in range(n))ⁿ where the x_i are the roots of the polynomial and z would be a primitive n-th root of unity. For that we'd need to know the roots and the cyclic order for them. It's easy to compute the roots in Sage but they don't come with the cyclic structure. On the other hand, PARI computes the roots and gives the automorphisms in terms of these, but its roots are modulo p^e, so they contain more information than GF(p) but I don't know how to match them to the algebraic roots in Sage, or how to perform a computation on them and turn the result into a number field element. I'm still not sure how to do the recusrion once this is resolved, but at the moment, I'm really wondering whether any of you knows how to match these roots.

Alternative approaches to compute this resolvent from the coefficients of the polynomial look very ugly for order 4 and I know of no feasible approach for higher orders.

[8d15854a-f726-4f6b-88e7-82ec1970fbba]

comment:12

I guess I hadn't grasped the full significance of the term p-adic integer in the docs. The way I see it, we could use Zp(gal[1][0], prec=gal[1][1], type="fixed-mod") to represent the roots returned by PARI in the Sage world, perform arithmetic on them and then use the rational_reconstruction if some result is known to be rational. Will have to look more closely, but I fear that some other results will belong to some field extension instead, and I'm not sure whether we could reconstruct these as well. Will need more reading and thinking.

[slel]

comment:14

See a related discussion in this Sage devel thread. In particular, this message suggests one could go further than radical expressions, when radical expressions are not enough.

[slel]

Changed keywords from radical galois symbolic to radical, galois, symbolic