Solving Solvable Polynomials with Galois Theory (Part 1)

06 Aug 2021 - Tags: sage , solving-solvable-polynomials

I’m super excited to be writing this post up! It’s been haunting me for almost exactly a month now, and it feels good to be close enough to done that I can finally share my hard work with the world ^_^.

I’ve been spending a lot of time studying galois theory in the past little while, since it’s easily my weakest area in the standard algebra curriculum. I’ve been meaning to really sit down and learn it for a few years now, but there’s always been something more pressing. I’m not a fan of the qual system overall, but I guess nothing is entirely without merit, and if I’m being honest with myself? I don’t know when I would have gotten around to learning galois theory if it weren’t for the upcoming quals.

Now, one of the most famous theorems in galois theory is that the roots of a polynomial $f$ (in $Q [x]$ , say) can be expressed with nested radicals if and only if the galois group $G$ of $f$ is solvable. So, if I were to give you a polynomial with solvable galois group, would you know how to actually… solve it?

The proof that’s given in most books is actually constructive, but enough details are left out to make it mostly unimplementable. I eventually found Gaal’s Classical Galois Theory With Examples, which is very explicit about how to do the solving process. In this blog post, we’re going to focus in on the case where most of the work is done: that of a cyclic group of prime order¹. You’ll remember that every solvable group is an interated extension of these groups (in fact, this is basically the definition of solvable), so by tackling this case, we’ll get the full solvability case by iteration. I’ll write up a follow up blog post soon where we talk about that process.

For now, though let’s see the algorithm:

Let $f$ be a monic polynomial of degree $n$ with cyclic galois group, generated by $σ$ .

Let $α_{0}, \dots, α_{n - 1}$ be the roots of $f$ – these are just symbols for now. But without loss of generality order them so that $σ α_{j} = α_{j + 1}$ .
Let $ω$ be an $n$ th root of unity
Look at the equations $θ_{k} = \sum_{0 \leq j < n} ω^{j k} α_{j}$ . Notice $σ θ_{k} = ω^{k} θ_{k}$ (do you see why?). This means

$σ θ_{k}^{n} = (σ θ_{k})^{n} = ω^{n k} θ_{k}^{n} = θ_{k}^{n}$

so $θ_{k}^{n}$ is fixed by $σ$ (and thus the whole galois group).
Since $θ_{k}^{n}$ is fixed by $G$ , it must lie in the base field. So it’s just a number, $ψ_{k}$ . That means $θ_{k} = \sqrt[n]{ψ_{k}}$ .
Now we recover $α_{0} = \frac{1}{n} \sum \sqrt[n]{ψ_{k}}$ . In fact, we can recover each of the $α_{j}$ as a weighted average of the $ψ_{k}$ . Do you see how?²

This is great and all, but there’s a bit of magic that happens in step $4$ … How do we actually figure out what the $ψ_{k}$ s are in the base field? Getting $ψ_{0}$ is pretty easy: It’s the first elementary symmetric polynomial, and as such, it’s just (the negation of) the $x^{n - 1}$ coefficient of $f$ !

This is actually a clue for how we might in theory get each of the $ψ_{k}$ . Even though the other $ψ_{k}$ aren’t symmetric in the $α_{j}$ , we can still leverage the theorem of symmetric polynomials (with some work) to figure out which numbers the $ψ_{k}$ really are. This is where Gaal’s book was first helpful – it introduced me to the following idea:

Consider the polynomial (which I now know is called the galois resolvent)

L (t) = \prod {t - \tilde{ψ} ∣ \tilde{ψ} \in S_{n} \cdot ψ_{1}} .

This product ranges over the orbit of $ψ_{1}$ under the action of the symmetric group (which permutes the $α_{j}$ inside $ψ_{1}$ ). Notice that each of the $ψ_{k}$ for $k \neq 0$ are in this orbit (do you see why?) so each nonzero $ψ_{k}$ is a root of $L$ .

A priori, $L$ lives in $Q (ω) [α_{0}, \dots, α_{n - 1}] [t]$ , since the $ψ$ s are all polynomials in the $α_{j}$ . But the coefficients of $L$ are all symmetric in the orbit of $ψ_{1}$ , and thus, symmetric in the $α_{j}$ (again, do you see why?). So $L$ is actually an element of $Q (ω) [t]$ ! But we know how to find the roots of a polynomial with constant coefficients (or rather, sage does), and these roots are exactly the $ψ_{k}$ !

I coded it up, and… it crashed my desktop. In lieu of downloading more ram, I tried to optimize my code (see here), which still didn’t work. I also tried to find a more effective procedure (see here), but it seems like this is really how it’s done. I know it’s been done before (in the gap radiroot package, for instance), but I didn’t want to have to reverse engineer someone else’s code unless I absolutely had to³. So, I kept looking.

Eventually I learned that $L$ is called the resolvent, and I spent some time learning more about resolvents and why they’re interesting and useful⁴. So my next step was to find an efficient algorithm for evaluating resolvents, and I found one! My code here is an implmentation of the algorithm described in Lehobey’s Resolvent Computations By Resultants Without Extraneous Powers, which is super readable if you’re interested!

