# How to efficiently compute the minimal polynomial of a number expressed in radicals?

Preamble:

I want to calculate the minimal polynomial of a number of the form

$$x=\sum_{i=1}^k \pm a_i^{1/k_i}$$

Where the $$a_i$$ are algebraic numbers also of this form with a finite expression and $$k_i$$ are positive integers. I can think of three approaches, all with serious problems:

1. Using the resultant repeatedly:

Apply the algorithm recursively until you arrive at a root of a rational number, therefore we can assume that we have the minimal polynomial of all $$a_i$$. Then compute the minimal polynomial of $$\pm a_1^{1/k_1}\pm a_2^{1/k_2}$$ using the resultant of their minimal polynomials then repeat for that with $$a_3, a_4...a_k$$. Let $$d_i$$ be the degree of the minimal polynomial of $$a_i^{1/k_i}$$. In general this method requires computing the determinant of square matrices of orders $$d_1+d_2$$, $$d_1+d_2+d_3,\ldots$$ up to $$d_1+d_2+\ldots+d_k$$. According to wikipedia calculating the determinant of a matrix of order $$n$$ has complexity of the order $$O(n^{2.373})$$ but practically it's faster to use the LU decomposition method (this is what, for example, the blaze c++ library does) which has a complexity of $$O(n^3)$$. If we have $$d_1=d_2=...=d$$ then the complexity is $$O((2d)^3+(3d)^3...(kd)^3)=O(d^3k^4)$$.

1. Using an Integer relation algorithm:

If we know that the degree of $$x$$ is $$n$$ then we can run for example the LLL algorithm on $${1,x,x^2...x^n}$$. Problems with that include:

• Determination of $$n$$: in the previous example we can get an upper bound of $$d^k$$ (yikes) so running the LLL algorithm on this would have complexity $$O(d^k)$$ assuming I'm interpreting wikipedia correctly. If $$n$$ does not match an upper bound then we either need to run the algorithm again or factor the polynomial.

• Precision: if we calculate these values outwards starting at the rational constants (for example $$\sqrt{2-\sqrt{2}}\approx\sqrt{0.5858}\approx 0.8367$$) we run into the problem that each operation has a tendency to exacerbate the errors of the previous ones. This is particularly problematic when subtracting two large numbers whose difference is small. So we'd need bounds on the error $$x^i\pm\epsilon$$ where the algorithm still returns the correct result. Edit: since writing this, I found this blog post by Jack Coughlin which details how to estimate this error. This greatly mitigates the problem.

Afterwards we'd need to check if the polynomial found really has $$x$$ as a root and whether it's irreducible (this is only necessary if we don't know $$n$$ in advance, then we'd have to use an upper bound).

1. Trying some symbolic manipulation, especially as a middle step of another algorithm:

In some cases isolating the term with the highest root might be helpful. For example

$$x=\sqrt{a}+\sqrt{b}$$

$$x^3-3\sqrt a x^2+3ax-a\sqrt a = b$$

$$(x^3+3ax-b)^2=9ax^4-6a^2x^2+a^3$$

I get the impression that this isn't efficient, so maybe we can use a lookup table for small exponents with little nesting.

Applications

Of course I want to implement this in a computer so it does all the number crunching for me. For example, recently I was trying to do a search over graphs composed of triangular lattices so what I did is set the vertices of the graph as points in $$\mathbb C$$ so I had to check whether several numbers of the form $$\eta(a+b\omega)+\mu$$ where $$w=\frac{1-\sqrt{-3}}2$$ are equal to each other with sympy but since it stores expressions "as is" there's no trivial way to check equality, so I needed to check if the minimal polynomial of $$\eta_1(a_1+b_1\omega)+\mu_1-\eta_2(a_2+b_2\omega)-\mu_2$$ was equal to the identity polynomial, which takes forever. Sympy uses the first algorithm by default and alternatively something involving Gröbner bases which I did not understand.

Of course the object-oriented-programming approach to solving this problem is declaring a class which represents an algebraically closed structure that holds all the numbers you're gonna need (and you better foresee all of them!). For example in c++:

// a member of a real or imaginary quadratic field
// a+b*sqrt(n) where a,b,n are integers
int a;
int b;
int n;
}
// ...
};


