Computation of a minimal polynomial

Question

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)$?

Answer 1

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

Answer 2

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

Answer 3

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)

Answer 4

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.

Answer 5

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

Answer 6

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