# Computation of a minimal polynomial

It is relatively easy (but sometimes quite cumbersome) to compute the minimal polynomial of an algebraic number $\alpha$ when $\alpha$ is expressible in radicals. For example, the simple query

"minimal polynomial 2^(1/5)*(1-exp(2*pi*i/5))"

to Wolfram Alpha will compute the minimal polynomial of $\sqrt[5]{2}\left(1-\exp\left(\frac{2\pi i}{5}\right)\right)$. However, I do not know how to compute a minimal polynomial of an algebraic number that is not expressible in radicals. Moreover, I want to compute a minimal polynomial of linear combinations of algebraic numbers.

For example, take $f(x) = x^5 - x + 1$. This polynomial is irreducible, has Galois group $S_5$, which is not solvable, so by Abel-Ruffini Theorem the roots of $f(x)$ are not expressible in radicals. I want to take two distinct roots of $f(x)$, say $\alpha_1$ and $\alpha_2$, and compute the minimal polynomial $g(x)$ of $\alpha_1 - \alpha_2$. I know how to do this algebraically, because

$$g(x) = \prod\limits_{\substack{1 \leq i, j \leq 5,\\i\neq j}}\left(x - (\alpha_i - \alpha_j)\right),$$

but the computation just seems too painful. Is there a function in Sage, or Mathematica, or Maple, or PARI/GP, that allows to find $g(x)$?

• Sage has a field "QQbar" for exact algebraic numbers, you can take the roots of your original polynomials over that field (e.g., change ring and factor), then perform whatever computations you need, and use the ".minpoly" method to get the minimal polynomial. – Gro-Tsen Sep 3 '17 at 18:29
• math.stackexchange.com/a/155153/39797 is a good general point of view. – Dima Pasechnik Sep 5 '17 at 12:10

To compute the minimal polynomial of integer multiple of an algebraic integer is easy, so the only thing you need for linear combinations is the minimal polynomials of sums. Now, note that if $A$ is the companion matrix of $\alpha$ and $B$ is the companion matrix of $\beta,$ then $A\otimes I + I\otimes B$ a companion matrix of $\alpha + \beta.$ (Of course, $A\otimes B$ is a companion matrix of $\alpha \beta.$)

Mathematica code:

frobeniusCompanion[poly_, x_] /; PolynomialQ[poly, x] :=
Module[{n = Exponent[poly, x], coef},
coef = CoefficientList[poly, x]; coef = -Most[coef]/Last[coef];
SparseArray[{{1, j_} :> coef[[-j]], Band[{2, 1}] -> 1}, {n, n}]]

AA= frobeniusCompanion[x^5 -x -1, x]
CC = KroneckerProduct[AA, IdentityMatrix[5]] + KroneckerProduct[IdentityMatrix[5], AA]
cp = CharacteristicPolynomial[Normal[CC], x]
Factor[cp]

gives $$-\left(x^5-16 x-32\right) \left(x^{10}+3 x^6+11 x^5-4 x^2+4 x-1\right)^2$$

• Thank you, Igor! I think this is precisely what I was looking for. Can you give an example of how this works w.r.t. two distinct roots of x^5-x+1? I don't quite understand how this works when two distinct algebraic numbers have the same companion matrix. – Anton Sep 3 '17 at 16:26
• @Anton In that case, the characteristic polynomial will be reducible (because the roots of the characteristic polynomials will also include things like $2\alpha.$). I will add mathematica code. – Igor Rivin Sep 3 '17 at 16:40
• Thanks, I appreciate that! Though I don't quite understand what do you mean by "characteristic polynomial will be reducible". Say, I know that if a and b are two distinct roots of x^5-x+1, then the minimal polynomial of a-b is x^20-10*x^16-...+5000x^2+2869, a.k.a. the characteristic polynomial of a companion matrix of a-b, is irreducible. What's the problem here? – Anton Sep 3 '17 at 16:59
• @Anton The problem is that you cannot algebraically distinguish one from another. So, some differences $\alpha - \beta$ of roots of your favorite polynomial are equal to zero, so their minimal polynomial is $x.$ The characteristic polynomial of the rational number $\alpha - \beta$ will be thus divisible by powers of $x.$ – Igor Rivin Sep 3 '17 at 17:49
• Mathematica would also allow you to do things like MinimalPolynomial[Root[#^5 - # + 1 &, 1] - Root[#^5 - # + 1 &, 5], x], to use @Anton's example. – J. M. is not a mathematician Sep 3 '17 at 19:59

In Sage you can do something like:

x=polygen(QQbar)
f=x^5 - x +1
roots = f.roots(QQbar)
a1 = roots[1][0]
a2 = roots[-1][0]
(a1 - a2).minpoly()
x^20 - 10*x^16 - 95*x^12 + 625*x^10 - 40*x^8 + 3750*x^6 + 400*x^4 + 5000*x^2 + 2869

Using resultants can greatly help. For example, knowing that $\sqrt[5]{2}$ is a zero of $P(x)=x^5-2$ and $1-\exp\frac{2\pi i}{5}$ is a zero of $Q(x)=(x-1)^5+1$, we conclude $\sqrt[5]{2}\cdot (1-\exp\frac{2\pi i}{5})$ is a zero of $\mathrm{Res}_y(P(y),y^5Q(\frac{x}{y}))$. In PARI/GP:

? polresultant(y^5-2,y^5*((x/y-1)^5+1),y)
%1 = x^25 + 2500*x^15 + 50000*x^5
? factor(%)
%2 =
[                       x 5]

[x^20 + 2500*x^10 + 50000 1]

That is, we conclude that the minimal polynomial is $x^{20} + 2500x^{10} + 50000$. Numerical verification:

? subst(x^20 + 2500*x^10 + 50000, x, 2^(1/5)*(1-exp(2*Pi*I/5)) )
%3 = 5.380530270144733809 E-34 - 2.104034275277802617 E-34*I

(the default precision is 38 decimal digits)

Although exact calculations are to be preferred where possible, if all you have is a high-accuracy numeric value of a number which you know to be algebraic then you can use an integer relation finding algorithm (e.g. LLL or PSLQ) to find a relation between the first $n$ powers of your number. LLL, at least, is implemented in Sage, Mathematica, Maple, and PARI/GP.

• Mathematica, for instance, has the RootApproximant[] function. – J. M. is not a mathematician Sep 4 '17 at 13:30
• In PARI/GP, this can be done with algdep(). – Max Alekseyev Sep 4 '17 at 15:26

PARI / GP code :

p=x^5-x+1
py=subst(p,x,y)

K = nfinit(py)

L0 = rnfinit(K,nffactor(K,p)[2,1])
p2y=subst(L0.polabs,x,y)
L = nfinit(p2y)

f(w,i) = Mod(- subst(nffactor(L,w)[i,1],x,0), p2y)

a1=f(p,1)
a2=f(p,2)
minpoly(a1-a2)

gives :

x^20 - 10*x^16 - 95*x^12 + 625*x^10 - 40*x^8 + 3750*x^6 + 400*x^4 + 5000*x^2 + 2869

SymPy code :

from sympy import *

x,y = symbols('x y')
roots = solve(x**5 - x + 1)
print(minimal_polynomial(roots[0] - roots[1], y))

gives :

y**20 - 10*y**16 - 95*y**12 + 625*y**10 - 40*y**8 + 3750*y**6 + 400*y**4 + 5000*y**2 + 2869