This will be fast to perform operations on, but if you have a number which is not of the form $$a+b\sqrt n$$ you need to rewrite everything from scratch, plus if $$n\equiv 1 \pmod 4$$, $$a$$ and $$b$$ need to be half integers, and if you wanted to perform operations on the field $$\mathbb Q[\sqrt 2]$$ instead the class should probably use a template type instead of ints...all of this complexity means that by the end the class will be over a hundred lines long.

If we could compute minimal polynomials fast though, we could have a class called algebraic which internally represents a number as its minimal polynomial together with a way of distinguishing it from the other roots. Then its usage could be something like:

algebraic a {"x^2-2", 0}; // the 0th root of x^2-2 i.e. sqrt(2)
algebraic b {"x^3-3", 0}; // in 3^(1/3)

// x = 2^(1/2)+3^(1/3)
// so x is an algebraic number of degree 6
// the internal representation of x holds the polynomial
// t^6-6t^4-6t^3+12t^2-36t+1
// as a list of 7 integers (1,0,-6,-6,12,-36,1) plus the number of the root
// which is an integer from 0 to 5
auto x = a + b;


Where the 0th, 1st, etc roots are ordered say by absolute value and then by complex argument.

This way we still get fast arithmetic without having to rewrite our class every time we need to work on a different algebraic structure.

Statement of the question:

What is the fastest algorithm to calculate the minimal polynomial of a number expressed in radicals? It should have good asymptotic complexity in the number of terms/the degree of the input. Performance for small entries is not so important because a look-up table can be implemented.

• You might want to take a look at this paper. Also Mathematica and Maple have a minimal polynomial function (which are likely heavily optimized), so it might be worth it to take a look at the algorithms those use. – Varun Vejalla Sep 4 '20 at 2:14
• @VarunVejalla, I prefer free software. But to be fair I don't use any computer algebra systems very often. If you know that mathematica or whatever does this fast and how it works and what algorithm it uses please post it as an answer. – Sophie Sep 4 '20 at 15:38
• Detalized version is ready. – Yuri Negometyanov Sep 11 '20 at 6:08
• New approach added – Yuri Negometyanov Sep 12 '20 at 11:37

$$\color{brown}{\textbf{The new version.}}$$

Each radical $$\sqrt[k]a$$ forms an array of the variables $$(1, p, p^2,\dots, p^{k-1}),$$ wherein there is not a problem to exclude them via a homogenius linear system.

From the expression $$\color{brown}{\mathbf{\ x=\sqrt[k] a}}\$$ should

$$\begin{cases} x-p_j = 0\\ p_j = \sqrt[k]a\, \omega_k^j\\ \omega_k = e^{^{\Large\frac{2\pi i}k}} \end{cases}\Rightarrow \begin{cases} \prod\limits_{j=0}^{m-1}(x-p_j) = 0\\ p_j = \sqrt[k]a\, \omega_k^j\\ j=0\dots k-1\\ \omega_k = e^{^{\Large\frac{2\pi i}k}} \end{cases}\Rightarrow x^k-a = 0.$$

From the expression $$\color{brown}{\mathbf{\ x=\sqrt a + \sqrtb}}\$$ should $$\begin{cases} x-q \pm \sqrt a = 0\\ q^3 = b \end{cases}\Rightarrow \begin{cases} (x-q+\sqrt a)(x-q+\sqrt a) = 0\\ q^3 = b \end{cases}\Rightarrow \begin{cases} (x-q_j)^2 - a = 0\\ q_j = \sqrtb\omega_3^j\\ \omega_3 = e^{^{\Large\frac{2\pi i}3}}, \end{cases}$$\begin{align} &\Big((x-\sqrt{b})^2-a\Big)\Big((x-\sqrt{b}\,\omega_3)^2-a\Big) \Big((x-\sqrt{b}\,\omega_3^2)^2-a\Big) = 0,\\[4pt] &\Big((x-\sqrt{b})(x-\sqrt{b}\,\omega_3)(x-\sqrt{b}\,\omega_3^2)\Big)^2\\[4pt] &-a\Big((x-\sqrt{b})^2(x-\sqrt{b}\,\omega_3)^2 +(x-\sqrt{b}\,\omega_3)^2(x-\sqrt{b}\,\omega_3^2)^2+ (x-\sqrt{b}\,\omega_3^2)^2 (x-\sqrt{b})^2\Big)\\[4pt] &+ a^2\Big((x-\sqrt{b})^2 + (x-\sqrt{b}\,\omega_3)^2 + (x-\sqrt{b}\,\omega_3^2)^2\Big) -a^3 =0,\\[4pt] &(x^3-b)^2-3ax(x^3+2b)+3a^2x^2-a^3 =0,\\[4pt] &x^6 - 3ax^4 - 2bx^3 + 3a^2x^2 - 6abx + b^2 - a^3 = 0 \end{align} (see also WA result).