It’s still kind of slow, but it’s less memory intensive, and it definitely gets the job done!

xxxxxxxxxx
 
def interpolating_functions(f):
  """
  Build a list of interpolating functions for f
​
  We require f : K[x1,...,xn][x]
  """
  R = f.parent().base_ring()
  xs = R.gens()
​
  n = f.degree()
​
  out = [f]  
  for k in range(f.degree()):
    fk = (out[-1] - out[-1].subs(x=xs[n-k-1]))/(x - xs[n-k-1])
    out += [f.parent()(fk.simplify_full())]
​
  return out[::-1]
​
​
def stabilizer(G,p):
  """
  Compute the stabilizer of p by G
  """
  elems = []
  for g in G:
    if p * g == p:
      elems += [g]
  return G.subgroup(elems)
​
​
​
def truncated_root(p,r,d):
  """
  Compute q so that q^r = p (working mod x^d)
​
  Assumes p : A[t] has constant term 1 
  and that such a q : A[t] actually exists!
  """
​
  r = int(r)
  t = p.variables()[0]
  n = p.degree()
  p = p.truncate(d)
​
  ps = p.coefficients(sparse=False)
​
  # these will be the coefficients of q
  qs = [1]
​
  for k in range(n // r):
    qs += [1/(k+1) * sum([(k+1 - (r+1)*j) * qs[j] * ps[k+1-j] / r for j in range(k+1)])]
​
  return sum([qs[j] * t^j for j in range(len(qs))])
​
​
​
def resolvent(f,Theta):
  """
  Compute mathcal{L}_{Theta,f} as per the paper
​
  Assumes f : K[x1,...,xn][x] and Theta : K[x1,...,xn]
  """
  R = Theta.parent()
  K = R.base_ring()
  xs = R.gens()
​
  SIterated.<t> = PolynomialRing(R)
  S = SIterated.flattening_morphism().codomain()
​
  T = K[t]
​
  Rj = (t - SIterated(Theta)).reverse()
  Rj = S(Rj)
​
  fs = interpolating_functions(f)
​
  HprevOrder = 1
​
  n = f.degree()
​
  for j in range(1,n+1):
    print(j, "/", n)
    Sj = SymmetricGroup(j)
    Hj = stabilizer(Sj,Theta)
    dj = factorial(j) / Hj.order()
    mj = Hj.order() / HprevOrder
​
    # update the previous order for the next cycle
    HprevOrder = Hj.order()
​
    # there's an annoying off-by-one error with the variable names
    # compared to everything else
    fj = S(fs[j].subs(x=xs[j-1]))
​
    res = fj.resultant(Rj, S(xs[j-1]))
​
    Rj = truncated_root(SIterated(res),mj,dj+1)
    Rj = S(Rj)
​
  return T(Rj).reverse()
​
​
def solveByRadicals(f):
  """
  Compute a root of f using radicals
​
  f(x) is assumed to be symbolic
  """
​
  n = int(f.degree(x))
  K.<w> = CyclotomicField(n)
​
  R = PolynomialRing(K,n,'x')
  xs = R.gens()
​
  R1 = R[x]
  f = R1(f)
​
  Theta = sum(xs[k] * w^(k) for k in range(n))
​
  # Theta^n is preserved under the action of the galois group,
  # while Theta itself is an eigenvector with eigenvalue w
  L = resolvent(f,Theta^n)
​
  psis = L.roots(multiplicities=False)
  thetas = [psi^(1/n) for psi in psis]
​
  # we need to choose the ~correct~ nth root for each psi.
  # I don't actually know how you're supposed to know which 
  # one is right, so we just try them all...
  #
  # There must be a better way to do this, but I want to start
  # working on other things.
​
  from itertools import product
  for es in product([w^k for k in range(n)], repeat=n-2):
    r = (-list(f)[-2] + thetas[0] + sum(es[k-1] * thetas[k] for k in range(1,n-1)))/n
    # there's definitely a better way to do this too...
    if abs(f(r).n()) < 0.000000001:
      return r
​
  # if we never found a root
  print("Uh oh!")
​
R = QQ[x]
​
deg3s = [x^3 - x^2 - a*x + b for (a,b) in [(26,-41), (32,79), (34,61), (36,4), (42,-80), (46,-103)]]
​
deg5s = [x^5 + x^4 - 4*x^3 - 3*x^2 + 3*x + 1, 
         x^5 + x^4 - 12*x^3 - 21*x^2 + 1*x + 5, 
         x^5 + x^4 - 16*x^3 + 5*x^2 + 21*x - 9, 
         x^5 + x^4 - 24*x^3 - 17*x^2 + 41*x - 13]
​
deg7s = [x^7 + x^6 - 12*x^5 - 7*x^4 + 28*x^3 + 14*x^2 - 9*x + 1]
​
fs = deg3s + deg5s + deg7s
​
@interact
def _ (f=selector(fs, label="$f$"), auto_update=False):
  show(solveByRadicals(f))

f

Help | Powered by SageMath

Accepted: {"data":{"application/sage-clear":{"changed":["f","auto_update"]},"text/plain":"Clear display"},"source":"sagecell"}

Accepted: {"data":{"application/sage-interact":{"new_interact_id":"f7e022c2-5961-4bad-944a-6cda5a6ae5af","controls":{"f":{"control_type":"selector","subtype":"list","values":11,"value_labels":["x^3 - x^2 - 26*x - 41","x^3 - x^2 - 32*x + 79","x^3 - x^2 - 34*x + 61","x^3 - x^2 - 36*x + 4","x^3 - x^2 - 42*x - 80","x^3 - x^2 - 46*x - 103","x^5 + x^4 - 4*x^3 - 3*x^2 + 3*x + 1","x^5 + x^4 - 12*x^3 - 21*x^2 + x + 5","x^5 + x^4 - 16*x^3 + 5*x^2 + 21*x - 9","x^5 + x^4 - 24*x^3 - 17*x^2 + 41*x - 13","x^7 + x^6 - 12*x^5 - 7*x^4 + 28*x^3 + 14*x^2 - 9*x + 1"],"default":0,"nrows":null,"ncols":null,"raw":true,"width":"None","label":"$f$","update":false},"auto_update":{"control_type":"button","width":"","text":"Update","raw":true,"label":"","update":true}},"layout":[[["f",1]],[["auto_update",1]],[["_output",1]]],"locations":null,"readonly":false},"text/plain":"Sage Interact"},"source":"sagecell"}

Messages

The sagemath online server times these out for polynomials of degree $5$ and $7$ , but you can run this code locally to see that it does give the right answer⁵. The degree $5$ case is slow, but manageable. The degree $7$ case is… really impressively slow. You’ve been warned.

As a slightly tricky exercise, we’re assuming throughout that we have access to roots of unity… But how do we know that we can write roots of unity in terms of radicals?

Show that roots of unity can always be written in terms of radicals⁶.

I’m fairly confident this will work on cyclic groups of composite order too, and I’ve even tested it on a few polynomials of degree $4$ . But the prime case is all we need for solvability, and I can say for sure that it always works in that case, so that’s what we’re going with. ↩
As a little hint, you’ll only need to scale the $ψ_{k}$ by powers of $ω$ . You should try writing these sums out, say in the $n = 3$ case, to try and get a handle for what’s going on. You might also want to watch Richard E Borcherds’ (characteristically excellent) video here

As a massive hint (and a cute connection with the rest of mathematics), the $θ_{k}$ are actually the Discrete Fourier Transform of the $α_{j}$ ! So once we know the $θ_{k}$ s, we can recover the $α_{j}$ s using the inverse DFT! ↩
Especially since the package is based on Andreas Distler’s thesis, and my German is…. nicht so gut. ↩
There’s not much sense writing up a post about it, though, since there’s actually a ton of great resources on this! Healey’s Resultants, Resolvents, and the Computation of Galois Groups (available here) certainly comes to mind. ↩
The answers are pretty unwieldy, though. So maybe it’s for the best you need to run it locally.

For instance, if $ζ$ is an $11$ th root of unity, then $Q (ζ)$ is a cyclic extension of degree $10$ . This means $ζ + ζ^{- 1}$ (which is of degree $2$ under $Q (ζ)$ ) generates an extension of degree $5$ over $Q$ .

Sage tells us the minimal polynomial of $ζ + ζ^{- 1}$ is $x^{5} + x^{4} - 4 x^{3} - 3 x^{2} + 3 x + 1$ , and so we can use this code to write $ζ + ζ^{- 1}$ as:
$\begin{array}{r} - \frac{α^{3} + 4 α^{2} + 16 α + 64}{320} {(- \frac{165}{64} α^{3} - \frac{385}{16} α^{2} - \frac{275}{4} α - 451)}^{\frac{1}{5}} \\ + \frac{α}{20} {(- \frac{55}{32} α^{3} + \frac{55}{8} α^{2} + \frac{275}{4} α - 176)}^{\frac{1}{5}} \\ + \frac{1}{10} {(\frac{385 α^{3} + 440 α^{2} + 3520 α - 4224}{2})}^{\frac{1}{5}} \\ + \frac{1}{5} {(- \frac{55}{32} α^{3} + \frac{165}{16} α^{2} - 55 α - 286)}^{\frac{1}{5}} \\ - \frac{1}{5} \end{array}$
where $α = \sqrt{5} + i \sqrt{2 \sqrt{5} + 10} - 1$ .

NB: I made the substitution for $α$ myself, so there might be a minor arithmetic error in the above expression… Though I don’t think that’s actually very important :P ↩
As a hint, let $ζ$ be an $n$ th root of unity for $n$ odd. First show that you can recover $ζ$ from $ζ + ζ^{- 1}$ , and that $ζ + ζ^{- 1}$ satisfies a (cyclic) equation of degree $\frac{n - 1}{2}$ . Then, inductively, we can solve this by radicals, and then get $ζ$ with one more square root.

Do you see how to do the $n$ even case as well? ↩