From the expression $$\color{brown}{\mathbf{\ z=\sqrt a + \sqrtb + \sqrtc + \sqrtd,\ }}$$ should

$$\begin{cases} z-p-q-r-s = 0\\ p^2=a\\ q^3 = b\\ r^4 = c\\ s^5 = d \end{cases}\Rightarrow \begin{cases} x = p+q = z-r-s\\ f_2(x) = x^6 - 3ax^4 - 2bx^3 + 3a^2x^2 - 6abx + b^2 - a^3 = 0\\ r^4 = c\\ s^5 = d \end{cases}\Rightarrow \begin{cases} y = x+r = z-s\\ f_3(y) = \prod\limits_{j=0}^3 f_2(y-r_j) = 0\\ r_j = \sqrtc\, i^j\\ s^5 = d, \end{cases}$$

$$f_3(y)=$$ (see WA result).

Alternative approaches are considered below.

$$\color{brown}{\textbf{The old version.}}$$

Each radical $$\sqrt[k+1]a$$ forms an array of the variables $$(p, p^2,\dots, p^k),$$ wherein there is not a problem to exclude them via a homogenius linear system.

In particular, for the expression $$\color{brown}{\mathbf{\ x=\sqrt a + \sqrtb}}\$$ easily to get

$$\begin{cases} p+q-x = 0\\ p^2=a\\ q^3 = b \end{cases}\Rightarrow \begin{cases} p+q-x = 0\\ a + p(q-x) = 0\\ q^3 = b \end{cases}\Rightarrow \begin{cases} \begin{vmatrix} 1 & q-x \\ q-x & a \end{vmatrix} = 0\\ q^3 =b\end{cases}$$$$\begin{cases} q^2-2xq+x^2 - a = 0\\ b-2xq^2+(x^2-a)q = 0\\ bq - 2bx + (x^2-a) q^2 = 0 \end{cases}\Rightarrow \begin{vmatrix} 1 & - 2x & x^2 - a\\ -2x & x^2-a & b\\ x^2-a & b & -2bx \end{vmatrix} = 0,$$$$x^6 - 3ax^4 - 2bx^3 + 3a^2x^2 - 6abx + b^2 - a^3 = 0$$ (see also WA result).

"Brite force" way is $$\begin{cases} p+q-x = 0\\ a+pq-px = 0\\ pq+q^2-qx = 0\\ aq+pq^2-pqx = 0 \\ pq^2+b-q^2x = 0\\ aq^2+bp-pq^2x = 0 \end{cases}\Rightarrow \begin{pmatrix} -x & 1 & 1 & 0 & 0 & 0\\ a & -x & 0 & 1 & 0 & 0\\ 0 & 0 & -x & 1 & 1 & 0\\ 0 & 0 & a & -x & 0 & 1\\ b & 0 & 0 & 0 & -x & 1\\ 0 & b & 0 & 0 & a & -x\ \end{pmatrix} \begin{pmatrix} 1\\p\\q\\pq\\q^2\\pq^2 \end{pmatrix}= 0,$$$$\begin{vmatrix} -x & 1 & 1 & 0 & 0 & 0\\ a & -x & 0 & 1 & 0 & 0\\ 0 & 0 & -x & 1 & 1 & 0\\ 0 & 0 & a & -x & 0 & 1\\ b & 0 & 0 & 0 & -x & 1\\ 0 & b & 0 & 0 & a & -x\ \end{vmatrix} =\begin{vmatrix} 1 & 0 & 0 & 0 &-x & 1 \\ 0 & 1 & 0 & 0 & a &-x \\ -x & 1 & 1 & 0 & 0 & 0 \\ a &-x & 0 & 1 & 0 & 0 \\ 0 & 0 &-x & 1 & b & 0 \\ 0 & 0 & a &-x & 0 & b \end{vmatrix}$$$$=\begin{vmatrix} 1 & 0 & 0 & 0 &-x & 1 \\ 0 & 1 & 0 & 0 & a &-x \\ 0 & 1 & 1 & 0 &-x^2 & x \\ 0 &-x & 0 & 1 & ax &-a \\ 0 & 0 &-x & 1 & b & 0 \\ 0 & 0 & a &-x & 0 & b \end{vmatrix} =\begin{vmatrix} 1 & 0 & 0 & a &-x \\ 0 & 1 & 0 &-x^2-a & 2x \\ 0 & 0 & 1 & 2ax &-x^2-a \\ 0 &-x & 1 & b & 0 \\ 0 & a &-x & 0 & b \end{vmatrix}$$$$=\begin{vmatrix} 1 & 0 &-x^2-a & 2x \\ 0 & 1 & 2ax &-x^2-a \\ 0 & 1 & b-x(x^2+a) & -2x^2 \\ 0 &-x & a(x^2+a) & b-2ax \end{vmatrix} =\begin{vmatrix} 1 & 2ax &-x^2-a \\ 0 & b-x(x^2+3a) & -x^2+a \\ 0 & a(3x^2+a) & b-x(x^2+3a) \end{vmatrix}$$$$=\begin{vmatrix} b-x(x^2+3a) & -x^2+a \\ a(3x^2+a) & b-x(x^2+3a) \end{vmatrix}=0,$$$$x^6 - 3ax^4 - 2bx^3 + 3a^2x^2 - 6abx + b^2 - a^3 = 0.$$

If $$\color{brown}{\mathbf{\ x=\sqrt a + \sqrtb + \sqrtc + \sqrtd,\ }}$$ then

$$\begin{cases} p+q+r+s-x = 0\\ p^2=a\\ q^3 = b\\ r^4 = c\\ s^5 = d \end{cases}\Rightarrow \begin{cases} x_1 = x-q-r-s\\ p-x_1 = 0\\ a - px_1 = 0\\ q^3 = b\\ r^4 = c\\ s^5 = d \end{cases}\Rightarrow \begin{cases} \begin{vmatrix} 1 & -x_1 \\ -x_1 & a \end{vmatrix} = 0\\ q^3 = b\\ r^4 = c\\ s^5 = d, \end{cases}$$$$\begin{cases} x_2 = x - r - s = x_1+q\\ q^2-2x_2q+(x_2^2 - a) = 0\\ b-2x_2q^2+(x_2^2-a)q = 0\\ bq - 2bx_2 + (x_2^2-a) q^2 = 0\\ r^4 = c\\ s^5 = d \end{cases}\Rightarrow \begin{cases} f_2(x_2) = \begin{vmatrix} 1 & - 2x_2 & x_2^2 - a\\ -2x_2 & x_2^2-a & b\\ x_2^2-a & b & -2bx_2 \end{vmatrix} = 0\\ r^4 = c\\ s^5 = d, \end{cases}$$$$\begin{cases} x_3 = x - s = x_2+r\\ f_2(x_2) = x_2^6 - 3ax_2^4 - 2bx_2^3 + 3a^2x_2^2 - 6abx_2 + b^2 - a^3 = 0\\ f_2(x_3-r) = 0\\ r^4 = c\\ s^5 = d, \end{cases}$$$$\begin{cases} x_3 = x - s = x_2+r\\[4pt] r^6 + f_{25} r^5 + f_{24} r^4 + f_{23} r^3 + f_{22} r^2 + f_{21} r + f_{20} = 0\\[4pt] f_{25} = \dfrac1{5!}\dfrac{\partial^5}{\partial r^5}f_2(x_3-r)\bigg|_{r=0} = -6x_3\\[4pt] f_{24} = \dfrac1{4!}\dfrac{\partial^4}{\partial r^4}f_2(x_3-r)\bigg|_{r=0} = 15x_3^2-3a\\[4pt] f_{23} = \dfrac1{3!}\dfrac{\partial^3}{\partial r^3}f_2(x_3-r)\bigg|_{r=0} = -\dbinom63 x_3^3 + 3a\dbinom43 x_3 + 2b = -20x^3 + 12ax +2b\\[4pt] f_{22} = \dfrac1{2!}\dfrac{\partial^2}{\partial r^2}f_2(x_3-r)\bigg|_{r=0} = \dbinom62x_3^4 - 3a\dbinom42x_3^2 - 2b\dbinom32x_3 +3a^2\\[4pt] = 15x_3^4 - 18a x_3^2 - 6bx_3 + 3a^2\\[4pt] f_{21} = -\dfrac{\partial}{\partial r}f_2(x_3-r)\bigg|_{r=0} = - 6x_3^5 + 12ax_3^2 + 6bx_3^2 - 6a^2x_3 + 6ab\\[4pt] f_{20} = f_2(x_3) = x_3^6 - 3ax_3^4 - 2bx_3^3 + 3a^2x_3^2 - 6abx_3 + b^2 - a^3\\[4pt] r^4 = c\\[4pt] s^5 = d, \end{cases}$$$$\begin{cases} x_3 = x - s = x_2+r\\[4pt] f_{23} r^3 + (f_{22}+c) r^2 + (f_{21}+cf_{25}) r + (f_{20} + cf_{24}) = 0\\[4pt] cf_{23} + (f_{22}+c) r^3 + (f_{21}+cf_{25}) r^2 + (f_{20} + cf_{24}) r = 0\\[4pt] cf_{23} r + c(f_{22}+c) + (f_{21}+cf_{25}) r^3 + (f_{20} + cf_{24}) r^2 = 0\\[4pt] cf_{23} r^2 + c(f_{22}+c) r + c(f_{21}+cf_{25}) + (f_{20} + cf_{24}) r^3 = 0\\[4pt] s^5 = d, \end{cases}$$

and this leads to the determinant of a Toeplitz matrix $$f_3(x) = \begin{vmatrix} f_{20} + cf_{24} & f_{21} + cf_{25} & f_{22}+c &f_{23} \\ cf_{23} & f_{20} + cf_{24} & f_{21}+cf_{25} & f_{22}+c \\ cf_{22}+c^2 & cf_{23} & f_{20}+cf_{24} & f_{21}+cf_{25} \\ cf_{21}+c^2f_{25} & cf_{22}+c^2 & cf_{23} & f_{20}+cf_{24} \end{vmatrix} = 0,\tag1$$ where

$${\small\begin{cases} f_{25} = -6(x-s) \\[4pt] f_{24} = 15(x-s)^2 - 3a \\[4pt] f_{23} =-20(x-s)^3 + 12a(x-s) + 2b \\[4pt] f_{22} = 15(x-s)^4 - 18a(x-s)^2 - 6b(x-s) + 3a^2 \\[4pt] f_{21} = -6(x-s)^5 + 12a(x-s)^2 + 6b(x-s)^2 - 6a^2(x-s) + 6ab \\[4pt] f_{20} = (x-s)^6 - 3a(x-s)^4 - 2b(x-s)^3 + 3a^2(x-s)^2 - 6ab(x-s) + b^2 - a^3. \end{cases}}\tag2$$

Computation complexity of this approach can be defined as $$\sum\limits_{i} D(k_i),$$ where $$D(k)$$ is the computation complexity of the determinant $$k\times k.$$

It is less than the OP estimation.

Since the Groebner basis algorithm has a great computational complexity and the manipulations approach cannot be easily formalized, this approach looks the fastest of the proposed variants.

• I'm confused about your steps in the Gauss-Jordan elemination process (with $\det$ equations). Could you elaborate on the reason in each step? – ExercisingMathematician Sep 8 '20 at 22:13
• @AlgebraicGeometryStudent The full bases are $(1,p)$ in the first case and $(1,q,q^2)$ in the second case. Multiplying to $p$ or $q$ allows to get homogenius linear systems, which has non-zero solution only if the main determinant is zero. Therefore, proposed approach is not the Gauss-Jordan process. – Yuri Negometyanov Sep 8 '20 at 22:41
• Can you say a little more about the general case? – Integrand Sep 9 '20 at 19:59
• @Integrand Thanks foor advice, elaborated. – Yuri Negometyanov Sep 11 '20 at 6:06
• @AlgebraicGeometryStudent Detalized version. – Yuri Negometyanov Sep 11 '20 at 6:07